ROCsoft: beamtime/cern-oct12/go4/AnalysisMacros/spadic_PIDpurity

Go to the documentation of this file.
00001 #include "spadic_noise_Analysis.cxx"
00002 #include "TH3.h"
00003 
00004 
00005 
00006 Bool_t A(false), B(false), C(false), D(false), E(false), F(false), G(false), H (false), Rest(false);
00007 
00008 void conditionTest(Double_t Ch1, Double_t Ch2, Double_t Pb){
00009   A = false;
00010   B = false;
00011   C = false;
00012   D = false;
00013   E = false;
00014   F = false;
00015   G = false;
00016   H = false;
00017   Rest = false;
00018 
00019   Bool_t Ch1El = false;
00020   Bool_t Ch1Pi = false;
00021   Bool_t Ch2El = false;
00022   Bool_t Ch2Pi = false;
00023   Bool_t PbEl  = false;
00024   Bool_t PbPi  = false;
00025 
00026   if (Ch1 > pidcuts[0]  && Ch1 < pidcuts[1])  Ch1El = true;
00027   if (Ch2 > pidcuts[2]  && Ch2 < pidcuts[3])  Ch2El = true;
00028   if (Pb  > pidcuts[4]  && Pb  < pidcuts[5])  PbEl  = true;
00029   if (Ch1 > pidcuts[6]  && Ch1 < pidcuts[7])  Ch1Pi = true;
00030   if (Ch2 > pidcuts[8]  && Ch2 < pidcuts[9])  Ch2Pi = true;
00031   if (Pb  > pidcuts[10] && Pb  < pidcuts[11]) PbPi  = true;
00032 
00033   if      (Ch1Pi && Ch2Pi && PbPi)   A = true; // Pions 
00034   else if (Ch1Pi && Ch2Pi && PbEl)   B = true;
00035   else if (Ch1El && Ch2Pi && PbPi)   C = true;
00036   else if (Ch1El && Ch2El && PbEl)   D = true; // Electrons 
00037   else if (Ch1El && Ch2El && PbPi)   E = true;
00038   else if (Ch1Pi && Ch2El && PbEl)   F = true;
00039   else if (Ch1El && Ch2Pi && PbEl)   G = true;
00040   else if (Ch1Pi && Ch2El && PbPi)   H = true;
00041   else                            Rest = true;
00042   //else if (/*!Ch1Pi && !Ch2Pi && !PbPi && !Ch1El && !Ch2El && !PbEl*/!A && !B && !C && !D && !E && !F && !G && !H) Rest = true; 
00043 }
00044 
00045 void spadic_PIDpurity_Analysis(TString filename="data2012/TTrees/merged/3GeV_1775V_500V_B_K_3xEF700.root", Bool_t debug=false){
00046 
00047   gROOT->Reset(); 
00048   gROOT->SetStyle("Plain");
00049   gStyle->SetPalette(1,0);  
00050   gStyle->SetOptTitle(kFALSE);
00051   gStyle->SetPadTickX(1);                         //<-- tic marks on all axes
00052   gStyle->SetPadTickY(1);  
00053   gStyle->SetFillStyle(4000);
00054   gStyle->SetFrameFillStyle(0);
00055   if (!debug){
00056     gStyle->SetOptStat(kFALSE);
00057     gStyle->SetStatStyle(0);
00058   }
00059   gStyle->SetTitleStyle(0);
00060   gStyle->SetCanvasBorderSize(0);
00061   gStyle->SetFrameBorderSize(0);
00062   gStyle->SetLegendBorderSize(0);
00063   gStyle->SetStatBorderSize(0);
00064   gStyle->SetTitleBorderSize(0);
00065   const Int_t maxADC = 4000;  
00066 
00067   TString name, title;
00068   GetPidCuts(filename);
00069   //pidcuts[4] = 350;
00070 
00071   if (debug)
00072     printf("\nEl C1 min: %6i\nEl C1 max: %6i\nEl C2 min: %6i\nEl C2 max: %6i\nEl Pb min: %6i\nEl Pb max: %6i\nPi C1 min: %6i\nPi C1 max: %6i\nPi C2 min: %6i\nPi C2 max: %6i\nPi Pb min: %6i\nPi Pb max: %6i\n\n",
00073            pidcuts[0],pidcuts[1],pidcuts[2],pidcuts[3],pidcuts[4],pidcuts[5], // El
00074            pidcuts[6],pidcuts[7],pidcuts[8],pidcuts[9],pidcuts[10],pidcuts[11]// Pi
00075            );
00076 
00077   TFile *inputFile = new TFile(filename,"READ");
00078   if (!inputFile->IsOpen()){
00079     cout << "file not found: " << filename << endl;
00080   } else {
00081     cout << "file found: " << filename << endl;
00082     TString outpic = filename;
00083     outpic.ReplaceAll("data2012/TTrees/merged/","");
00084     outpic.ReplaceAll(".root","");
00085     TString outfilename;
00086     outfilename.Form("data2012/root/Purity/Spectra_%s.root",outpic.Data());
00087     //printf("  %s\n  %s\n  %s\n",filename.Data(), outpic.Data(), outfilename.Data());
00088     TFile *output = new TFile(outfilename,"RECREATE");
00089 
00090 
00091     TH3I* allPIDs = new TH3I("allPIDs","allPIDs",250,0,maxADC,250,0,4000,250,0,maxADC);
00092     allPIDs->SetXTitle("Cherenkov 1 [a.u.]");
00093     allPIDs->SetYTitle("Pb-glass calorimeter [a.u.]");
00094     allPIDs->SetZTitle("Cherenkov 2 [a.u.]");
00095 
00096     TH3I* APIDs = new TH3I("APIDs","A: Ch1Pi && Ch2Pi && PbPi",250,0,maxADC,250,0,4000,250,0,maxADC);
00097     APIDs->SetXTitle("Cherenkov 1 [a.u.]");
00098     APIDs->SetYTitle("Pb-glass calorimeter [a.u.]");
00099     APIDs->SetZTitle("Cherenkov 2 [a.u.]");
00100 
00101     TH3I* BPIDs = new TH3I("BPIDs","B: Ch1Pi && Ch2Pi && PbEl",250,0,maxADC,250,0,4000,250,0,maxADC);
00102     BPIDs->SetXTitle("Cherenkov 1 [a.u.]");
00103     BPIDs->SetYTitle("Pb-glass calorimeter [a.u.]");
00104     BPIDs->SetZTitle("Cherenkov 2 [a.u.]");
00105 
00106     TH3I* CPIDs = new TH3I("CPIDs","C: Ch1El && Ch2Pi && PbPi",250,0,maxADC,250,0,4000,250,0,maxADC);
00107     CPIDs->SetXTitle("Cherenkov 1 [a.u.]");
00108     CPIDs->SetYTitle("Pb-glass calorimeter [a.u.]");
00109     CPIDs->SetZTitle("Cherenkov 2 [a.u.]");
00110 
00111     TH3I* DPIDs = new TH3I("DPIDs","D: Ch1El && Ch2El && PbEl",250,0,maxADC,250,0,4000,250,0,maxADC);
00112     DPIDs->SetXTitle("Cherenkov 1 [a.u.]");
00113     DPIDs->SetYTitle("Pb-glass calorimeter [a.u.]");
00114     DPIDs->SetZTitle("Cherenkov 2 [a.u.]");
00115 
00116     TH3I* EPIDs = new TH3I("EPIDs","E: Ch1El && Ch2El && PbPi",250,0,maxADC,250,0,4000,250,0,maxADC);
00117     EPIDs->SetXTitle("Cherenkov 1 [a.u.]");
00118     EPIDs->SetYTitle("Pb-glass calorimeter [a.u.]");
00119     EPIDs->SetZTitle("Cherenkov 2 [a.u.]");
00120 
00121     TH3I* FPIDs = new TH3I("FPIDs","F: Ch1Pi && Ch2El && PbEl",250,0,maxADC,250,0,4000,250,0,maxADC);
00122     FPIDs->SetXTitle("Cherenkov 1 [a.u.]");
00123     FPIDs->SetYTitle("Pb-glass calorimeter [a.u.]");
00124     FPIDs->SetZTitle("Cherenkov 2 [a.u.]");
00125 
00126     TH3I* GPIDs = new TH3I("GPIDs","G: Ch1El && Ch2Pi && PbEl",250,0,maxADC,250,0,4000,250,0,maxADC);
00127     GPIDs->SetXTitle("Cherenkov 1 [a.u.]");
00128     GPIDs->SetYTitle("Pb-glass calorimeter [a.u.]");
00129     GPIDs->SetZTitle("Cherenkov 2 [a.u.]");
00130 
00131     TH3I* HPIDs = new TH3I("HPIDs","H: Ch1Pi && Ch2El && PbPi",250,0,maxADC,250,0,4000,250,0,maxADC);
00132     HPIDs->SetXTitle("Cherenkov 1 [a.u.]");
00133     HPIDs->SetYTitle("Pb-glass calorimeter [a.u.]");
00134     HPIDs->SetZTitle("Cherenkov 2 [a.u.]");
00135 
00136     TH3I* RestPIDs = new TH3I("RestPIDs","RestPIDs",250,0,maxADC,250,0,4000,250,0,maxADC);
00137     RestPIDs->SetXTitle("Cherenkov 1 [a.u.]");
00138     RestPIDs->SetYTitle("Pb-glass calorimeter [a.u.]");
00139     RestPIDs->SetZTitle("Cherenkov 2 [a.u.]");
00140 
00141     //TH2I *allC1C2 = new TH2I();
00142 
00143     TH1D* allC1 = new TH1D("allC1","allC1",500,0,4000);
00144     allC1->SetXTitle("Cherenkov 1 [a.u.]");
00145     allC1->SetLineColor(8);
00146     TH1D* allC2 = new TH1D("allC2","allC2",500,0,4000);
00147     allC2->SetXTitle("Cherenkov 2 [a.u.]");
00148     allC2->SetLineColor(8);
00149     TH1D* allPb = new TH1D("allPb","allPb",500,0,4000);
00150     allPb->SetXTitle("Pb-glass calorimeter [a.u.]");
00151     allPb->SetLineColor(8);
00152     TH1D* restC1 = new TH1D("restC1","restC1",500,0,4000);
00153     restC1->SetXTitle("Cherenkov 1 [a.u.]");
00154     restC1->SetLineColor(16);
00155     TH1D* restC2 = new TH1D("restC2","restC2",500,0,4000);
00156     restC2->SetXTitle("Cherenkov 2 [a.u.]");
00157     restC2->SetLineColor(16);
00158     TH1D* restPb = new TH1D("restPb","restPb",500,0,4000);
00159     restPb->SetXTitle("Pb-glass calorimeter [a.u.]");
00160     restPb->SetLineColor(16);
00161     TH1D* C1_e = new TH1D("C1_e","C1_e",500,0,4000);
00162     C1_e->SetXTitle("Cherenkov 1 [a.u.]");
00163     C1_e->SetLineColor(9);
00164     TH1D* C1_p = new TH1D("C1_p","C1_p",500,0,4000);
00165     C1_p->SetXTitle("Cherenkov 1 [a.u.]");
00166     TH1D* C1_m = new TH1D("C1_m","C1_m",500,0,4000);
00167     C1_m->SetXTitle("Cherenkov 1 [a.u.]");
00168     C1_m->SetLineColor(7);
00169     TH1D* C2_e = new TH1D("C2_e","C2_e",500,0,4000);
00170     C2_e->SetXTitle("Cherenkov 2 [a.u.]");
00171     C2_e->SetLineColor(9);
00172     TH1D* C2_p = new TH1D("C2_p","C2_p",500,0,4000); 
00173     C2_p->SetXTitle("Cherenkov 2 [a.u.]");
00174     TH1D* C2_m = new TH1D("C2_m","C2_m",500,0,4000); 
00175     C2_m->SetXTitle("Cherenkov 2 [a.u.]");
00176     C2_m->SetLineColor(7);
00177     TH1D* Pb_e = new TH1D("Pb_e","Pb_e",500,0,4000);
00178     Pb_e->SetXTitle("Pb-glass calorimeter [a.u.]");
00179     Pb_e->SetLineColor(9);
00180     TH1D* Pb_p = new TH1D("Pb_p","Pb_p",500,0,4000);
00181     Pb_p->SetXTitle("Pb-glass calorimeter [a.u.]");
00182     TH1D* Pb_m = new TH1D("Pb_m","Pb_m",500,0,4000);
00183     Pb_m->SetXTitle("Pb-glass calorimeter [a.u.]");
00184     Pb_m->SetLineColor(7);
00185     TStopwatch timer;
00186     timer.Start();
00187 
00188     TTree* theTree=NULL;
00189     TKey* kee=NULL;
00190     TIter iter(inputFile->GetListOfKeys());
00191     while ( ( kee=dynamic_cast<TKey*>(iter()) ) !=0 ) {
00192       theTree = dynamic_cast<TTree*>(kee->ReadObj());
00193       if (theTree)
00194         break; // we take first Tree in file, disregarding its name...
00195     }
00196     if(theTree==NULL) {
00197       cout <<"Error: no Tree in file "<< filename.Data() << endl;
00198       return;
00199     }
00200     TCernOct12UnpackEvent* evnt = new TCernOct12UnpackEvent;
00201     TGo4EventElement* theBase=evnt;
00202     evnt->synchronizeWithTree(theTree, &theBase);
00203     const Int_t entries = theTree->GetEntries();
00204     printf("%i Events found in TTree\n",entries);
00205     TMbsCrateEvent* fCrateInputEvent = NULL;
00206     Double_t Ch1(0), Ch2(0), Pb(0);
00207     for(Int_t i=0;i<entries;++i) {
00208       Statusbar(i,entries);
00209       theTree->GetEntry(i);
00210       
00211       fCrateInputEvent=dynamic_cast<TMbsCrateEvent*>(evnt->GetSubEvent("MBSCRATE")); 
00212       if (fCrateInputEvent->fIsPulser) continue;
00213      
00214       Ch1 = fCrateInputEvent->fData1182[1][0];
00215       Ch2 = fCrateInputEvent->fData1182[0][0];
00216       Pb  = fCrateInputEvent->fData1182[1][1];
00217       if (debug) printf("%8i Ch1: %6f    Ch2: %6f    Pb: %6f\n",i,Ch1,Ch2,Pb);  
00218       getPID(Ch1, Ch2, Pb);
00219 
00220       allPIDs->Fill(Ch1, Pb, Ch2);
00221       allC1->Fill(Ch1);
00222       allC2->Fill(Ch2);
00223       allPb->Fill(Pb);
00224       conditionTest(Ch1, Ch2, Pb);
00225       if (isMyon){
00226         C1_m->Fill(Ch1);
00227         C2_m->Fill(Ch2);
00228         Pb_m->Fill(Pb);
00229       }
00230       if (A){ //pion
00231         APIDs->Fill(Ch1, Pb, Ch2);
00232         C1_p->Fill(Ch1);
00233         C2_p->Fill(Ch2);
00234         Pb_p->Fill(Pb);
00235       }
00236       if (B){
00237         BPIDs->Fill(Ch1, Pb, Ch2);
00238       }
00239       if (C){
00240         CPIDs->Fill(Ch1, Pb, Ch2);
00241       }
00242       if (D){ //electron
00243         DPIDs->Fill(Ch1, Pb, Ch2);
00244         C1_e->Fill(Ch1);
00245         C2_e->Fill(Ch2);
00246         Pb_e->Fill(Pb);
00247       }
00248       if (E){
00249         EPIDs->Fill(Ch1, Pb, Ch2);
00250       }
00251       if (F){
00252         FPIDs->Fill(Ch1, Pb, Ch2);
00253       }
00254       if (G){
00255         GPIDs->Fill(Ch1, Pb, Ch2);
00256       }
00257       if (H){
00258         HPIDs->Fill(Ch1, Pb, Ch2);
00259       }
00260       if (Rest && !isMyon){
00261         RestPIDs->Fill(Ch1, Pb, Ch2); 
00262         restC1->Fill(Ch1);
00263         restC2->Fill(Ch2);
00264         restPb->Fill(Pb);
00265       }
00266     }
00267     inputFile->Close();
00268     timer.Stop();
00269     Double_t rtime = timer.RealTime();
00270     Double_t ctime = timer.CpuTime();
00271 
00272     printf("\n\n*********************************************************\n");
00273     printf("    RealTime = %i min. %.1f s,  CpuTime = %i min. %.1f s\n",Int_t(rtime/60),(rtime/60 - Int_t(rtime/60))*60.,Int_t(ctime/60),(ctime/60 - Int_t(ctime/60))*60.);
00274     printf("             %.1f events/s,             %.1f events/s\n",Float_t(entries)/rtime,Float_t(entries)/ctime);
00275     printf("*********************************************************\n\n");
00276     Int_t all(allPIDs->GetEntries()),
00277       a(APIDs->GetEntries()),
00278       b(BPIDs->GetEntries()),
00279       c(CPIDs->GetEntries()),
00280       d(DPIDs->GetEntries()),
00281       e(EPIDs->GetEntries()),
00282       f(FPIDs->GetEntries()),
00283       g(GPIDs->GetEntries()),
00284       h(HPIDs->GetEntries()),
00285       rest(RestPIDs->GetEntries());
00286 
00287     printf(" all: %8i (%6.2f%%)\n   A: %8i (%6.2f%%) poins\n   B: %8i (%6.2f%%)\n   C: %8i (%6.2f%%)\n   D: %8i (%6.2f%%) electrons\n   E: %8i (%6.2f%%)\n   F: %8i (%6.2f%%)\n   G: %8i (%6.2f%%)\n   H: %8i (%6.2f%%)\nrest: %8i (%6.2f%%)\n sum: %8i (%6.2f%%)\n",
00288            all,100.*all/all,a,100.*a/all,b,100.*b/all,c,100.*c/all,d,100.*d/all,e,100.*e/all,f,100.*f/all,g,100.*g/all,h,100.*h/all,rest,100.*rest/all,a+b+c+d+e+f+g+h+rest,100.*(a+b+c+d+e+f+g+h+rest)/all);
00289 
00290     if (a == 0) a = 1;
00291     if (b == 0) b = 1;
00292     if (c == 0) c = 1;
00293     if (d == 0) d = 1;
00294     if (e == 0) e = 1;
00295     if (f == 0) f = 1;
00296     if (g == 0) g = 1;
00297     if (h == 0) h = 1;
00298 
00299     Double_t epsilon_Cal_e(0), epsilon_Cal_pi(0), epsilon_Ch1_e(0),  epsilon_Ch1_pi(0), epsilon_Ch2_e(0),  epsilon_Ch2_pi(0), epsilon_e(0), epsilon_pi(0);
00300     epsilon_Cal_e = 1./(1. + (d / e));
00301     epsilon_Ch1_e = 1./(1. + (d / f));
00302     epsilon_Ch2_e = 1./(1. + (d / g));
00303     epsilon_e = epsilon_Cal_e * epsilon_Ch1_e * epsilon_Ch2_e;
00304     epsilon_Cal_pi = 1./(1. + (a / b));
00305     epsilon_Ch1_pi = 1./(1. + (a / c));
00306     epsilon_Ch2_pi = 1./(1. + (a / h));
00307     epsilon_pi = epsilon_Cal_pi * epsilon_Ch1_pi * epsilon_Ch2_pi;
00308 
00309     printf("\n    Pb            Ch1           Ch2           epsilon\ne:  %E  %E  %E  %E  \npi: %E  %E  %E  %E  \n\n",epsilon_Cal_e,epsilon_Ch1_e,epsilon_Ch2_e,epsilon_e,epsilon_Cal_pi,epsilon_Ch1_pi,epsilon_Ch2_pi,epsilon_pi);
00310   
00311     if (output->IsOpen()){
00312       if (debug)
00313         {
00314           TCanvas *cC = new TCanvas("cC","cC",1.8*800,1.8*600);
00315           cC->Divide(3,3);
00316           cC->cd(1);
00317           //RestPIDs->Draw("surf");
00318           allPIDs->Draw("LEGO");
00319           cC->cd(2);
00320           APIDs->Draw("LEGO");
00321           cC->cd(3);
00322           BPIDs->Draw("LEGO");
00323           cC->cd(4);
00324           CPIDs->Draw("LEGO");
00325           cC->cd(5);
00326           DPIDs->Draw("LEGO");
00327           cC->cd(6);
00328           EPIDs->Draw("LEGO");
00329           cC->cd(7);
00330           FPIDs->Draw("LEGO");
00331           cC->cd(8);
00332           GPIDs->Draw("LEGO");
00333           cC->cd(9);
00334           HPIDs->Draw("LEGO");
00335           output->cd();  
00336           cC->Write("", TObject::kOverwrite);
00337         } 
00338 
00339       Double_t allEvents = allC1->Integral();
00340       printf(" all: %8i (%6.2f%%)\n",(Int_t)allEvents,100.*allEvents/all);
00341       allC1->Scale(1./allEvents);//allC1->Integral());
00342       allC2->Scale(1./allEvents);//allC2->Integral());
00343       allPb->Scale(1./allEvents);//allPb->Integral());
00344       restC1->Scale(1./allEvents);//restC1->Integral());
00345       restC2->Scale(1./allEvents);//restC2->Integral());
00346       restPb->Scale(1./allEvents);//restPb->Integral());
00347       C1_e->Scale(1./allEvents);//C1_e->Integral());
00348       C1_p->Scale(1./allEvents);//C1_p->Integral());
00349       C1_m->Scale(1./allEvents);//C1_p->Integral());
00350       C2_e->Scale(1./allEvents);//C2_e->Integral());
00351       C2_p->Scale(1./allEvents);//C2_p->Integral());
00352       C2_m->Scale(1./allEvents);//C1_p->Integral());
00353       Pb_e->Scale(1./allEvents);//Pb_e->Integral());
00354       Pb_p->Scale(1./allEvents);//Pb_p->Integral());
00355       Pb_m->Scale(1./allEvents);//Pb_p->Integral());
00356 
00357       output->cd(); 
00358       //Double_t max_e(C1_e->GetBinContent(C1_e->GetMaximumBin())), max_p(C1_p->GetBinContent(C1_p->GetMaximumBin()));
00359       TF1 *gC1_e = new TF1("gC1_e","gaus",pidcuts[0],C1_e->GetBinCenter(C1_e->GetMaximumBin()));
00360       gC1_e->SetLineWidth(1);
00361       gC1_e->SetLineColor(2);
00362       TF1 *gC2_e = new TF1("gC2_e","gaus",pidcuts[2],C2_e->GetBinCenter(C2_e->GetMaximumBin()));
00363       gC2_e->SetLineWidth(1);
00364       gC2_e->SetLineColor(2);
00365       TF1 *gPb1_e = new TF1("gPb1_e","gaus",pidcuts[4],Pb_e->GetBinCenter(Pb_e->GetMaximumBin()));
00366       gPb1_e->SetLineWidth(1);
00367       gPb1_e->SetLineColor(2);
00368       TF1 *gPb2_e = new TF1("gPb2_e","gaus",Pb_e->GetBinCenter(Pb_e->GetMaximumBin()),pidcuts[5]);
00369       gPb2_e->SetLineWidth(1);
00370       gPb2_e->SetLineColor(2);
00371       TF1 *gC1_p = new TF1("gC1_p","gaus",pidcuts[6],pidcuts[7]);
00372       gC1_p->SetLineWidth(1);
00373       gC1_p->SetLineColor(2);
00374       TF1 *gC2_p = new TF1("gC2_p","gaus",pidcuts[8],pidcuts[9]);//C2_p->GetBinCenter(C2_p->GetMaximumBin()));//
00375       gC2_p->SetLineWidth(1);
00376       gC2_p->SetLineColor(2);
00377       TF1 *gPb_p = new TF1("gPb_p","gaus(0)+gaus(3)+gaus(6)",pidcuts[10],pidcuts[11]);//Pb_e->GetBinCenter(Pb_e->GetMaximumBin()));//pidcuts[11]);//2700);
00378       gPb_p->SetLineWidth(1);
00379       gPb_p->SetLineColor(2);
00380       gPb_p->SetParameter(0, 0.006);//3496);
00381       gPb_p->SetParameter(1, 699);//606);
00382       gPb_p->SetParameter(2, 146);//154);
00383       gPb_p->SetParameter(3, 0.035);//16270);
00384       gPb_p->SetParameter(4, 684);//694);
00385       gPb_p->SetParameter(5, 37);//39);
00386       gPb_p->SetParameter(6, 0.0018);//172);
00387       gPb_p->SetParameter(7, 373);//821);
00388       gPb_p->SetParameter(8, 795);//454);
00389 
00390 
00391       TF1 *g_C1 = new TF1("g_C1","gaus(0)+gaus(3)+gaus(6)",320,1.25*(C1_e->GetBinCenter(C1_e->GetMaximumBin())));
00392       g_C1->SetLineWidth(1);
00393       g_C1->SetLineColor(6);
00394 
00395       g_C1->SetParameter(0, 1.11777e-01);// pion
00396       g_C1->SetParLimits(0, 1e-5, 1);
00397       g_C1->SetParameter(1, C1_p->GetBinCenter(C1_p->GetMaximumBin()));//3.54611e+02);
00398       g_C1->SetParameter(2, 7.56330e+00);
00399 
00400       g_C1->SetParameter(3, 2.09163e-03);// myon
00401       g_C1->SetParLimits(3, 1e-5, 1);
00402       g_C1->SetParameter(4, C1_m->GetBinCenter(C1_m->GetMaximumBin()));//6.08804e+02);
00403       //g_C1->SetParLimits(4, 1.25*(C1_p->GetBinCenter(C1_p->GetMaximumBin())), 0.70*(C1_e->GetBinCenter(C1_e->GetMaximumBin())));      
00404       g_C1->SetParameter(5, 1.07958e+02);
00405       g_C1->SetParLimits(5, 10, 300);
00406 
00407       g_C1->SetParameter(6, 1.130092-02);// electron
00408       g_C1->SetParLimits(6, 1e-5, 1);
00409       g_C1->SetParameter(7, C1_e->GetBinCenter(C1_e->GetMaximumBin()));//5.85493e+02);
00410       g_C1->SetParLimits(7, 0.95*C1_e->GetBinCenter(C1_e->GetMaximumBin()) , 1.05*C1_e->GetBinCenter(C1_e->GetMaximumBin()));
00411       g_C1->SetParameter(8, 1.73518e+02);
00412       g_C1->SetParLimits(8, 0,4000);
00413 
00414       TF1 *g_C2 = new TF1("g_C2","gaus(0)+gaus(3)+gaus(6)",300,1.25*(C2_e->GetBinCenter(C2_e->GetMaximumBin())));
00415       g_C2->SetLineWidth(1);
00416       g_C2->SetLineColor(6);
00417 
00418       g_C2->SetParameter(0, 1.26312e-01);// pion
00419       g_C2->SetParLimits(0, 1e-5, 1);
00420       g_C2->SetParameter(1, C2_p->GetBinCenter(C2_p->GetMaximumBin()));//3.28660e+02);
00421       g_C2->SetParameter(2, 8.06925e+00);
00422 
00423       g_C2->SetParameter(3, 1.30409e-03);// myon
00424       g_C2->SetParLimits(3, 1e-5, 1);
00425       g_C2->SetParameter(4, C2_m->GetBinCenter(C2_m->GetMaximumBin()));
00426       //g_C2->SetParLimits(4, 1.10*(C2_p->GetBinCenter(C2_p->GetMaximumBin())), 0.60*(C2_e->GetBinCenter(C2_e->GetMaximumBin())));      
00427       g_C2->SetParameter(5, 2.41815e+01);
00428       g_C2->SetParLimits(5, 10,50);
00429 
00430       g_C2->SetParameter(6, 8.375564-03);// electron
00431       g_C2->SetParLimits(6, 1e-5, 1);
00432       g_C2->SetParameter(7, C2_e->GetBinCenter(C2_e->GetMaximumBin()));//4.16179e+02);
00433       g_C2->SetParLimits(7, 0.95*C2_e->GetBinCenter(C2_e->GetMaximumBin()) , 1.05*C2_e->GetBinCenter(C2_e->GetMaximumBin()));
00434       g_C2->SetParameter(8, 2.10801e+02);
00435       g_C2->SetParLimits(8, 0,4000);
00436 
00437       if (debug) printf("C1 myon [%6.2f, %6.2f]  [%6.2f, %6.2f]\n",C1_p->GetBinCenter(C1_p->GetMaximumBin()), C1_e->GetBinCenter(C1_e->GetMaximumBin()),
00438                         1.25*(C1_p->GetBinCenter(C1_p->GetMaximumBin())), 0.70*(C1_e->GetBinCenter(C1_e->GetMaximumBin())));
00439       if (debug) printf("C2 myon [%6.2f, %6.2f]  [%6.2f, %6.2f]\n",C2_p->GetBinCenter(C2_p->GetMaximumBin()), C2_e->GetBinCenter(C2_e->GetMaximumBin()),
00440                         1.10*(C2_p->GetBinCenter(C2_p->GetMaximumBin())), 0.60*(C2_e->GetBinCenter(C2_e->GetMaximumBin())));
00441 
00442       TF1 *g_draw = new TF1("g_draw","gaus",0,4000);
00443       g_draw->SetLineWidth(1);
00444       g_draw->SetLineColor(kOrange);
00445       /*
00446         3GeV/c
00447         C1 fit on myons and electrons
00448 
00449         EXT PARAMETER                                   STEP         FIRST   
00450         NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
00451 
00452         1  const        1.18555e-02   1.89777e-02   1.98348e-05  -1.06042e-02 // myon
00453         2  mean         1.35578e+03   2.02666e+02   3.64636e-04   3.94911e-04
00454         3  sigma        1.73336e+02   1.81602e+02   6.69301e-02   7.75157e-07
00455         4  const        1.88695e-03   1.04824e-02   2.32852e-05  -4.16375e-03 // electron
00456         5  mean         5.80714e+02   1.39558e+02   3.35215e-03  -4.59047e-05
00457         6  sigma        1.42798e+02   7.53364e+02   2.50860e-01  -4.43909e-07
00458 
00459         1  const        1.11519e-01   8.21435e-01   1.06118e-03  -1.65024e-05 // pion
00460         2  mean         3.54611e+02   1.55122e+01   2.14058e-04   2.30851e-04
00461         3  sigma        7.51765e+00   1.06722e+01   5.02180e-04  -6.53331e-05
00462         4  const        1.18917e-02   1.88549e-02   6.72532e-05   1.67643e-03 // electron
00463         5  mean         1.35514e+03   2.02753e+02   7.95256e-04  -2.20698e-04
00464         6  sigma        1.72854e+02   2.44815e+01   3.87382e-03   4.04079e-05
00465         7  const        1.92235e-03   1.02686e-02   7.88927e-05  -1.50732e-03 // myon
00466         8  mean         5.85493e+02   8.19587e+01   5.31307e-03  -6.07572e-06
00467         9  sigma        1.37099e+02   2.09178e+01   1.69806e-02  -1.71586e-05
00468         __________________________________________________________________________________
00469         C2 fit on myons and electrons
00470 
00471         EXT PARAMETER                                   STEP         FIRST   
00472         NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
00473 
00474         1  const        8.30548e-03   1.73549e-02   6.86730e-06  -1.11522e-03 // myon
00475         2  mean         9.22404e+02   1.05339e+03   1.95778e-01  -3.81617e-07
00476         3  sigma        1.93678e+02   7.25385e+02   1.29769e-01   4.10184e-07
00477         4  const        1.34654e-03   2.19041e-02   7.39886e-06  -3.13219e-03 // electron
00478         5  mean         4.17995e+02   3.70401e+02   1.70182e-01  -4.37387e-07
00479         6  sigma        2.59027e+01   4.71255e+02   1.44508e-01   2.04406e-07
00480 
00481         1  const        1.26312e-01   2.42052e-01   3.31324e-04   2.13792e-04 // pion
00482         2  mean         3.28660e+02   1.38632e+01   2.86049e-04   6.41990e-05
00483         3  sigma        8.06925e+00   1.23596e+00   8.06849e-03   7.93309e-06
00484         4  const        8.34995e-03   1.16594e-02   6.22369e-05   2.09342e-03 // electron
00485         5  mean         9.60960e+02   1.07833e+02   2.12445e-03  -3.90537e-05
00486         6  sigma        2.20000e+02   4.27189e+01   2.97556e-02   8.92212e-07
00487         7  const        1.25924e-03   6.77964e-02   2.23469e-04  -3.68309e-04 // myon
00488         8  mean         4.16179e+02   5.72496e+01   4.15726e-03   1.47764e-05
00489         9  sigma        2.34271e+01   3.11161e+00   1.71880e-01   1.83821e-06
00490       */
00491       TLegend *leg = new TLegend(0.65,0.85-(0.04*7),0.85,0.85);
00492       leg->SetFillStyle(0);
00493       leg->SetTextSize(0.035);
00494       leg->SetLineStyle(0);
00495       leg->SetLineColor(0);
00496       leg->AddEntry(allC1, "full spectrum","l");
00497       leg->AddEntry(C1_e,  "e spectrum","l");
00498       leg->AddEntry(C1_p,  "#pi spectrum","l");
00499       leg->AddEntry(C1_m,  "#mu spectrum","l");
00500       leg->AddEntry(restC1,"rest spectrum","l");
00501       leg->AddEntry(g_C1,  "3 gaussian fit","l");
00502       //leg->AddEntry(g_draw,"3 single gaus","l");
00503       TCanvas * c2 = new TCanvas("c2","c2",2*800,1*600);
00504       c2->Divide(2,1);
00505       c2->cd(1)->SetLogy(1);
00506       allC1->DrawCopy();
00507       allC1->Fit(g_C1,"0QR");
00508       allC1->Fit(g_C1,"0QR");  
00509       C1_p->DrawCopy("same");
00510       //C1_p->Fit(gC1_p,"0QR");
00511       //gC1_p->DrawCopy("same");
00512       C1_e->DrawCopy("same");
00513       //C1_e->Fit(gC1_e,"0QR");
00514       //gC1_e->DrawCopy("same");
00515       C1_m->DrawCopy("same");
00516       restC1->DrawCopy("same"); 
00517       allC1->DrawCopy("same");
00518       allC1->DrawCopy("AXIS,same");
00519       g_C1->DrawCopy("same");
00520       /*
00521         g_draw->SetParameter(0, g_C1->GetParameter(0)); // pion
00522         g_draw->SetParameter(1, g_C1->GetParameter(1));
00523         g_draw->SetParameter(2, g_C1->GetParameter(2));
00524         g_draw->DrawCopy("same");
00525         g_draw->SetParameter(0, g_C1->GetParameter(3)); // myon
00526         g_draw->SetParameter(1, g_C1->GetParameter(4));
00527         g_draw->SetParameter(2, g_C1->GetParameter(5));
00528         g_draw->DrawCopy("same");
00529         g_draw->SetParameter(0, g_C1->GetParameter(6)); // electron
00530         g_draw->SetParameter(1, g_C1->GetParameter(7));
00531         g_draw->SetParameter(2, g_C1->GetParameter(8));
00532         g_draw->DrawCopy("same");
00533       */
00534       leg->Draw("same");
00535       c2->cd(2)->SetLogy(1);
00536       allC2->DrawCopy();
00537       allC2->Fit(g_C2,"0QR");
00538       allC2->Fit(g_C2,"0QR");
00539       C2_p->DrawCopy("same");
00540       //C2_p->Fit(gC2_p,"0QR");
00541       //gC2_p->DrawCopy("same");
00542       C2_e->DrawCopy("same");
00543       //C2_e->Fit(gC2_e,"0QR");
00544       //gC2_e->DrawCopy("same");
00545       C2_m->DrawCopy("same");
00546       restC2->DrawCopy("same");
00547       allC2->DrawCopy("same");
00548       allC2->DrawCopy("AXIS,same");
00549       g_C2->DrawCopy("same");
00550       /*
00551         g_draw->SetParameter(0, g_C2->GetParameter(0)); // pion
00552         g_draw->SetParameter(1, g_C2->GetParameter(1));
00553         g_draw->SetParameter(2, g_C2->GetParameter(2));
00554         g_draw->DrawCopy("same");
00555         g_draw->SetParameter(0, g_C2->GetParameter(3)); // myon
00556         g_draw->SetParameter(1, g_C2->GetParameter(4));
00557         g_draw->SetParameter(2, g_C2->GetParameter(5));
00558         g_draw->DrawCopy("same");
00559         g_draw->SetParameter(0, g_C2->GetParameter(6)); // electron
00560         g_draw->SetParameter(1, g_C2->GetParameter(7));
00561         g_draw->SetParameter(2, g_C2->GetParameter(8));
00562         g_draw->DrawCopy("same");
00563       */
00564       leg->Draw("same");
00565       c2->Update();
00566       outpic.ReplaceAll(".",",");
00567       c2->SaveAs("data2012/pics/Purity/Cherenkovs_"+outpic+".pdf");
00568       c2->SaveAs("data2012/pics/Purity/Cherenkovs_"+outpic+".png");
00569 
00570       TLegend *leg_pb = new TLegend(0.65,0.85-(0.04*6),0.85,0.85);
00571       leg_pb->SetFillStyle(0);
00572       leg_pb->SetTextSize(0.035);
00573       leg_pb->SetLineStyle(0);
00574       leg_pb->SetLineColor(0);
00575       leg_pb->AddEntry(allPb, "full spectrum","l");
00576       leg_pb->AddEntry(Pb_e,  "e spectrum","l");
00577       leg_pb->AddEntry(Pb_p,  "#pi spectrum","l");
00578       leg_pb->AddEntry(Pb_m,  "#mu spectrum","l");
00579       leg_pb->AddEntry(restPb,"rest spectrum","l");
00580       TCanvas * c3 = new TCanvas("c3","c3",1*800,1*600);
00581       c3->cd()->SetLogy(1);
00582       allPb->DrawCopy();
00583       Pb_p->DrawCopy("same");
00584       //Pb_p->Fit(gPb_p,"0QR");
00585       //gPb_p->DrawCopy("same");
00586       Pb_e->DrawCopy("same");    
00587       //Pb_e->Fit(gPb1_e,"0QR");
00588       //gPb1_e->DrawCopy("same");
00589       //Pb_e->Fit(gPb2_e,"0QR");
00590       //gPb2_e->DrawCopy("same");
00591       Pb_m->DrawCopy("same");
00592       restPb->DrawCopy("same");
00593       allPb->DrawCopy("AXIS,same");
00594       leg_pb->Draw("same");
00595       c3->Update();
00596       outpic.ReplaceAll(".",",");
00597       c3->SaveAs("data2012/pics/Purity/Calorimeter_"+outpic+".pdf");
00598       c3->SaveAs("data2012/pics/Purity/Calorimeter_"+outpic+".png");
00599 
00600       TCanvas * c1 = new TCanvas("c1","c1",2*800,2*600);
00601       c1->Divide(2,2);
00602       c1->cd(1)->SetLogy(1);
00603       allC1->DrawCopy();
00604       allC1->Fit(g_C1,"0QR");
00605       g_C1->DrawCopy("same");
00606       C1_p->DrawCopy("same");
00607       C1_p->Fit(gC1_p,"0QR");
00608       gC1_p->DrawCopy("same");
00609       C1_e->DrawCopy("same");
00610       C1_e->Fit(gC1_e,"0QR");
00611       gC1_e->DrawCopy("same");
00612       C1_m->DrawCopy("same");
00613       restC1->DrawCopy("same");
00614       printf("\n const          mean           sigma    \n%.4E     %.4E     %.4E pion\n%.4E     %.4E     %.4E myon\n%.4E     %.4E     %.4E electron\nChi^2/NDF %E\n",
00615              g_C1->GetParameter(0),
00616              g_C1->GetParameter(1),
00617              g_C1->GetParameter(2),
00618              g_C1->GetParameter(3),
00619              g_C1->GetParameter(4),
00620              g_C1->GetParameter(5),
00621              g_C1->GetParameter(6),
00622              g_C1->GetParameter(7),
00623              g_C1->GetParameter(8),
00624              Double_t(g_C1->GetChisquare())/Double_t(g_C1->GetNDF()));
00625       g_draw->SetParameter(0, g_C1->GetParameter(0)); // pion
00626       g_draw->SetParameter(1, g_C1->GetParameter(1));
00627       g_draw->SetParameter(2, g_C1->GetParameter(2));
00628       g_draw->DrawCopy("same");
00629       g_draw->SetParameter(0, g_C1->GetParameter(3)); // myon
00630       g_draw->SetParameter(1, g_C1->GetParameter(4));
00631       g_draw->SetParameter(2, g_C1->GetParameter(5));
00632       g_draw->DrawCopy("same");
00633       g_draw->SetParameter(0, g_C1->GetParameter(6)); // electron
00634       g_draw->SetParameter(1, g_C1->GetParameter(7));
00635       g_draw->SetParameter(2, g_C1->GetParameter(8));
00636       g_draw->DrawCopy("same");
00637       c1->cd(2)->SetLogy(1);
00638       allC2->DrawCopy();
00639       allC2->Fit(g_C2,"0QR");
00640       g_C2->DrawCopy("same");
00641       C2_p->DrawCopy("same");
00642       C2_p->Fit(gC2_p,"0QR");
00643       gC2_p->DrawCopy("same");
00644       C2_e->DrawCopy("same");
00645       C2_e->Fit(gC2_e,"0QR");
00646       gC2_e->DrawCopy("same");
00647       C2_m->DrawCopy("same");
00648       restC2->DrawCopy("same");
00649       printf("\n%.4E     %.4E     %.4E pion\n%.4E     %.4E     %.4E myon\n%.4E     %.4E     %.4E electron\nChi^2/NDF %E\n",
00650              g_C2->GetParameter(0),
00651              g_C2->GetParameter(1),
00652              g_C2->GetParameter(2),
00653              g_C2->GetParameter(3),
00654              g_C2->GetParameter(4),
00655              g_C2->GetParameter(5),
00656              g_C2->GetParameter(6),
00657              g_C2->GetParameter(7),
00658              g_C2->GetParameter(8),
00659              Double_t(g_C2->GetChisquare())/Double_t(g_C2->GetNDF()));
00660       g_draw->SetParameter(0, g_C2->GetParameter(0)); // pion
00661       g_draw->SetParameter(1, g_C2->GetParameter(1));
00662       g_draw->SetParameter(2, g_C2->GetParameter(2));
00663       g_draw->DrawCopy("same");
00664       g_draw->SetParameter(0, g_C2->GetParameter(3)); // myon
00665       g_draw->SetParameter(1, g_C2->GetParameter(4));
00666       g_draw->SetParameter(2, g_C2->GetParameter(5));
00667       g_draw->DrawCopy("same");
00668       g_draw->SetParameter(0, g_C2->GetParameter(6)); // electron
00669       g_draw->SetParameter(1, g_C2->GetParameter(7));
00670       g_draw->SetParameter(2, g_C2->GetParameter(8));
00671       g_draw->DrawCopy("same");
00672       c1->cd(3)->SetLogy(1);
00673       allPb->DrawCopy();
00674       Pb_p->DrawCopy("same");
00675       Pb_p->Fit(gPb_p,"0QR");
00676       gPb_p->DrawCopy("same");
00677       Pb_e->DrawCopy("same");    
00678       Pb_e->Fit(gPb1_e,"0QR");
00679       gPb1_e->DrawCopy("same");
00680       Pb_e->Fit(gPb2_e,"0QR");
00681       gPb2_e->DrawCopy("same");
00682       Pb_m->DrawCopy("same");
00683       restPb->DrawCopy("same");
00684 
00685       printf("\n$\\pi$ & %8.2f & %8.2f & %8.2f &           & %8.2f & %8.2f & %8.2f \\\\ \n$\\mu$ & %8.2f & %8.2f & %8.2f &           & %8.2f & %8.2f & %8.2f \\\\ \ne     & %8.2f & %8.2f & %8.2f &           & %8.2f & %8.2f & %8.2f \\\\ \n\n",
00686              fabs(g_C1->GetParameter(1) - g_C1->GetParameter(1)) / g_C1->GetParameter(2),//C1 pi
00687              fabs(g_C1->GetParameter(4) - g_C1->GetParameter(1)) / g_C1->GetParameter(5),
00688              fabs(g_C1->GetParameter(7) - g_C1->GetParameter(1)) / g_C1->GetParameter(8),
00689              fabs(g_C2->GetParameter(1) - g_C2->GetParameter(1)) / g_C2->GetParameter(2),//C2
00690              fabs(g_C2->GetParameter(4) - g_C2->GetParameter(1)) / g_C2->GetParameter(5),
00691              fabs(g_C2->GetParameter(7) - g_C2->GetParameter(1)) / g_C2->GetParameter(8),
00692              fabs(g_C1->GetParameter(1) - g_C1->GetParameter(4)) / g_C1->GetParameter(2),//C1 my
00693              fabs(g_C1->GetParameter(4) - g_C1->GetParameter(4)) / g_C1->GetParameter(5),
00694              fabs(g_C1->GetParameter(7) - g_C1->GetParameter(4)) / g_C1->GetParameter(8),
00695              fabs(g_C2->GetParameter(1) - g_C2->GetParameter(4)) / g_C2->GetParameter(2),//C2
00696              fabs(g_C2->GetParameter(4) - g_C2->GetParameter(4)) / g_C2->GetParameter(5),
00697              fabs(g_C2->GetParameter(7) - g_C2->GetParameter(4)) / g_C2->GetParameter(8),
00698              fabs(g_C1->GetParameter(1) - g_C1->GetParameter(7)) / g_C1->GetParameter(2),//C1 el
00699              fabs(g_C1->GetParameter(4) - g_C1->GetParameter(7)) / g_C1->GetParameter(5),
00700              fabs(g_C1->GetParameter(7) - g_C1->GetParameter(7)) / g_C1->GetParameter(8),
00701              fabs(g_C2->GetParameter(1) - g_C2->GetParameter(7)) / g_C2->GetParameter(2),//C2
00702              fabs(g_C2->GetParameter(4) - g_C2->GetParameter(7)) / g_C2->GetParameter(5),
00703              fabs(g_C2->GetParameter(7) - g_C2->GetParameter(7)) / g_C2->GetParameter(8)
00704              );
00705       /*
00706         if (Ch1 > pidcuts[0]  && Ch1 < pidcuts[1])  Ch1El = true;
00707         if (Ch2 > pidcuts[2]  && Ch2 < pidcuts[3])  Ch2El = true;
00708         if (Pb  > pidcuts[4]  && Pb  < pidcuts[5])  PbEl  = true;
00709         if (Ch1 > pidcuts[6]  && Ch1 < pidcuts[7])  Ch1Pi = true;
00710         if (Ch2 > pidcuts[8]  && Ch2 < pidcuts[9])  Ch2Pi = true;
00711         if (Pb  > pidcuts[10] && Pb  < pidcuts[11]) PbPi  = true;
00712       */
00713       printf("$\\pi$ & %8.2f & %8.2f & %8.2f &           & %8.2f & %8.2f & %8.2f \\\\ \n$\\mu$ & %8.2f & %8.2f & %8.2f &           & %8.2f & %8.2f & %8.2f \\\\ \ne     & %8.2f & %8.2f & %8.2f &           & %8.2f & %8.2f & %8.2f \\\\ \n\n",
00714              0.5 * (fabs(pidcuts[6]  - g_C1->GetParameter(1)) + fabs(pidcuts[7] - g_C1->GetParameter(1))) / g_C1->GetParameter(2),//C1 pi
00715              fabs(pidcuts[12] - g_C1->GetParameter(1)) / g_C1->GetParameter(5),
00716              fabs(pidcuts[0]  - g_C1->GetParameter(1)) / g_C1->GetParameter(8),
00717              0.5 * (fabs(pidcuts[8]  - g_C2->GetParameter(1)) + fabs(pidcuts[9] - g_C2->GetParameter(1))) / g_C2->GetParameter(2),//C2
00718              fabs(pidcuts[14] - g_C2->GetParameter(1)) / g_C2->GetParameter(5),
00719              fabs(pidcuts[2]  - g_C2->GetParameter(1)) / g_C2->GetParameter(8),
00720              fabs(pidcuts[7]  - g_C1->GetParameter(4)) / g_C1->GetParameter(2),//C1 my
00721              0.5 * (fabs(pidcuts[12] - g_C1->GetParameter(4)) + fabs(pidcuts[13] - g_C1->GetParameter(4))) / g_C1->GetParameter(5),
00722              fabs(pidcuts[0]  - g_C1->GetParameter(4)) / g_C1->GetParameter(8),
00723              fabs(pidcuts[9]  - g_C2->GetParameter(4)) / g_C2->GetParameter(2),//C2
00724              0.5 * (fabs(pidcuts[14] - g_C2->GetParameter(4)) + fabs(pidcuts[15] - g_C2->GetParameter(4))) / g_C2->GetParameter(5),
00725              fabs(pidcuts[2]  - g_C2->GetParameter(4)) / g_C2->GetParameter(8),
00726              fabs(pidcuts[7]  - g_C1->GetParameter(7)) / g_C1->GetParameter(2),//C1 el
00727              fabs(pidcuts[13] - g_C1->GetParameter(7)) / g_C1->GetParameter(5),
00728              0.5 * (fabs(pidcuts[0]  - g_C1->GetParameter(7)) + fabs(pidcuts[1] - g_C1->GetParameter(7))) / g_C1->GetParameter(8),
00729              fabs(pidcuts[9]  - g_C2->GetParameter(7)) / g_C2->GetParameter(2),//C2
00730              fabs(pidcuts[15] - g_C2->GetParameter(7)) / g_C2->GetParameter(5),
00731              0.5 * (fabs(pidcuts[2]  - g_C2->GetParameter(7)) + fabs(pidcuts[3] - g_C2->GetParameter(7))) / g_C2->GetParameter(8)
00732              );
00733 
00734       Double_t C1_e_in_p = gC1_e->Integral(  pidcuts[6],  pidcuts[7]);
00735       Double_t C2_e_in_p = gC2_e->Integral(  pidcuts[8],  pidcuts[9]);
00736       Double_t Pb_e_in_p = gPb1_e->Integral(pidcuts[10],  Pb_e->GetBinCenter(Pb_e->GetMaximumBin())) + gPb2_e->Integral( Pb_e->GetBinCenter(Pb_e->GetMaximumBin()), pidcuts[11]);
00737       if (Pb_e_in_p > 1.0) Pb_e_in_p = 1.0;
00738       Double_t C1_p_in_e = gC1_p->Integral(  pidcuts[0],  pidcuts[1]);
00739       if (C1_p_in_e == 0.0) C1_p_in_e = 1.0/(allEvents+1);
00740       Double_t C2_p_in_e = gC2_p->Integral(  pidcuts[2],  pidcuts[3]);
00741       Double_t Pb_p_in_e = gPb_p->Integral(  pidcuts[4],  pidcuts[5]);
00742       if (Pb_p_in_e > 1.0) Pb_p_in_e = 1.0;
00743       Double_t e_in_p = Pb_e_in_p * C1_e_in_p * C2_e_in_p;
00744       Double_t p_in_e = Pb_p_in_e * C1_p_in_e * C2_p_in_e;
00745       printf("\n    Pb            Ch1           Ch2           epsilon\ne:  %E  %E  %E  %E  \npi: %E  %E  %E  %E  \n\n",Pb_e_in_p,C1_e_in_p,C2_e_in_p,e_in_p,Pb_p_in_e,C1_p_in_e,C2_p_in_e,p_in_e);
00746 
00747     
00748 
00749       c1->cd(4);
00750       TString text;
00751       TText *t = NULL;
00752       TPaveText *results = new TPaveText(0.10,0.95-0.04*8,0.10,0.95);
00753       results->SetLineColor(0);
00754       results->SetShadowColor(0);
00755       results->SetFillStyle(0);
00756       results->SetTextSize(0.035);
00757       results->AddText("");  
00758       t = results->AddText("           Pb                  Ch1                 Ch2                    #varepsilon");
00759       t->SetTextAlign(11); 
00760       results->AddText("");
00761       text.Form("e:    %.4E    %.4E    %.4E    %.4E  ",epsilon_Cal_e,epsilon_Ch1_e,epsilon_Ch2_e,epsilon_e);
00762       t = results->AddText(text.Data());
00763       t->SetTextAlign(11);
00764       text.Form("#pi:    %.4E    %.4E    %.4E    %.4E  ",epsilon_Cal_pi,epsilon_Ch1_pi,epsilon_Ch2_pi,epsilon_pi);
00765       t = results->AddText(text.Data());
00766       t->SetTextAlign(11);
00767       results->AddText("");
00768       //results->AddLine(.0,.5,1.,.5);
00769       text.Form("e:    %.4E    %.4E    %.4E    %.4E  ",Pb_e_in_p,C1_e_in_p,C2_e_in_p,e_in_p);
00770       t = results->AddText(text.Data());
00771       t->SetTextAlign(11);
00772       text.Form("#pi:    %.4E    %.4E    %.4E    %.4E  ",Pb_p_in_e,C1_p_in_e,C2_p_in_e,p_in_e);
00773       t = results->AddText(text.Data());
00774       t->SetTextAlign(11);
00775       results->AddText("");
00776       results->AddText("");
00777       results->Draw();
00778       c1->Update();
00779       TPaveText *results2 = new TPaveText(0.10,0.65-0.04*10,0.10,0.65);
00780       results2->SetLineColor(0);
00781       results2->SetShadowColor(0);
00782       results2->SetFillStyle(0);
00783       results2->SetTextSize(0.035);
00784       t = results2->AddText("all: ");
00785       t->SetTextAlign(31);
00786       t = results2->AddText("A: ");
00787       t->SetTextAlign(31);
00788       t = results2->AddText("B: ");
00789       t->SetTextAlign(31);
00790       t = results2->AddText("C: ");
00791       t->SetTextAlign(31);
00792       t = results2->AddText("D: ");
00793       t->SetTextAlign(31);
00794       t = results2->AddText("E: ");
00795       t->SetTextAlign(31);
00796       t = results2->AddText("F: ");
00797       t->SetTextAlign(31);
00798       t = results2->AddText("G: ");
00799       t->SetTextAlign(31);
00800       t = results2->AddText("H: ");
00801       t->SetTextAlign(31);
00802       t = results2->AddText("rest: ");
00803       t->SetTextAlign(31);
00804       results2->Draw("same");
00805       c1->Update();
00806       TPaveText *results3 = new TPaveText(0.21,0.65-0.04*10,0.21,0.65);
00807       results3->SetLineColor(0);
00808       results3->SetShadowColor(0);
00809       results3->SetFillStyle(0);
00810       results3->SetTextSize(0.035);
00811       text.Form("%8i ",all);
00812       t = results3->AddText(text.Data());
00813       t->SetTextAlign(31);
00814       text.Form("%8i ",a);
00815       t = results3->AddText(text.Data());
00816       t->SetTextAlign(31);
00817       text.Form("%8i ",b);
00818       t = results3->AddText(text.Data());
00819       t->SetTextAlign(31);
00820       text.Form("%8i ",c);
00821       t = results3->AddText(text.Data());
00822       t->SetTextAlign(31);
00823       text.Form("%8i ",d);
00824       t = results3->AddText(text.Data());
00825       t->SetTextAlign(31);
00826       text.Form("%8i ",e);
00827       t = results3->AddText(text.Data());
00828       t->SetTextAlign(31);
00829       text.Form("%8i ",f);
00830       t = results3->AddText(text.Data());
00831       t->SetTextAlign(31);
00832       text.Form("%8i ",g);
00833       t = results3->AddText(text.Data());
00834       t->SetTextAlign(31);
00835       text.Form("%8i ",h);
00836       t = results3->AddText(text.Data());
00837       t->SetTextAlign(31);
00838       text.Form("%8i ",rest);
00839       t = results3->AddText(text.Data());
00840       t->SetTextAlign(31);
00841       results3->Draw("same");
00842       c1->Update();
00843       TPaveText *results4 = new TPaveText(0.34,0.65-0.04*10,0.34,0.65);
00844       results4->SetLineColor(0);
00845       results4->SetShadowColor(0);
00846       results4->SetFillStyle(0);
00847       results4->SetTextSize(0.035);
00848       text.Form("%6.2f%%",100.*all/all);
00849       t = results4->AddText(text.Data());
00850       t->SetTextAlign(31);
00851       text.Form("%6.2f%%",100.*a/all);
00852       t = results4->AddText(text.Data());
00853       t->SetTextAlign(31);
00854       text.Form("%6.2f%%",100.*b/all);
00855       t = results4->AddText(text.Data());
00856       t->SetTextAlign(31);
00857       text.Form("%6.2f%%",100.*c/all);
00858       t = results4->AddText(text.Data());
00859       t->SetTextAlign(31);
00860       text.Form("%6.2f%%",100.*d/all);
00861       t = results4->AddText(text.Data());
00862       t->SetTextAlign(31);
00863       text.Form("%6.2f%%",100.*e/all);
00864       t = results4->AddText(text.Data());
00865       t->SetTextAlign(31);
00866       text.Form("%6.2f%%",100.*f/all);
00867       t = results4->AddText(text.Data());
00868       t->SetTextAlign(31);
00869       text.Form("%6.2f%%",100.*g/all);
00870       t = results4->AddText(text.Data());
00871       t->SetTextAlign(31);
00872       text.Form("%6.2f%%",100.*h/all);
00873       t = results4->AddText(text.Data());
00874       t->SetTextAlign(31);
00875       text.Form("%6.2f%%",100.*rest/all);
00876       t = results4->AddText(text.Data());
00877       t->SetTextAlign(31);
00878       results4->Draw("same");
00879       /*
00880         text.Form("   all: %08i (%06.2f%%)",all,100.*all/all);
00881         t = results->AddText(text.Data());
00882         t->SetTextAlign(11);
00883         text.Form("    A: %08i (%06.2f%%)",a,100.*a/all);
00884         t = results->AddText(text.Data());
00885         t->SetTextAlign(11);
00886         text.Form("    B: %08i (%06.2f%%)",b,100.*b/all);
00887         t = results->AddText(text.Data());
00888         t->SetTextAlign(11);
00889         text.Form("    C: %08i (%06.2f%%)",c,100.*c/all);
00890         t = results->AddText(text.Data());
00891         t->SetTextAlign(11);
00892         text.Form("    D: %08i (%06.2f%%)",d,100.*d/all);
00893         t = results->AddText(text.Data());
00894         t->SetTextAlign(11);
00895         text.Form("    E: %08i (%06.2f%%)",e,100.*e/all);
00896         t = results->AddText(text.Data());
00897         t->SetTextAlign(11);
00898         text.Form("    F: %08i (%06.2f%%)",f,100.*f/all);
00899         t = results->AddText(text.Data());
00900         t->SetTextAlign(11);
00901         text.Form("    G: %08i (%06.2f%%)",g,100.*g/all);
00902         t = results->AddText(text.Data());
00903         t->SetTextAlign(11);
00904         text.Form("    H: %08i (%06.2f%%)",h,100.*h/all);
00905         t = results->AddText(text.Data());
00906         t->SetTextAlign(11);
00907         text.Form("rest: %08i (%06.2f%%)",rest,100.*rest/all);
00908         t = results->AddText(text.Data());
00909         t->SetTextAlign(11);
00910       */
00911  
00912    
00913       c1->Update();
00914       outpic.ReplaceAll(".",",");
00915       c1->SaveAs("data2012/pics/Purity/Spectra_"+outpic+".pdf");
00916       c1->SaveAs("data2012/pics/Purity/Spectra_"+outpic+".png");
00917 
00918       /*
00919         printf("\nPb_p_min: %i, Pb_e_center: %.0f, Pb_e_center: %.0f, Pb_p_max: %i\n   Int1: %f   Int2: %f  Int1+2: %f\n",
00920         pidcuts[10],  
00921         Pb_e->GetBinCenter(Pb_e->GetMaximumBin()), 
00922         Pb_e->GetBinCenter(Pb_e->GetMaximumBin()), 
00923         pidcuts[11],
00924         gPb1_e->Integral(pidcuts[10], Pb_e->GetBinCenter(Pb_e->GetMaximumBin())),
00925         gPb2_e->Integral( Pb_e->GetBinCenter(Pb_e->GetMaximumBin()), pidcuts[11]),
00926         (gPb1_e->Integral(pidcuts[10],  Pb_e->GetBinCenter(Pb_e->GetMaximumBin())) + gPb2_e->Integral( Pb_e->GetBinCenter(Pb_e->GetMaximumBin()), pidcuts[11]))
00927         );
00928       */
00929       c2->Write("", TObject::kOverwrite);
00930       c1->Write("", TObject::kOverwrite);
00931       gC1_e->Write("", TObject::kOverwrite);
00932       gC2_e->Write("", TObject::kOverwrite);
00933       gPb1_e->Write("", TObject::kOverwrite);
00934       gPb2_e->Write("", TObject::kOverwrite);
00935       gC1_p->Write("", TObject::kOverwrite);
00936       gC2_p->Write("", TObject::kOverwrite);
00937       gPb_p->Write("", TObject::kOverwrite);
00938       allC1->Write("", TObject::kOverwrite);
00939       allC2->Write("", TObject::kOverwrite);
00940       allPb->Write("", TObject::kOverwrite);
00941       restC1->Write("", TObject::kOverwrite);
00942       restC2->Write("", TObject::kOverwrite);
00943       restPb->Write("", TObject::kOverwrite);
00944       C1_e->Write("", TObject::kOverwrite);
00945       C1_p->Write("", TObject::kOverwrite);
00946       C1_m->Write("", TObject::kOverwrite);
00947       C2_e->Write("", TObject::kOverwrite);
00948       C2_p->Write("", TObject::kOverwrite);
00949       C2_m->Write("", TObject::kOverwrite);
00950       Pb_e->Write("", TObject::kOverwrite);
00951       Pb_p->Write("", TObject::kOverwrite);
00952       Pb_m->Write("", TObject::kOverwrite);
00953       allPIDs->Write("", TObject::kOverwrite);
00954       RestPIDs->Write("", TObject::kOverwrite);
00955       APIDs->Write("", TObject::kOverwrite);
00956       BPIDs->Write("", TObject::kOverwrite);
00957       CPIDs->Write("", TObject::kOverwrite);
00958       DPIDs->Write("", TObject::kOverwrite);
00959       EPIDs->Write("", TObject::kOverwrite);
00960       FPIDs->Write("", TObject::kOverwrite);
00961       GPIDs->Write("", TObject::kOverwrite);
00962       HPIDs->Write("", TObject::kOverwrite);      
00963       output->Close();
00964     }
00965   }
00966 }
beamtime/cern-oct12/go4/AnalysisMacros/spadic_PIDpurity_Analysis.cxx (r4864/r4848)
