RICHAnalysis_macros/macros_analyze/bnl_rich_analyzer_pid.cpp

0001 /*
0002   Adjustments made to original:
0003   Fixed a |= somewhere that should have been !=
0004   Increased size of em[][] array
0005   Added if statement to make sure imom != 0
0006   Changed the i/o to append to existing output file
0007   Changed read() to read(string input_filename)
0008   Changed "tree_rich" to updated "eval_rich"
0009
0010   Personal preference: Changed the bounds on polar angle to [15,16] instead of [12,13]
0011 */
0012
0013 #include <Riostream.h>
0014 #include <stdio.h>
0015 #include <fstream>
0016 #include <sstream>
0017 #include <string>
0018
0019 #include <TMath.h>
0020 #include <TEnv.h>
0021 #include <TH2F.h>
0022 #include <TH1F.h>
0023 #include <TTree.h>
0024 #include <TDatabasePDG.h>
0025 #include <TStyle.h>
0026 #include <TLine.h>
0027 #include <TEllipse.h>
0028 #include <TCanvas.h>
0029 #include <TProfile.h>
0030 #include <TH1D.h>
0031 #include <TH3F.h>
0032 #include <TF1.h>
0033 #include <TText.h>
0034 #include <TLegend.h>
0035 #include <TChain.h>
0036 #include <TString.h>
0037 #include <TLorentzVector.h>
0038 #include <TRandom.h>
0039 #include <TMatrixTBase.h>
0040 #include <TMatrixT.h>
0041 #include "TPaveText.h"
0042 #include <TGraphErrors.h>
0043 #include <TFile.h>
0044 #include <TTree.h>
0045 #include <TSpline.h>
0046 #include <TFrame.h>
0047
0048 using namespace std;
0049
0050 class eic_bnl_rich {
0051
0052  private:
0053
0054  public:
0055
0056   Double_t ind_ray(double Ex, double Ey, double Ez, double Dx, double Dy, double Dz, double cx, double cy, double cz, double vx, double vy, double vz);
0057
0058   void acq(string input_filename, int ind, int pix_lab, int part);
0059   Double_t momentum, ch, ch_v[500];
0060   Double_t c_mean, c_rms, p_ang;
0061   Int_t sec;
0062
0063   double qe(double energy);
0064
0065   void read(string input_filename);
0066
0067 };
0068
0069 Double_t eic_bnl_rich::ind_ray(double Ex, double Ey, double Ez, double Dx, double Dy, double Dz, double cx, double cy, double cz, double vx, double vy, double vz){
0070
0071   //Double_t cx = 100.;
0072   //Double_t cy = 0.;
0073   //Double_t cz = 134.5;
0074
0075   //Double_t cx,cy,cz;
0076   //Double_t sx,sy,sz;
0077   LongDouble_t cex,cey,cez;
0078   LongDouble_t cdx,cdy,cdz;
0079
0080   Int_t i,iwhere;
0081
0082   LongDouble_t th,a,d;
0083      LongDouble_t x,dx;
0084
0085   LongDouble_t y,y1;
0086
0087   LongDouble_t eps = 0.00000000001;
0088   LongDouble_t R = 195.;  //rich
0089   //LongDouble_t R = 299.9;    //rich3
0090
0091   //eic_dual_rich *f = new eic_dual_rich();
0092
0093   Double_t sx,sy,sz;
0094   Double_t esx,esy,esz;
0095   Double_t es, theta_c;
0096
0097   cex = -cx+Ex;
0098   cey = -cy+Ey;
0099   cez = -cz+Ez;
0100
0101   cdx = -cx+Dx;
0102   cdy = -cy+Dy;
0103   cdz = -cz+Dz;
0104
0105   //cout<<"ce is: "<<cex<<"  "<<cey<<"  "<<cez<<endl;
0106   //cout<<"cd is: "<<cdx<<"  "<<cdy<<"  "<<cdz<<endl;
0107
0108   a = TMath::Sqrt(cex*cex+cey*cey+cez*cez);
0109   d = TMath::Sqrt(cdx*cdx+cdy*cdy+cdz*cdz);
0110   th = TMath::ACos((cdx*cex+cdy*cey+cdz*cez)/a/d);
0111
0112   //cout<<"a,d,th is: "<<a<<"  "<<d<<"  "<<th<<endl;
0113
0114   i = 0;
0115   x = th/2.;
0116   //cout<<"x, sinx, sin(th-x) is:  "<<x<<"  "<<sin(x)<<"  "<<sin(th-x)<<endl;
0117   y = R*(a*sin(x)-d*sin(th-x))+a*d*sin(th-2*x);
0118   y1 = R*(a*cos(x)+d*cos(th-x))-2*a*d*cos(th-2*x);
0119   //cout<<"y, y1 is:  "<<y<<"  "<<y1<<endl;
0120   //dx = -(f->newp(x, th, a, d, R)/f->newpp(x, th, a, d, R));
0121   dx = -y/y1;
0122
0123   while(TMath::Abs(dx)>eps && i<100){
0124
0125     x+=dx;
0126     y = R*(a*sin(x)-d*sin(th-x))+a*d*sin(th-2*x);
0127     y1 = R*(a*cos(x)+d*cos(th-x))-2*a*d*cos(th-2*x);
0128     //dx = -(f->newp(x, th, a, d, R)/f->newpp(x, th, a, d, R));
0129     dx = -y/y1;
0130     i++;
0131
0132   }
0133
0134   //if(i>=100) cout<<"Not convergent"<<endl;
0135
0136   if(i<100){
0137     sx = cx + (R*cos(x)/a-R*sin(x)/tan(th)/a)*cex + (R*sin(x)/d/sin(th))*cdx;
0138     sy = cy + (R*cos(x)/a-R*sin(x)/tan(th)/a)*cey + (R*sin(x)/d/sin(th))*cdy;
0139     sz = cz + (R*cos(x)/a-R*sin(x)/tan(th)/a)*cez + (R*sin(x)/d/sin(th))*cdz;
0140
0141   }
0142
0143   esx = sx - Ex;
0144   esy = sy - Ey;
0145   esz = sz - Ez;
0146
0147   //cout<<"S: "<<sx<<"  "<<sy<<"  "<<sz<<endl;
0148
0149   es = sqrt(esx*esx+esy*esy+esz*esz);
0150
0151   esx = esx/es;
0152   esy = esy/es;
0153   esz = esz/es;
0154
0155   theta_c = TMath::ACos((esx*vx+esy*vy+esz*vz));
0156
0157   return theta_c;
0158
0159 }
0160
0161 double eic_bnl_rich::qe(double energy){
0162
0163   Double_t qe = 0.;
0164
0165   qe = 12.5*(energy-6.2);
0166   if(energy>10.3) qe = 0.; // no photons below 120 nm = 10.3 eV
0167
0168   return qe;
0169
0170 }
0171
0172 void eic_bnl_rich::acq(string input_filename, int ind, int pix_lab, int part){
0173
0174   c_mean = 0.;
0175   c_rms = 0.;
0176   p_ang = 0.;
0177
0178   Double_t rtd = 180./TMath::Pi();
0179
0180   TFile *file=new TFile(input_filename.c_str());
0181   if (file->IsZombie()) {
0182     cout << "Error opening file" << input_filename << endl;
0183     exit(-1);
0184   }
0185   //else cout << "open file " << input_filename << endl;
0186
0187   TTree *eic_rich = (TTree*) file->Get("eval_rich");
0188
0189   Int_t eic_rich_event=0,eic_rich_bankid=0,eic_rich_volumeid=0,eic_rich_hitid=0,*eic_rich_pid=0,eic_rich_mpid=0,eic_rich_trackid=0,eic_rich_mtrackid=0,eic_rich_otrackid=0;
0190   Double_t eic_rich_hit_x=0,eic_rich_hit_y=0,eic_rich_hit_z=0, eic_rich_emi_x=0, eic_rich_emi_y=0, eic_rich_emi_z=0,eic_rich_px=0,eic_rich_py=0,eic_rich_pz=0,eic_rich_mpx=0,eic_rich_mpy=0,eic_rich_mpz=0,eic_rich_e=0,eic_rich_me=0,eic_rich_edep=0;
0191
0192   eic_rich->SetBranchAddress("event",&eic_rich_event);
0193   eic_rich->SetBranchAddress("hit_x",&eic_rich_hit_x);
0194   eic_rich->SetBranchAddress("hit_y",&eic_rich_hit_y);
0195   eic_rich->SetBranchAddress("hit_z",&eic_rich_hit_z);
0196   eic_rich->SetBranchAddress("emi_x",&eic_rich_emi_x);
0197   eic_rich->SetBranchAddress("emi_y",&eic_rich_emi_y);
0198   eic_rich->SetBranchAddress("emi_z",&eic_rich_emi_z);
0199   eic_rich->SetBranchAddress("px",&eic_rich_px);
0200   eic_rich->SetBranchAddress("py",&eic_rich_py);
0201   eic_rich->SetBranchAddress("pz",&eic_rich_pz);
0202   eic_rich->SetBranchAddress("mpx",&eic_rich_mpx);
0203   eic_rich->SetBranchAddress("mpy",&eic_rich_mpy);
0204   eic_rich->SetBranchAddress("mpz",&eic_rich_mpz);
0205   eic_rich->SetBranchAddress("e",&eic_rich_e);
0206   eic_rich->SetBranchAddress("me",&eic_rich_me);
0207   eic_rich->SetBranchAddress("edep",&eic_rich_edep);
0208   eic_rich->SetBranchAddress("bankid",&eic_rich_bankid);
0209   eic_rich->SetBranchAddress("volumeid",&eic_rich_volumeid);
0210   eic_rich->SetBranchAddress("hitid",&eic_rich_hitid);
0211   eic_rich->SetBranchAddress("pid",&eic_rich_pid);
0212   eic_rich->SetBranchAddress("mpid",&eic_rich_mpid);
0213   eic_rich->SetBranchAddress("trackid",&eic_rich_trackid);
0214   eic_rich->SetBranchAddress("mtrackid",&eic_rich_mtrackid);
0215   eic_rich->SetBranchAddress("otrackid",&eic_rich_otrackid);
0216
0217   TRandom rran;
0218   rran.SetSeed(0);
0219   Int_t ev_count[61];
0220   for(Int_t i=0;i<61;i++) ev_count[i]=0;
0221   Int_t ph_count = 0;
0222   Int_t ph_ch_count = 0;
0223
0224   TH2F *ph_det = new TH2F("ph_det","", 1666, -250., 250., 1666, -250., 250.);
0225   TH2F *ph_det1 = new TH2F("ph_det1","", 1666, -250., 250., 1666, -250., 250.);
0226   Double_t momvv_x[10000];
0227   Double_t momvv_y[10000];
0228   Double_t momvv_z[10000];
0229
0230   TH1F *h = new TH1F("","",1000,-70,70);
0231   TH1F *spectrum_ev = new TH1F("","",1000,4.,12.);
0232   TH1F *spectrum_nm = new TH1F("","",1000,100.,300.);
0233   TH1F *ch_ang = new TH1F("","",1000,0.,0.3);
0234   Double_t ang_cut = 0.;
0235   Double_t momv[3] = {0.,0.,0.};
0236   //Double_t momvv[10000][3];
0237   Double_t m_emi[3] = {0.,0.,0.};
0238   Double_t d_hit[3] = {0.,0.,0.};
0239   Double_t xbin = 0.;
0240   Double_t ybin = 0.;
0241   Double_t theta = 0.;
0242   Double_t phi = 0.;
0243   Double_t rr[3] = {0.,0.,0.};
0244   Double_t theta_cc = 0.;
0245   Double_t m[4]={0.000511, 0.13957018, 0.493677, 0.938272};
0246   Double_t em[3][10000];
0247   Double_t vv = 0.;
0248   Int_t imom = 0;
0249   Double_t mom_check = 999.;
0250   Int_t event_check = 1.;
0251   Int_t event_check1 = 1.;
0252   Int_t nph_count = 0;
0253
0254   //theta_cc = 0.02;
0255
0256     //cout<<theta_cc<<endl;
0257
0258   TH1F **ch_pi_h;
0259   ch_pi_h=new TH1F*[61];
0260   for(Int_t i=0;i<61;i++){
0261     ch_pi_h[i] = new TH1F(Form("ch_pi_h%d",i),"",10000,0,0.3);
0262   }
0263   TH1F *h_mom = new TH1F("","",100,0.,100.);
0264
0265
0266   for(Int_t i=0;(i<eic_rich->GetEntries());i++){
0267
0268     eic_rich->GetEntry(i);
0269
0270     if(eic_rich_event<=ind){
0271
0272       ph_count++;
0273
0274       em[0][ph_count-1] = eic_rich_emi_x;
0275       em[1][ph_count-1] = eic_rich_emi_y;
0276       em[2][ph_count-1] = eic_rich_emi_z;
0277
0278       //cout<<em[2][0]<<"  "<<em[2][11]<<endl;
0279
0280       vv = sqrt((em[0][11]-em[0][0])*(em[0][11]-em[0][0])+(em[1][11]-em[1][0])*(em[1][11]-em[1][0])+(em[2][11]-em[2][0])*(em[2][11]-em[2][0]));
0281
0282       if(ph_count>50){
0283     momvv_x[eic_rich_event-1]=(em[0][11]-em[0][0])/vv;
0284     momvv_y[eic_rich_event-1]=(em[1][11]-em[1][0])/vv;
0285     momvv_z[eic_rich_event-1]=(em[2][11]-em[2][0])/vv;
0286       }
0287       if(eic_rich_event!=event_check){
0288     ph_count=0;
0289     event_check = eic_rich_event;
0290     //cout<<event_check<<endl;
0291       }
0292     }
0293   }
0294
0295   for(Int_t i=0;(i<eic_rich->GetEntries());i++){
0296
0297     eic_rich->GetEntry(i);
0298     //if(eic_rich_event==1){
0299     //h->Fill(eic_rich_mpz);
0300     //}
0301
0302     ang_cut = TMath::ACos(eic_rich_mpz/sqrt(eic_rich_mpx*eic_rich_mpx+eic_rich_mpy*eic_rich_mpy+eic_rich_mpz*eic_rich_mpz))*rtd;
0303
0304
0305
0306     //if(eic_rich_event<=ind && rran.Uniform(0,100)<qe(eic_rich_edep*1000000000.)){
0307     if(eic_rich_event<=ind){
0308
0309       momentum = sqrt(eic_rich_mpx*eic_rich_mpx+eic_rich_mpy*eic_rich_mpy+eic_rich_mpz*eic_rich_mpz);
0310       imom = momentum;
0311       if (imom == 0)
0312     continue;
0313       theta_cc = acos(sqrt(imom*imom+m[part]*m[part])/imom/1.00054);
0314
0315       momv[0] = momvv_x[eic_rich_event-1];
0316       momv[1] = momvv_y[eic_rich_event-1];
0317       momv[2] = momvv_z[eic_rich_event-1];
0318
0319       //if(eic_rich_event==1)cout<<momv[2]<<"  "<<eic_rich_event<<endl;
0320
0321       //momv[0] = eic_rich_mpx/sqrt(eic_rich_mpx*eic_rich_mpx+eic_rich_mpy*eic_rich_mpy+eic_rich_mpz*eic_rich_mpz);
0322       //momv[1] = eic_rich_mpy/sqrt(eic_rich_mpx*eic_rich_mpx+eic_rich_mpy*eic_rich_mpy+eic_rich_mpz*eic_rich_mpz);
0323       //momv[2] = eic_rich_mpz/sqrt(eic_rich_mpx*eic_rich_mpx+eic_rich_mpy*eic_rich_mpy+eic_rich_mpz*eic_rich_mpz);
0324
0325       /*momv[0] = eic_rich_emi_x/sqrt(eic_rich_emi_x*eic_rich_emi_x+eic_rich_emi_y*eic_rich_emi_y+eic_rich_emi_z*eic_rich_emi_z);
0326       momv[1] = eic_rich_emi_y/sqrt(eic_rich_emi_x*eic_rich_emi_x+eic_rich_emi_y*eic_rich_emi_y+eic_rich_emi_z*eic_rich_emi_z);
0327       momv[2] = eic_rich_emi_z/sqrt(eic_rich_emi_x*eic_rich_emi_x+eic_rich_emi_y*eic_rich_emi_y+eic_rich_emi_z*eic_rich_emi_z);*/
0328
0329       //cout<<sqrt(eic_rich_mpx*eic_rich_mpx+eic_rich_mpy*eic_rich_mpy+eic_rich_mpz*eic_rich_mpz)<<endl;
0330
0331       theta = TMath::ACos(momv[2]);
0332       phi = TMath::ATan2(momv[1],momv[0]);
0333
0334       p_ang = theta*rtd;
0335
0336       //cout<<theta*rtd<<endl;
0337       //cout<<1240./(eic_rich_e*1000000000.)<<endl;
0338
0339       m_emi[0] = ((220.)/momv[2])*momv[0];  //mean emission point
0340       m_emi[1] = ((220.)/momv[2])*momv[1];
0341       m_emi[2] = ((220.)/momv[2])*momv[2];
0342
0343       /*m_emi[0] = eic_rich_emi_x; //real emission point
0344       m_emi[1] = eic_rich_emi_y;
0345       m_emi[2] = eic_rich_emi_z;*/
0346
0347       d_hit[0] = eic_rich_hit_x;
0348       d_hit[1] = eic_rich_hit_y;
0349       d_hit[2] = eic_rich_hit_z;
0350
0351       xbin = ph_det->GetXaxis()->FindBin(eic_rich_hit_x);
0352       ybin = ph_det->GetYaxis()->FindBin(eic_rich_hit_y);
0353
0354       if(pix_lab!=0)d_hit[0] = ph_det->GetXaxis()->GetBinCenter(xbin);
0355       if(pix_lab!=0)d_hit[1] = ph_det->GetYaxis()->GetBinCenter(ybin);
0356
0357       ph_det->Fill(d_hit[0],d_hit[1]);
0358
0359       //cout<<theta*rtd<<"  "<<phi*rtd<<endl;
0360       //cout<<"V: "<<momv[0]<<"  "<<momv[1]<<"  "<<momv[2]<<endl;
0361       //cout<<"E: "<<m_emi[0]<<"  "<<m_emi[1]<<"  "<<m_emi[2]<<endl;
0362       //cout<<"D: "<<d_hit[0]<<"  "<<d_hit[1]<<"  "<<d_hit[2]<<endl;
0363       //cout<<eic_rich_volumeid<<endl;
0364
0365       rr[0] = -18.5*TMath::Sin(eic_rich_volumeid*TMath::Pi()/4.);
0366       rr[1] = 18.5*TMath::Cos(eic_rich_volumeid*TMath::Pi()/4.);
0367       rr[2] = 75.;
0368       sec = eic_rich_volumeid;
0369
0370       //cout<<rr[0]<<"  "<<rr[1]<<"  "<<rr[2]<<endl;
0371
0372       ch = ind_ray(m_emi[0], m_emi[1], m_emi[2], d_hit[0], d_hit[1], d_hit[2], rr[0], rr[1], rr[2], momv[0], momv[1], momv[2]);
0373       //ch_v[ph_ch_count] = ch;
0374       ph_ch_count++;
0375       //cout<<ch<<endl;
0376       //if(ch<(theta_cc+0.01) && ch>(theta_cc-0.01))ch_ang->Fill(ch);
0377
0378       //if(ch<(theta_cc+0.008) && ch>(theta_cc-0.008) && theta*rtd>15. && theta*rtd<16. && (ph_det->GetBinContent(xbin,ybin) == 1)){ //here are the cuts: central polar angle [12,13] (changed to [15,16]) and the cuts on theta_cc are about 5 sigma (assuming the measured 1 p.e. resulution) interval around the central value of theta_cc, the last cut takes just 1 p.e. per pixel  ####################################################
0379     //if(rran.Uniform(0,100)<qe(eic_rich_edep*1000000000.))ch_pi_h[imom]->Fill(ch);
0380       if(ch<(theta_cc+0.008) && ch>(theta_cc-0.008) && theta*rtd>15. && theta*rtd<16.){
0381
0382     ch_pi_h[imom]->Fill(ch);
0383     spectrum_nm->Fill(1240./(eic_rich_edep*1000000000.));
0384     //nph_count++;
0385     ph_det1->Fill(d_hit[0],d_hit[1]);
0386
0387     //if((eic_rich_event-1)==0)cout<<nph_v[0]<<endl;
0388
0389     if(momentum!=mom_check){
0390       h_mom->Fill(imom);
0391       ev_count[imom]++;
0392       //if(imom==14) cout<<eic_rich_mpid<<endl;
0393       //if(imom==60)cout<<momentum<<"  "<<imom<<endl;
0394       //nph_count++;
0395       mom_check=momentum;
0396     }
0397       }
0398       //if(imom==14)cout<<ch<<"  "<<eic_rich_mpid<<"  "<<theta*rtd<<"  "<<phi*rtd<<endl;
0399       if(eic_rich_event!=event_check1){
0400     //if(eic_rich_event<100)cout<<event_check1<<endl;
0401     //cout<<ph_det_nosup->GetEntries()<<endl;
0402     ph_det->Reset();
0403     event_check1 = eic_rich_event;
0404     //nph_count=0;
0405     }
0406     }
0407
0408   }
0409
0410   //cout<<ph_count<<endl;
0411
0412   //ch_ang->Draw();
0413   //ch_pi_h[32]->Draw();
0414   //spectrum_nm->Draw();
0415   //cout<<ev_count[9]<<endl;
0416   //h_mom->Draw();
0417
0418   Int_t summm=0;
0419
0420   ofstream out("pid_test.txt", ios::app);
0421   streambuf *coutbuf = cout.rdbuf();
0422   cout.rdbuf(out.rdbuf());
0423
0424   for(Int_t i=0;i<61;i++){
0425     if(ch_pi_h[i]->GetEntries()!=0)cout<<ch_pi_h[i]->GetMean()<<"  "<<ch_pi_h[i]->GetRMS()<<"  "<<(int)(ch_pi_h[i]->GetEntries()/ev_count[i])<<"  "<<i<<"  "<<part<<endl;
0426     else if(i<60)cout<<ch_pi_h[i]->GetMean()<<"  "<<ch_pi_h[i]->GetRMS()<<"  "<<0<<"  "<<i<<"  "<<part<<endl;
0427     summm = summm + ev_count[i];
0428   }
0429
0430   //cout<<nph_count<<endl;
0431
0432   ph_det1->Draw("colz");
0433   ch_pi_h[45]->Draw();
0434
0435   //c_mean = ch_ang->GetMean();
0436   //c_rms = ch_ang->GetRMS();
0437
0438   //delete ch_ang;
0439 }
0440
0441
0442 void eic_bnl_rich::read(string input_filename){
0443
0444   ifstream inputFile;
0445   inputFile.open(input_filename);
0446   if (!inputFile) {
0447     cerr << "Error in opening the file";
0448     exit(1);
0449   }
0450
0451   Double_t theta_ch, sigma_theta, nph, momentum, ptype;
0452   Double_t ch_e_v[60],rms_e_v[60],nph_e_v[60],mom_e_v[60];
0453   Double_t ch_pi_v[60],rms_pi_v[60],nph_pi_v[60],mom_pi_v[60];
0454   Double_t ch_k_v[60],rms_k_v[60],nph_k_v[60],mom_k_v[60];
0455   Double_t ch_p_v[60],rms_p_v[60],nph_p_v[60],mom_p_v[60];
0456   Double_t n_sig_e_pi[60],n_sig_pi_k[60],n_sig_k_p[60],mm_e_pi[60],mm_pi_k[60],mm_k_p[60];
0457   Int_t idx[4] = {0,0,0,0};
0458
0459   while (inputFile>>theta_ch>>sigma_theta>>nph>>momentum>>ptype){
0460
0461     if(ptype==0){
0462       ch_e_v[idx[0]]=theta_ch;
0463       rms_e_v[idx[0]]=sigma_theta;
0464       nph_e_v[idx[0]]=nph*0.6;
0465       if(theta_ch==0) nph_e_v[idx[0]]=0;
0466       mom_e_v[idx[0]]=momentum;
0467       idx[0]++;
0468     }
0469     if(ptype==1){
0470       ch_pi_v[idx[1]]=theta_ch;
0471       rms_pi_v[idx[1]]=sigma_theta;
0472       nph_pi_v[idx[1]]=nph*0.6; //0.6 is a scale factor to normalize the MC value on expwerimantal data
0473       if(theta_ch==0) nph_pi_v[idx[1]]=0;
0474       mom_pi_v[idx[1]]=momentum;
0475       idx[1]++;
0476     }
0477     if(ptype==2){
0478       ch_k_v[idx[2]]=theta_ch;
0479       rms_k_v[idx[2]]=sigma_theta;
0480       nph_k_v[idx[2]]=nph*0.6;
0481       if(theta_ch==0) nph_k_v[idx[2]]=0;
0482       mom_k_v[idx[2]]=momentum;
0483       idx[2]++;
0484     }
0485     if(ptype==3){
0486       ch_p_v[idx[3]]=theta_ch;
0487       rms_p_v[idx[3]]=sigma_theta;
0488       nph_p_v[idx[3]]=nph*0.6;
0489       if(theta_ch==0) nph_p_v[idx[3]]=0;
0490       mom_p_v[idx[3]]=momentum;
0491       idx[3]++;
0492     }
0493     //cout<<theta_ch<<"  "<<sigma_theta<<"  "<<nph<<"  "<<momentum<<"  "<<ptype<<endl;
0494
0495   }
0496
0497   Int_t ii=0;
0498   Int_t iii=0;
0499   Int_t iiii=0;
0500
0501   for(Int_t i=0;i<idx[0];i++){
0502
0503     if(ch_e_v[i]!=0 && ch_pi_v[i]!=0 && i>2){
0504       n_sig_e_pi[ii] = (ch_e_v[i]-ch_pi_v[i])*sqrt((nph_e_v[i]+nph_pi_v[i])*0.5)/((rms_e_v[i]+rms_pi_v[i])*0.5);
0505       mm_e_pi[ii] = mom_e_v[i]+0.5;
0506       //cout<<"e/pi: "<<n_sig_e_pi[ii]<<"  "<<mm_e_pi[ii]<<endl;
0507       ii++;
0508     }
0509     if(ch_pi_v[i]!=0 && ch_k_v[i]!=0 && i>16){
0510       n_sig_pi_k[iii] = (ch_pi_v[i]-ch_k_v[i])*sqrt((nph_pi_v[i]+nph_k_v[i])*0.5)/((rms_pi_v[i]+rms_k_v[i])*0.5);
0511       mm_pi_k[iii] = mom_pi_v[i]+0.5;
0512       //cout<<ch_pi_v[i]<<"  "<<ch_k_v[i]<<endl;
0513       //cout<<"pi/k: "<<n_sig_pi_k[iii]<<"  "<<mm_pi_k[iii]<<endl;
0514       iii++;
0515     }
0516     if(ch_k_v[i]!=0 && ch_p_v[i]!=0 && i>30){
0517       n_sig_k_p[iiii] = (ch_k_v[i]-ch_p_v[i])*sqrt((nph_k_v[i]+nph_p_v[i])*0.5)/((rms_k_v[i]+rms_p_v[i])*0.5);
0518       mm_k_p[iiii] = mom_k_v[i]+0.5;
0519       //cout<<"k/p: "<<n_sig_k_p[iiii]<<"  "<<mm_k_p[iiii]<<endl;
0520       iiii++;
0521     }
0522
0523     //cout<<ch_e_v[i]<<"  "<<nph_e_v[i]<<endl;
0524
0525   }
0526
0527   TGraph *sig_pi = new TGraph(idx[1],mom_pi_v,rms_pi_v);
0528   TGraph *nph_pi = new TGraph(idx[1],mom_pi_v,nph_pi_v);
0529
0530   TGraph *gr_e_pi = new TGraph(ii,mm_e_pi,n_sig_e_pi);
0531   TGraph *gr_pi_k = new TGraph(iii,mm_pi_k,n_sig_pi_k);
0532   TGraph *gr_k_p = new TGraph(iiii,mm_k_p,n_sig_k_p);
0533
0534   TCanvas *c1 = new TCanvas("c1","",1000,600);
0535   c1->Divide(1,1);
0536   c1->cd(1);
0537   gr_e_pi->SetLineStyle(2);
0538   gr_e_pi->SetLineWidth(2);
0539   gr_e_pi->SetLineColor(kGreen);
0540   gr_e_pi->SetMarkerColor(kGreen);
0541   gr_e_pi->SetTitle("");
0542   gr_e_pi->GetXaxis()->SetTitle("Momentum [Gev/c]");
0543   gr_e_pi->GetYaxis()->SetTitle("N_{#sigma}^{Ring}");
0544   gr_e_pi->Draw("AC*");
0545   gr_pi_k->SetLineStyle(3);
0546   gr_pi_k->SetLineWidth(2);
0547   gr_pi_k->SetLineColor(kRed);
0548   gr_pi_k->SetMarkerColor(kRed);
0549   gr_pi_k->SetTitle("");
0550   gr_pi_k->Draw("C*");
0551   gr_k_p->SetLineStyle(4);
0552   gr_k_p->SetLineWidth(2);
0553   gr_pi_k->SetLineColor(kBlue);
0554   gr_pi_k->SetMarkerColor(kBlue);
0555   gr_k_p->SetTitle("");
0556   gr_k_p->Draw("C*");
0557   c1->Print("test1.pdf","pdf");
0558
0559   for(Int_t j=0;j<iiii;j++) cout<<n_sig_k_p[j]<<"  "<<mm_k_p[j]<<endl;
0560
0561   TCanvas *c2 = new TCanvas("c2","",1000,600);
0562   c2->Divide(1,1);
0563   c2->cd(1);
0564   sig_pi->SetTitle("");
0565   sig_pi->GetXaxis()->SetTitle("Momentum [Gev/c]");
0566   sig_pi->GetYaxis()->SetTitle("#sigma_{#theta_{C}}");
0567   sig_pi->Draw("AP*");
0568   c2->Print("test2.pdf","pdf");
0569
0570   TCanvas *c3 = new TCanvas("c3","",1000,600);
0571   c3->Divide(1,1);
0572   c3->cd(1);
0573   nph_pi->SetTitle("");
0574   nph_pi->GetXaxis()->SetTitle("Momentum [Gev/c]");
0575   nph_pi->GetYaxis()->SetTitle("N_{p.e.}");
0576   nph_pi->Draw("AP*");
0577   c3->Print("test3.pdf","pdf");
0578
0579 }