2 class AliSpectraBothHistoManager;
3 class AliSpectraBothEventCuts;
4 class AliSpectraBothTrackCuts;
5 TString Charge[]={"Pos","Neg"};
6 TString Sign[]={"Plus","Minus"};
7 TString Particle[]={"Pion","Kaon","Proton"};
8 TString symboles[]={"#pi^{+}","K^{+}","p","pi^{-}","K^{-}","#bar{p}"};
9 AliSpectraBothHistoManager* managerdata=0x0;
10 AliSpectraBothEventCuts* ecutsdata=0x0;
11 AliSpectraBothTrackCuts* tcutsdata=0x0;
13 AliSpectraBothHistoManager* managermc=0x0;
14 AliSpectraBothEventCuts* ecutsmc=0x0;
15 AliSpectraBothTrackCuts* tcutsmc=0x0;
17 Float_t TOFMatchingScalling[2]={-1,-1};
18 Int_t Color[3]={1,2,4};
19 Int_t Marker[6]={20,21,22,24,25,26};
20 Double_t Range[3]={0.3,0.3,0.5}; // LowPt range for pi k p
23 Double_t FitRange[2]={-2.0,2.0};
24 Double_t CutRange[2]={-3.0,3.0};
25 Double_t minptformaterial[6]={0.0,0.2,0.0,0.0,0.2,0.0};
26 Double_t maxptformaterial[6]={0.0,0.6,1.3,0.0,0.6,0.0};
27 Double_t minptforWD[6]={0.2,100.0,0.3,0.2,100.0,0.3};
28 Double_t maxptforWD[6]={1.5,-100.0,2.0,1.5,-100.0,2.0};
29 Double_t minRanges[3]={0.3,0.3,0.45};
30 Double_t maxRanges[3]={1.5,1.2,2.2};
31 Double_t TOFPIDsignalmatching[]={-1.0,-1.0,-1.0};
32 Double_t fMaxContaminationPIDMC=0.2;
33 TString filenames[]={"eff.root","pid.root","sec.root"};
41 kdodca=0x1, //dca fits are made
42 kgeantflukaKaon=0x2,// geant fluka correction is used for kaons
43 kgeantflukaProton=0x4, // geant fluka correction is used for protons and antiprotons
44 knormalizationtoeventspassingPhySel=0x8,// spectra are divided by number of events passing physic selection
45 kveretxcorrectionandbadchunkscorr=0x10, // correction for difference in z vertex distribution in data and MC and correction for bad chunks is applied
46 kmcisusedasdata=0x20, // the result of the looping over MC is used as data input
47 kdonotusedcacuts=0x40, // allows to use the constant dca cut for all pt bins not the pt dependet defined in stardrad track cuts 2011
48 kuseprimaryPIDcont=0x80, //pid contamination is calculated using only primiary particle in this case K should use dca fits
49 knormalizationwithbin0integralsdata=0x100, // the normalization factor is calcualte using integral over z vertex distributions (in this case reconstructed vertex disitrbution uses z vertex for data)
50 knormalizationwithbin0integralsMC=0x200, //in this case reconstructed vertex disitrbution uses z vertex for data, those to options will be use only if knormalizationtoeventspassingPhySel is not set
51 kuserangeonfigfile=0x400, // use of config file for dca fit settings
52 kskipconcutonspectra=0x800, //do not use conPID<02 cut useful for syst. studies
53 kuseTOFmatchingcorrection=0x1000, // if set tof matching correction is applied.
54 kuseTOFcorrforPIDsignalmatching=0x2000, // rescale the for spectra by the factor given in config files
55 kuseeffcorrectionfromfile=0x4000, //use the efficiency from the file specfied in config file
56 kusePIDcontaminatiofromfile=0x8000, //use the PID contamination from the file specfied in config file
57 kuseseccontaminatiofromfile=0x10000 //use the secondary contamination from the file specfied in config file
61 Bool_t OpenFile(TString dirname, TString outputname, Bool_t mcflag,Bool_t mcasdata=false);
62 void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TString outnamemc="",TString configfile="" )
64 gStyle->SetOptStat(0);
65 TH1::AddDirectory(kFALSE);
66 gSystem->Load("libCore");
67 gSystem->Load("libPhysics");
68 gSystem->Load("libTree");
69 gSystem->Load("libMatrix");
70 gSystem->Load("libSTEERBase");
71 gSystem->Load("libESD");
72 gSystem->Load("libAOD");
73 gSystem->Load("libANALYSIS");
74 gSystem->Load("libOADB");
75 gSystem->Load("libANALYSISalice");
76 gSystem->Load("libTender");
77 gSystem->Load("libCORRFW");
78 gSystem->Load("libPWGTools");
79 gSystem->Load("libPWGLFspectra");
81 gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/QAPlotsBoth.C");
83 mass[0] = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
84 mass[1] = TDatabasePDG::Instance()->GetParticle("K+")->Mass();
85 mass[2] = TDatabasePDG::Instance()->GetParticle("proton")->Mass();
87 TFormula* dcacutxy=0x0;
88 if(!(options&kdonotusedcacuts))
91 AliESDtrackCuts* esdtrackcuts= AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE);
92 TString formulastring(esdtrackcuts->GetMaxDCAToVertexXYPtDep());
93 formulastring.ReplaceAll("pt","x");
94 dcacutxy=new TFormula("dcacutxy",formulastring.Data());
96 if(options&kuserangeonfigfile)
97 if(!ReadConfigFile(configfile))
99 TList* lout=new TList();
102 TString indirname=Form("/output/train%s",outdate.Data());
103 //TString indirname("/output/train24072012");
104 if(outnamemc.Length()==0)
105 outnamemc=outnamedata;
106 cout<<indirname.Data()<<" "<<outnamemc.Data()<<endl;
110 OpenFile(indirname,outnamemc,true);
111 OpenFile(indirname,outnamedata,false,((Bool_t)(options&kmcisusedasdata)));
112 if(!managermc||!managerdata)
114 cout<<managermc<<" "<<managerdata<<endl;
117 TH1F* rawspectradata[6];
118 TH1F* rawspectramc[6];
131 TH1F* contMatfit[12];
132 TH1F* primaryfit[12];
137 TH1F* spectraLeonardo[6];
139 TH1F* corrLeonardo[6];
140 TH1F* doubleconuntsdata[3];
141 TH1F* doubleconuntsMC[3];
142 //GetSpectra(managerdata,rawspectradata,true);
143 //GetSpectra(managermc,rawspectramc,true,true);
145 GetPtHistFromPtDCAhisto("hHistPtRecSigma","SpectraMC",managermc,rawspectramc,dcacutxy);
146 GetPtHistFromPtDCAhisto("hHistPtRecSigma","SpectraDATA",managerdata,rawspectradata,dcacutxy);
147 GetPtHistFromPtDCAhisto("hHistPtRecTruePrimary","eff",managermc,eff,dcacutxy);
148 GetPtHistFromPtDCAhisto("hHistPtRecTrue","conPID",managermc,contPID,dcacutxy);
149 GetPtHistFromPtDCAhisto("hHistPtRecSigmaSecondaryWeakDecay","conWD",managermc,contWD,dcacutxy);
150 GetPtHistFromPtDCAhisto("hHistPtRecSigmaSecondaryMaterial","conMat",managermc,contMat,dcacutxy);
151 GetPtHistFromPtDCAhisto("hHistPtRecSigmaPrimary","conPIDprimary",managermc,contPIDpri,dcacutxy);
154 Double_t neventsmcall = 1 ; //if loop over MC is done after or befor events cuts this will be changed
155 Double_t neventsdata = 1;
156 Double_t neventsmc = 1;
158 //Normaliztion of MCtruth depends if the loop was done after of before ESD event cuts.
159 //In currect code this cannot be check on the level of macro.
160 //If the loop was done before MC should be done to all processed events (NumberOfProcessedEvents())
161 //If loop was done after MC should be normalized to all accepted events (NumberOfEvents())
162 // The option one will be alaways use.
164 neventsmcall= ecutsmc->NumberOfProcessedEvents();
165 if(options&knormalizationtoeventspassingPhySel)
167 //neventsmcall= ecutsmc->NumberOfProcessedEvents();
168 neventsdata=ecutsdata->NumberOfPhysSelEvents();
169 neventsmc=ecutsmc->NumberOfPhysSelEvents();
171 else if ((options&knormalizationwithbin0integralsdata)||(options&knormalizationwithbin0integralsMC))
173 neventsdata=Normaliztionwithbin0integrals(options);
174 neventsmc=ecutsmc->NumberOfPhysSelEvents();
178 neventsdata=ecutsdata->NumberOfEvents(); //number of accepted events
179 neventsmc=ecutsmc->NumberOfEvents();
180 neventsmcall= ecutsmc->NumberOfEvents();
184 cout<<neventsdata<<" Events"<<endl;
185 cout<< neventsmc<<" Events "<<endl;
188 cout<<"No events DATA"<<endl;
191 if(neventsmc<1&&(((options&kusePIDcontaminatiofromfile)!=kusePIDcontaminatiofromfile)||((options&kuseseccontaminatiofromfile)!=kuseseccontaminatiofromfile)||((options&kuseseccontaminatiofromfile)!=kuseseccontaminatiofromfile)))
193 cout<<"No events MC"<<endl;
198 neventsmc=1; //dumpy normalization
199 neventsmcall=1; //dumpy normalization
203 cout<<"A "<<neventsmc<<endl;
204 TH1F* allgen=((TH1F*)managermc->GetPtHistogram1D("hHistPtGen",1,1))->Clone();
205 allgen->SetName("AllGen");
206 TH1F* allrecMC=GetOneHistFromPtDCAhisto("hHistPtRec","rawallMC",managermc,dcacutxy);
207 TH1F* alleff=GetOneHistFromPtDCAhisto("hHistPtRecPrimary","effall",managermc,dcacutxy);
208 TH1F* allrecdata=GetOneHistFromPtDCAhisto("hHistPtRec","rawalldata",managerdata,dcacutxy);
213 TH1F* spectraall=(TH1F*)allrecdata->Clone("recNch");
214 spectraall->SetTitle("recNch");
215 TH1F* contall=(TH1F*)allrecMC->Clone("contall");
216 contall->SetTitle("contall");
217 //contall->Add(alleff,-1);
218 SubHistWithFullCorr(contall,alleff);
219 alleff->Divide(alleff,allgen,1,1,"B");
220 contall->Divide(contall,allrecMC,1,1,"B");
222 GetCorrectedSpectra(spectraall,allrecdata,alleff,contall);
223 Divideby2pipt(spectraall);
225 allrecdata->Scale(1./neventsdata,"width");
226 allgen->Scale(1./neventsmcall,"width");
227 allrecMC->Scale(1./neventsmc,"width");
228 spectraall->Scale(1./neventsdata,"width");
234 lout->Add(allrecdata);
235 lout->Add(spectraall);
238 for (int i=0;i<6;i++)
242 TString tmpname(rawspectramc[i]->GetTitle());
243 tmpname.ReplaceAll("SpectraMC","%s");
244 contallMC[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contallMC"));
245 contallMC[i]->SetTitle(Form(tmpname.Data(),"contallMC"));
246 contfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contfit"));
247 contfit[i]->SetTitle(Form(tmpname.Data(),"contfit"));
248 contWDfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contWDfit"));
249 contWDfit[i]->SetTitle(Form(tmpname.Data(),"contWDfit"));
250 contMatfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contMatfit"));
251 contMatfit[i]->SetTitle(Form(tmpname.Data(),"contMatfit"));
252 primaryfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"primaryfit"));
253 primaryfit[i]->SetTitle(Form(tmpname.Data(),"primaryfit"));
254 contfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contfitonMC"));
255 contfit[i+6]->SetTitle(Form(tmpname.Data(),"contfitonMC"));
256 contWDfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contWDfitonMC"));
257 contWDfit[i+6]->SetTitle(Form(tmpname.Data(),"contWDfitonMC"));
258 contMatfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contMatfitonMC"));
259 contMatfit[i+6]->SetTitle(Form(tmpname.Data(),"contMatfitonMC"));
260 primaryfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"primaryfitMC"));
261 primaryfit[i+6]->SetTitle(Form(tmpname.Data(),"primaryfitMC"));
267 contWDfit[i]->Reset();
268 contMatfit[i]->Reset();
269 primaryfit[i]->Reset();
272 contfit[i+6]->Reset();
273 contWDfit[i+6]->Reset();
274 contMatfit[i+6]->Reset();
275 primaryfit[i+6]->Reset();
277 SetBintoOne(primaryfit[i]);
278 SetBintoOne(primaryfit[i+6]);
279 spectra[i]=(TH1F*)rawspectradata[i]->Clone(Form(tmpname.Data(),"SpectraFinal"));
280 spectra[i]->SetTitle(Form(tmpname.Data(),"SpectraFinal"));
281 spectraLeonardo[i]=(TH1F*)rawspectradata[i]->Clone(Form(tmpname.Data(),"SpectraFinalLeonardo"));
282 spectraLeonardo[i]->SetTitle(Form(tmpname.Data(),"SpectraFinalLeonardo"));
284 corrLeonardo[i]=(TH1F*)MCTruth[i]->Clone(Form(tmpname.Data(),"CorrFactLeonardo"));
285 corrLeonardo[i]->SetTitle(Form(tmpname.Data(),"CorrFactLeonardo"));
286 corrLeonardo[i]->Divide(corrLeonardo[i],rawspectramc[i],1,1,"B");
289 if((options&kusePIDcontaminatiofromfile)==kusePIDcontaminatiofromfile)
291 CopyCorrectionFromFile(filenames[1],"conPID",contPID);
292 CopyCorrectionFromFile(filenames[1],"conPIDprimary",contPIDpri);
297 //contallMC[i]->Add(eff[i],-1.0);
298 SubHistWithFullCorr(contallMC[i],eff[i]);
299 //RecomputeErrors(contallMC[i]);
300 contallMC[i]->Sumw2();
301 contallMC[i]->Divide(contallMC[i],rawspectramc[i],1,1,"B");
302 // contamintaion from PID but only primaries
303 //contPIDpri[i]->Add(eff[i],-1.0);
304 SubHistWithFullCorr(contPIDpri[i],eff[i]);
305 //RecomputeErrors(contPIDpri[i]);
306 contPIDpri[i]->Divide(contPIDpri[i],rawspectramc[i],1,1,"B");
309 if((options&kuseeffcorrectionfromfile)==kuseeffcorrectionfromfile)
310 CopyCorrectionFromFile(filenames[0],"eff",eff);
312 eff[i]->Divide(eff[i],MCTruth[i],1,1,"B");
316 rawspectramc[i]->Sumw2();
317 //contPID[i]->Add(contPID[i],rawspectramc[i],-1,1);
318 SubHistWithFullCorr(contPID[i],rawspectramc[i]);
319 contPID[i]->Scale(-1.0);
321 //RecomputeErrors(contPID[i]);
322 contPID[i]->ResetStats();
324 contPID[i]->Divide(contPID[i],rawspectramc[i],1,1,"B");
326 if((options&kuseprimaryPIDcont)==kuseprimaryPIDcont)
327 confinal[i]=(TH1F*)contPIDpri[i]->Clone(Form(tmpname.Data(),"confinal"));
329 confinal[i]=(TH1F*)contPID[i]->Clone(Form(tmpname.Data(),"confinal"));
330 confinal[i]->SetTitle(Form(tmpname.Data(),"confinal"));
332 contSecMC[i]=(TH1F*)contWD[i]->Clone(Form(tmpname.Data(),"conSecMC"));
333 contSecMC[i]->Add(contMat[i]);
334 contWD[i]->Divide(contWD[i],rawspectramc[i],1,1,"B");
335 contMat[i]->Divide(contMat[i],rawspectramc[i],1,1,"B");
336 contSecMC[i]->Divide(contSecMC[i],rawspectramc[i],1,1,"B");
339 doubleconuntsdata[i-3]=(TH1F*)rawspectradata[i]->Clone(Form("DoublecountsDATA%s",Particle[i-3].Data()));
340 doubleconuntsdata[i-3]->Reset();
341 doubleconuntsMC[i-3]=(TH1F*)rawspectramc[i]->Clone(Form("DoublecountsMC%s",Particle[i-3].Data()));
342 doubleconuntsMC[i-3]->Reset();
344 CalculateDoubleCounts(doubleconuntsdata[i-3],rawspectradata,i-3,kTRUE);
345 CalculateDoubleCounts(doubleconuntsMC[i-3],rawspectramc,i-3,kFALSE);
353 rawspectradata[i]->Scale(1./neventsdata,"width");
354 rawspectramc[i]->Scale(1./neventsmc,"width");
355 MCTruth[i]->Scale(1./neventsmcall,"width");
356 spectraLeonardo[i]->Scale(1./neventsdata,"width");
360 lout->Add(rawspectradata[i]);
361 lout->Add(rawspectramc[i]);
362 lout->Add(MCTruth[i]);
364 lout->Add(contallMC[i]);
365 lout->Add(contPID[i]);
366 lout->Add(contWD[i]);
367 lout->Add(contMat[i]);
368 lout->Add(contfit[i]);
369 lout->Add(contWDfit[i]);
370 lout->Add(contMatfit[i]);
371 lout->Add(primaryfit[i]);
372 lout->Add(contfit[i+6]);
373 lout->Add(contWDfit[i+6]);
374 lout->Add(contMatfit[i+6]);
375 lout->Add(primaryfit[i+6]);
376 lout->Add(spectra[i]);
377 lout->Add(spectraLeonardo[i]);
378 lout->Add(confinal[i]);
379 lout->Add(contPIDpri[i]);
380 lout->Add(contSecMC[i]);
383 lout->Add(doubleconuntsdata[i-3]);
384 lout->Add(doubleconuntsMC[i-3]);
387 outdate.ReplaceAll("/","_");
388 configfile.ReplaceAll(".","_");
390 if(configfile.Length()>0&&(options&kuserangeonfigfile))
391 fout=new TFile(Form("./results/ResMY_%s_%s_%#X_%s.root",outnamemc.Data(),outdate.Data(),options,configfile.Data()),"RECREATE");
393 fout=new TFile(Form("./results/ResMY_%s_%s_%#X.root",outnamemc.Data(),outdate.Data(),options),"RECREATE");
394 if (((options&kdodca)==kdodca)&&((options&kuseseccontaminatiofromfile)!=kuseseccontaminatiofromfile))
395 DCACorrectionMarek(managerdata,managermc,dcacutxy,fout,contfit,contWDfit,contMatfit,primaryfit);
396 else if ((options&kuseseccontaminatiofromfile)==kuseseccontaminatiofromfile)
398 CopyCorrectionFromFile(filenames[2],"contfit",contfit);
399 CopyCorrectionFromFile(filenames[2],"contWDfit",contWDfit);
400 CopyCorrectionFromFile(filenames[2],"contMatfit",contMatfit);
401 CopyCorrectionFromFile(filenames[2],"primaryfit",primaryfit);
405 cout<<"Secondary from MC"<<endl;
408 for (int i=0;i<6;i++)
410 if(((options&kdodca)==kdodca)||((options&kuseseccontaminatiofromfile)==kuseseccontaminatiofromfile))
412 if((options&kuseprimaryPIDcont)||(i!=1&&i!=4)) //if we do not use cont PId only for primary and this is a kaon that do not use fit
413 confinal[i]->Add(contfit[i]);
414 GetCorrectedSpectra(spectra[i],rawspectradata[i],eff[i],confinal[i]);
418 GetCorrectedSpectra(spectra[i],rawspectradata[i],eff[i],contallMC[i]);
420 GetCorrectedSpectraLeonardo(spectraLeonardo[i],corrLeonardo[i],primaryfit[i],primaryfit[i+6]);
421 if(options&kskipconcutonspectra)
423 if(options&kuseprimaryPIDcont)
425 CleanHisto(spectra[i],-1,100,contPIDpri[i]);
426 CleanHisto(spectraLeonardo[i],-1,100,contPIDpri[i]);
430 CleanHisto(spectra[i],-1,100,contPID[i]);
431 CleanHisto(spectraLeonardo[i],-1,100,contPID[i]);
433 // Apply correction for wrongly simulated TOF signal in MC
434 if(options&kuseTOFcorrforPIDsignalmatching)
435 TOFPIDsignalmatchingApply(spectra[i],TOFPIDsignalmatching[i%3]);
439 GFCorrection(spectra,tcutsdata->GetPtTOFMatching(),options);
440 GFCorrection(spectraLeonardo,tcutsdata->GetPtTOFMatching(),options);
442 Double_t ycut=tcutsdata->GetY();
443 Double_t etacut=tcutsdata->GetEtaMax()-tcutsdata->GetEtaMin();
446 cout<<"using eta window 1.6"<<endl;
451 TH1F* allch=GetSumAllCh(spectra,mass,etacut);
453 if(options&kuseTOFmatchingcorrection)
455 MatchingTOFEff(spectra,lout);
456 MatchingTOFEff(spectraLeonardo);
457 TOFMatchingForNch(spectraall);
461 TList* listqa=new TList();
462 TList* canvaslist=new TList();
463 Float_t vertexcorrection=1.0;
464 Float_t corrbadchunksvtx=1.0;
465 if ((options&knormalizationwithbin0integralsdata)||(options&knormalizationwithbin0integralsMC))
466 corrbadchunksvtx=QAPlotsBoth(managerdata,managermc,ecutsdata,ecutsmc,tcutsdata,tcutsmc,listqa,canvaslist,0);
468 corrbadchunksvtx=QAPlotsBoth(managerdata,managermc,ecutsdata,ecutsmc,tcutsdata,tcutsmc,listqa,canvaslist,1);
469 if (options&kveretxcorrectionandbadchunkscorr)
470 vertexcorrection=corrbadchunksvtx;
471 cout<<" VTX corr="<<vertexcorrection<<endl;
472 if(TMath::Abs(ycut)>0.0)
473 vertexcorrection=vertexcorrection/(2.0*ycut);
474 for (int i=0;i<6;i++)
476 spectra[i]->Scale(vertexcorrection);
477 spectraLeonardo[i]->Scale(vertexcorrection);
478 if(TMath::Abs(ycut)>0.0)
480 rawspectradata[i]->Scale(1.0/(2.0*ycut));
481 rawspectramc[i]->Scale(1.0/(2.0*ycut));
482 MCTruth[i]->Scale(1.0/(2.0*ycut));
485 allch->Scale(vertexcorrection);
486 spectraall->Scale(vertexcorrection/etacut);
488 //spectraall->Scale(1.0/1.6);
489 lout->Write("output",TObject::kSingleKey);
490 listqa->Write("outputQA",TObject::kSingleKey);
491 canvaslist->Write("outputcanvas",TObject::kSingleKey);
494 //Normaliztionwithbin0integrals();
498 Bool_t OpenFile(TString dirname,TString outputname, Bool_t mcflag, Bool_t mcasdata)
502 TString nameFile = Form("./%s/AnalysisResults%s.root",dirname.Data(),(mcflag?"MC":"DATA"));
503 TFile *file = TFile::Open(nameFile.Data());
506 cout<<"no file"<<endl;
509 TString sname=Form("OutputBothSpectraTask_%s_%s",(mcflag?"MC":"Data"),outputname.Data());
512 cout<<"using MC as data "<<endl;
513 sname=Form("OutputBothSpectraTask_%s_%s","MC",outputname.Data());
516 TDirectoryFile *dir=(TDirectoryFile*)file->Get(sname.Data());
519 // cout<<"no dir "<<sname.Data()<<endl;
522 cout<<"using MC as data "<<endl;
523 sname=Form("OutputAODSpectraTask_%s_%s","MC",outputname.Data());
526 sname=Form("OutputAODSpectraTask_%s_%s",(mcflag?"MC":"Data"),outputname.Data());
527 // cout<<"trying "<<sname.Data()<<endl;
528 dir=(TDirectoryFile*)file->Get(sname.Data());
531 cout<<"no dir "<<sname.Data()<<endl;
535 cout << " -- Info about " <<(mcflag?"MC":"DATA") <<" -- "<< endl;
538 managermc= (AliSpectraBothHistoManager*) dir->Get("SpectraHistos");
539 ecutsmc = (AliSpectraBothEventCuts*) dir->Get("Event Cuts");
540 tcutsmc = (AliSpectraBothTrackCuts*) dir->Get("Track Cuts");
541 ecutsmc->PrintCuts();
542 tcutsmc->PrintCuts();
543 if(!managermc||!ecutsmc||!tcutsmc)
545 if(managermc->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()!=ecutsmc->GetHistoCuts()->GetBinContent(3))
546 cout<<"Please check MC file possible problem with merging"<<" "<<managermc->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()<<" "<<ecutsmc->GetHistoCuts()->GetBinContent(3)<<endl;
550 managerdata= (AliSpectraBothHistoManager*) dir->Get("SpectraHistos");
551 ecutsdata = (AliSpectraBothEventCuts*) dir->Get("Event Cuts");
552 tcutsdata = (AliSpectraBothTrackCuts*) dir->Get("Track Cuts");
553 ecutsdata->PrintCuts();
554 tcutsdata->PrintCuts();
555 if(!managerdata||!ecutsdata||!tcutsdata)
557 if(managerdata->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()!=ecutsdata->GetHistoCuts()->GetBinContent(3))
558 cout<<"Please check DATA file possible problem with merging"<<" "<<managerdata->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()<<" "<<ecutsdata->GetHistoCuts()->GetBinContent(3)<<endl;
564 void GetMCTruth(TH1F** MCTruth)
566 for(Int_t icharge=0;icharge<2;icharge++)
568 for(Int_t ipart=0;ipart<3;ipart++)
570 Int_t index=ipart+3*icharge;
571 TString hname=Form("hHistPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
572 MCTruth[index]=(TH1F*)((TH1F*)managermc->GetPtHistogram1D(hname.Data(),1,1))->Clone();
573 MCTruth[index]->SetName(Form("MCTruth_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
574 MCTruth[index]->SetTitle(Form("MCTruth_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
575 MCTruth[index]->Sumw2();
580 TH1F* GetOneHistFromPtDCAhisto(TString name,TString hnameout,AliSpectraBothHistoManager* hman,TFormula* dcacutxy)
582 histo =(TH1F*)((TH1F*) hman->GetPtHistogram1D(name.Data(),-1,-1))->Clone();
583 histo->SetName(hnameout.Data());
584 histo->SetTitle(hnameout.Data());
588 for(int ibin=1;ibin<histo->GetNbinsX();ibin++)
590 Double_t lowedge=histo->GetBinLowEdge(ibin);
591 Float_t cut=dcacutxy->Eval(lowedge);
592 TH1F* dcahist=(TH1F*)hman->GetDCAHistogram1D(name.Data(),lowedge,lowedge);
593 // Float_t inyield=dcahist->Integral(dcahist->GetXaxis()->FindBin(-1.0*cut),dcahist->GetXaxis()->FindBin(cut));
594 Float_t testyield=0.0;
595 Float_t testerror=0.0;
596 for (int itest=dcahist->GetXaxis()->FindBin(-1.0*cut);itest<=dcahist->GetXaxis()->FindBin(cut);itest++)
598 testyield+=dcahist->GetBinContent(itest);
599 testerror+=dcahist->GetBinError(itest)*dcahist->GetBinError(itest);
601 // cout<<"corr data "<<histo->GetBinContent(ibin)<<" "<<inyield<<" "<<dcahist->Integral()<<" "<<hnameout.Data()<<endl;
602 // cout<<"test dca "<<lowedge<<" "<<dcacutxy->Eval(lowedge)<<" "<<dcacutxy->Eval(histo->GetXaxis()->GetBinUpEdge(ibin))<<" "<<dcahist->GetBinLowEdge(dcahist->GetXaxis()->FindBin(-1.0*cut))<<" "<<dcahist->GetXaxis()->GetBinUpEdge(dcahist->GetXaxis()->FindBin(-1.0*cut))<<endl;
604 // cout<<testyield<<" "<<TMath::Sqrt(testerror)<<" error2 "<<inyield<<" "<<TMath::Sqrt(inyield)<<endl;
606 //histo->SetBinContent(ibin,inyield);
607 //histo->SetBinError(ibin,TMath::Sqrt(inyield));
608 histo->SetBinContent(ibin,testyield);
609 histo->SetBinError(ibin,TMath::Sqrt(testerror));
620 void GetPtHistFromPtDCAhisto(TString hnamein, TString hnameout, AliSpectraBothHistoManager* hman,TH1F** histo,TFormula* dcacutxy)
622 //Float_t min[3]={0.3,0.3,0.45};
623 //Float_t max[3]={1.5,1.2,2.2};
624 for(Int_t icharge=0;icharge<2;icharge++)
626 for(Int_t ipart=0;ipart<3;ipart++)
628 Int_t index=ipart+3*icharge;
629 Printf("Getting %s",hnamein.Data());
630 TString nameinfinal=Form("%s%s%s",hnamein.Data(),Particle[ipart].Data(),Sign[icharge].Data());
631 TString nameoutfinal=Form("%s%s%s",hnameout.Data(),Particle[ipart].Data(),Sign[icharge].Data());
634 histo[index]=GetOneHistFromPtDCAhisto(nameinfinal,nameoutfinal,hman,dcacutxy);
635 CleanHisto(histo[index],minRanges[ipart],maxRanges[ipart]);
639 void CleanHisto(TH1F* h, Float_t minV, Float_t maxV,TH1* contpid=0x0)
641 for (int i=0;i<=h->GetNbinsX();i++)
643 if(h->GetXaxis()->GetBinCenter(i)<minV||h->GetXaxis()->GetBinCenter(i)>maxV)
645 h->SetBinContent(i,0);
650 if(contpid->GetBinContent(i)>fMaxContaminationPIDMC)
652 h->SetBinContent(i,0);
660 void DCACorrectionMarek(AliSpectraBothHistoManager* hman_data, AliSpectraBothHistoManager* hman_mc,TFormula* fun,TFile *fout,TH1F** hcon,TH1F** hconWD,TH1F** hconMat,TH1F** hprimary)
662 printf("\n\n-> DCA Correction");
666 TString sample[2]={"data","mc"};
667 ofstream debug("debugDCA.txt");
668 TList* listofdcafits=new TList();
669 for(Int_t icharge=0;icharge<2;icharge++)
671 for(Int_t ipart=0;ipart<3;ipart++)
673 Int_t index=ipart+3*icharge;
674 for(Int_t isample=0;isample<2;isample++)
677 for(Int_t ibin_data=7;ibin_data<40;ibin_data++)
680 Double_t lowedge=hcon[index]->GetBinLowEdge(ibin_data);
681 Double_t binwidth=hcon[index]->GetBinWidth(ibin_data);
682 debug<<"NEW "<<Particle[ipart].Data()<<" "<<Sign[icharge].Data()<<" "<<lowedge<<endl;
685 CutRange[1]=fun->Eval(lowedge);
686 CutRange[0]=-1.0*CutRange[1];
688 debug<<"cut "<<CutRange[1]<<" "<<CutRange[0]<<endl;
689 Short_t fitsettings=DCAfitsettings(lowedge+0.5*binwidth,index);
690 debug<<"settings "<< fitsettings<<endl;
694 TCanvas *cDCA=new TCanvas(Form("cDCA%d%s%s%sbin%d",index,sample[isample].Data(),Particle[ipart].Data(),Sign[icharge].Data(),ibin_data),Form("cDCA%d%s%s%sbin%d",index,sample[isample].Data(),Particle[ipart].Data(),Sign[icharge].Data(),ibin_data),1700,1500);
695 cDCA->SetMargin(0.1,0.02,0.1,0.02);
696 TLegend* Leg1 = new TLegend(0.6,0.6,0.85,0.85,"","NDC");
697 Leg1->SetFillStyle(kFALSE);
698 Leg1->SetLineColor(kWhite);
699 Leg1->SetBorderSize(0);
703 TH1F *hToFit =(TH1F*) ((TH1F*)hman_data->GetDCAHistogram1D(Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
705 TH1F *hToFit =(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
706 debug<<Particle[ipart].Data()<<" "<<Sign[icharge].Data()<<" "<<lowedge<<endl;
707 TH1F *hmc1=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
708 TH1F *hmc2=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaSecondaryWeakDecay%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
709 TH1F *hmc3=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaSecondaryMaterial%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
710 Double_t minentries=100;
711 debug<<hToFit->GetEntries()<<" "<<hmc1->GetEntries()<<" "<<hmc2->GetEntries()<<" "<<hmc3->GetEntries()<<endl;
712 debug<<((fitsettings&0x1)&&hmc2->GetEntries()<=minentries)<<" "<<((fitsettings&0x2)&&hmc3->GetEntries()<=minentries)<<endl;
713 cout<<hToFit->GetEntries()<<" "<<hmc1->GetEntries()<<" "<<hmc2->GetEntries()<<" "<<hmc3->GetEntries()<<endl;
714 cout<<((fitsettings&0x1)&&hmc2->GetEntries()<=minentries)<<" "<<((fitsettings&0x2)&&hmc3->GetEntries()<=minentries)<<endl;
715 if(hToFit->GetEntries()<=minentries || hmc1->GetEntries()<=minentries || ((fitsettings&0x1)&&hmc2->GetEntries()<=minentries) || ((fitsettings&0x2)&&hmc3->GetEntries()<=minentries))
722 Float_t corrforrebinning[4]={1.0,1.0,1.0,1.0};
723 Int_t binCutRange[]={hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1])};
726 if(hmc3->GetNbinsX()>300)
729 corrforrebinning[0]=hToFit->Integral(binCutRange[0],binCutRange[1]);
730 corrforrebinning[1]=hmc1->Integral(binCutRange[0],binCutRange[1]);
731 corrforrebinning[2]=hmc2->Integral(binCutRange[0],binCutRange[1]);
732 corrforrebinning[3]=hmc3->Integral(binCutRange[0],binCutRange[1]);
739 binCutRange[0]=hmc1->GetXaxis()->FindBin(CutRange[0]);
740 binCutRange[1]=hmc1->GetXaxis()->FindBin(CutRange[1]);
743 //after rebbing we lose resolution of the dca this correction also us to do obtain inside used dca
745 if(hToFit->Integral(binCutRange[0],binCutRange[1])>0.0)
746 corrforrebinning[0]=corrforrebinning[0]/hToFit->Integral(binCutRange[0],binCutRange[1]);
748 corrforrebinning[0]=1.0;
749 if(hmc1->Integral(binCutRange[0],binCutRange[1])>0.0)
750 corrforrebinning[1]=corrforrebinning[1]/hmc1->Integral(binCutRange[0],binCutRange[1]);
752 corrforrebinning[1]=1.0;
753 if(hmc2->Integral(binCutRange[0],binCutRange[1])>0.0)
754 corrforrebinning[2]=corrforrebinning[2]/hmc2->Integral(binCutRange[0],binCutRange[1]);
756 corrforrebinning[2]=1.0;
757 if(hmc3->Integral(binCutRange[0],binCutRange[1])>0.0)
758 corrforrebinning[3]=corrforrebinning[3]/hmc3->Integral(binCutRange[0],binCutRange[1]);
760 corrforrebinning[3]=1.0;
774 mc = new TObjArray(3); // MC histograms are put in this array
777 mc = new TObjArray(2);
786 TFractionFitter* fit = new TFractionFitter(hToFit,mc); // initialise
787 fit->Constrain(0,0.0,1.0); // constrain fraction 1 to be between 0 and 1
788 Double_t defaultStep = 0.01;
789 for (int iparm = 0; iparm < Npar; ++iparm)
791 TString name("frac"); name += iparm;
793 fit->GetFitter()->SetParameter(iparm,name.Data(),0.85,0.01,0.0,1.0);
795 fit->GetFitter()->SetParameter(iparm,name.Data(),0.15,0.01,0.0,1.0);
797 fit->GetFitter()->SetParameter(iparm,name.Data(),0.075,0.01,0.0,1.0);
801 fit->Constrain(1,0.0,1.0); // constrain fraction 1 to be between 0 and 1
803 fit->Constrain(1+(fitsettings&0x1),0.0,1.0); // constrain fraction 1 to be between 0 and 1
805 Int_t binFitRange[]={hmc1->GetXaxis()->FindBin(FitRange[0]),hmc1->GetXaxis()->FindBin(FitRange[1])};
806 fit->SetRangeX(binFitRange[0],binFitRange[1]);
807 hToFit->GetXaxis()->SetRange(binFitRange[0],binFitRange[1]);
808 // hToFit->SetTitle(Form("DCA distr - %s %s %s %lf",sample[isample].Data(),Particle[ipart].Data(),Sign[icharge].Data(),lowedge));
809 hToFit->SetTitle("");
810 Int_t status = fit->Fit(); // perform the fit
811 cout << "fit status: " << status << endl;
812 debug<<"fit status: " << status << endl;
816 Double_t v1=0.0,v2=0.0,v3=0.0;
817 Double_t ev1=0.0,ev2=0.0,ev3=0.0;
820 // check on fit status
821 TH1F* result = (TH1F*) fit->GetPlot();
822 TH1F* PrimMCPred=(TH1F*)fit->GetMCPrediction(0);
824 TH1F* secStMCPred=0X0;
828 //first method, use directly the fit result
829 fit->GetResult(0,v1,ev1);
833 fit->GetResult(1,v2,ev2);
834 secStMCPred=(TH1F*)fit->GetMCPrediction(1);
839 fit->GetResult(1+(fitsettings&0x1),v3,ev3);
840 secMCPred=(TH1F*)fit->GetMCPrediction(1+(fitsettings&0x1));
842 cov=fit->GetFitter()->GetCovarianceMatrixElement(1,2);
848 debug<<v1<<" "<<ev1<<" "<<v2<<" "<<ev2<<" "<<v3<<" "<<ev3<<" "<<" "<<cov<<endl;
850 // becuase dca cut range is not a fit range the results from TFractionFitter should be rescale
851 // This can be done in two ways or use input histograms or output histograms
852 // The difference between those two methods should be on the level of statistical error
853 // I use output histograms
855 // Method 1 input histo
857 Float_t normalizationdata=hToFit->Integral(hToFit->GetXaxis()->FindBin(CutRange[0]),hToFit->GetXaxis()->FindBin(CutRange[1]))/hToFit->Integral(binFitRange[0],binFitRange[1]);
858 normalizationdata*=corrforrebinning[0];
860 Float_t normalizationmc1=(hmc1->Integral(binCutRange[0],binCutRange[1])/hmc1->Integral(binFitRange[0],binFitRange[1]))/normalizationdata;
861 Float_t normalizationmc2=0.0;
863 normalizationmc2=(hmc2->Integral(binCutRange[0],binCutRange[1])/hmc2->Integral(binFitRange[0],binFitRange[1]))/normalizationdata;
864 Float_t normalizationmc3=0.0;
866 normalizationmc3=(hmc3->Integral(binCutRange[0],binCutRange[1])/hmc3->Integral(binFitRange[0],binFitRange[1]))/normalizationdata;
868 normalizationmc1*=corrforrebinning[1];
869 normalizationmc2*=corrforrebinning[2];
870 normalizationmc3*=corrforrebinning[3];
872 debug<<"After Nor"<<endl;
873 debug<<v1*normalizationmc1<<" "<<ev1*normalizationmc1<<" "<<v2*normalizationmc2<<" "<<ev2*normalizationmc2<<" "<<v3*normalizationmc3<<" "<<ev3*normalizationmc3<<" "<<endl;
874 debug<<1.0-v1*normalizationmc1<<" "<<ev1*normalizationmc1<<" "<<v2*normalizationmc2+v3*normalizationmc3<<" "<<TMath::Sqrt(ev2*ev2*normalizationmc2*normalizationmc2+ev3*ev3*normalizationmc3*normalizationmc3+cov*normalizationmc3*normalizationmc2)<<endl;
875 debug<<"addtional info"<<endl;
878 //Method 2 output histo
880 Float_t normalizationdata1=result->Integral(binCutRange[0],binCutRange[1])/result->Integral(binFitRange[0],binFitRange[1]);
883 // if the cut range is bigger the fit range we should calculate the normalization factor for data using the data histogram
884 // because result histogram has entries only in fits range
885 if(FitRange[0]>CutRange[0]||FitRange[1]<CutRange[1])
886 normalizationdata1=normalizationdata;
888 normalizationdata1*=corrforrebinning[0];
891 Float_t normalizationmc11=(PrimMCPred->Integral(binCutRange[0],binCutRange[1])/PrimMCPred->Integral(binFitRange[0],binFitRange[1]))/normalizationdata1;
892 Float_t normalizationmc21=0.0;
894 normalizationmc21=(secStMCPred->Integral(binCutRange[0],binCutRange[1])/secStMCPred->Integral(binFitRange[0],binFitRange[1]))/normalizationdata1;
895 Float_t normalizationmc31=0.0;
897 normalizationmc31=(secMCPred->Integral(binCutRange[0],binCutRange[1])/secMCPred->Integral(binFitRange[0],binFitRange[1]))/normalizationdata1;
899 normalizationmc11*=corrforrebinning[1];
900 normalizationmc21*=corrforrebinning[2];
901 normalizationmc31*=corrforrebinning[3];
903 debug<<"After Nor 2"<<endl;
904 debug<<v1*normalizationmc11<<" "<<ev1*normalizationmc11<<" "<<v2*normalizationmc21<<" "<<ev2*normalizationmc21<<" "<<v3*normalizationmc31<<" "<<ev3*normalizationmc31<<endl;
906 debug<<1.0-v1*normalizationmc11<<" "<<ev1*normalizationmc11<<" "<<v2*normalizationmc21+v3*normalizationmc31<<" "<<TMath::Sqrt(ev2*ev2*normalizationmc21*normalizationmc21+ev3*ev3*normalizationmc31*normalizationmc31+cov*normalizationmc31*normalizationmc21)<<endl;
909 hconWD[index+6*isample]->SetBinContent(ibin_data,v2*normalizationmc21);
910 hconWD[index+6*isample]->SetBinError(ibin_data,ev2*normalizationmc21);
911 hconMat[index+6*isample]->SetBinContent(ibin_data,v3*normalizationmc31);
912 hconMat[index+6*isample]->SetBinError(ibin_data,ev3*normalizationmc31);
913 hprimary[index+6*isample]->SetBinContent(ibin_data,v1*normalizationmc11);
914 hprimary[index+6*isample]->SetBinError(ibin_data,ev1*normalizationmc11);
915 hcon[index+6*isample]->SetBinContent(ibin_data,v2*normalizationmc21+v3*normalizationmc31);
916 hcon[index+6*isample]->SetBinError(ibin_data,TMath::Sqrt(ev2*ev2*normalizationmc21*normalizationmc21+ev3*ev3*normalizationmc31*normalizationmc31+cov*normalizationmc31*normalizationmc21));
921 result->Scale(1.0/result->Integral(result->GetXaxis()->FindBin(FitRange[0]),result->GetXaxis()->FindBin(FitRange[1])));
922 hToFit->Scale(1.0/hToFit->Integral(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1])));
923 PrimMCPred->Scale(v1/PrimMCPred->Integral(PrimMCPred->GetXaxis()->FindBin(FitRange[0]),PrimMCPred->GetXaxis()->FindBin(FitRange[1])));
925 hToFit->SetMinimum(0.0001);
926 hToFit->GetXaxis()->SetTitle("DCA_{xy} (cm)");
927 hToFit->GetYaxis()->SetTitle("N_{counts}/N_{counts}(-3cm;3cm)");
928 hToFit->GetYaxis()->SetTitleOffset(1.3);
929 hToFit->DrawClone("E1x0");
930 Leg1->AddEntry(hToFit,"data","p");
931 result->SetTitle("Fit result");
932 result->SetLineColor(kBlack);
933 Leg1->AddEntry(result,"fit result","l");
934 result->DrawClone("histsame");
936 PrimMCPred->SetLineColor(kGreen+2);
937 PrimMCPred->SetLineStyle(2);
938 PrimMCPred->SetLineWidth(3.0);
939 Leg1->AddEntry(PrimMCPred,"primaries","l");
940 PrimMCPred->DrawClone("histsame");
944 secStMCPred->Scale(v2/secStMCPred->Integral(secStMCPred->GetXaxis()->FindBin(FitRange[0]),secStMCPred->GetXaxis()->FindBin(FitRange[1])));
945 secStMCPred->SetLineColor(kRed);
946 secStMCPred->SetLineWidth(3.0);
948 secStMCPred->SetLineStyle(3);
949 Leg1->AddEntry(secStMCPred,"weak decays","l");
950 secStMCPred->DrawClone("histsame");
956 secMCPred->Scale(v3/secMCPred->Integral(secMCPred->GetXaxis()->FindBin(FitRange[0]),secMCPred->GetXaxis()->FindBin(FitRange[1])));
957 secMCPred->SetLineColor(kBlue);
958 secMCPred->SetLineWidth(3.0);
960 secMCPred->SetLineStyle(4);
961 Leg1->AddEntry(secMCPred,"material","l");
962 secMCPred->DrawClone("histsame");
968 hconWD[index+6*isample]->SetBinContent(ibin_data,0.0);
969 hconWD[index+6*isample]->SetBinError(ibin_data,0.0);
970 hconMat[index+6*isample]->SetBinContent(ibin_data,0.0);
971 hconMat[index+6*isample]->SetBinError(ibin_data,0.0);
972 hcon[index+6*isample]->SetBinContent(ibin_data,0.0);
973 hcon[index+6*isample]->SetBinError(ibin_data,0.0);
974 hprimary[index+6*isample]->SetBinContent(ibin_data,1.0);
975 hprimary[index+6*isample]->SetBinError(ibin_data,0.0);
978 TLatex* texttitle=new TLatex();
980 texttitle->SetTextSize(0.04);
981 texttitle->DrawLatex(0.12,0.92,Form("%s %.2f<#it{p}_{T} < %.2f (GeV/#it{c})",symboles[index].Data(),lowedge,binwidth+lowedge));
982 listofdcafits->Add(cDCA);
994 listofdcafits->Write("DCAfits",TObject::kSingleKey);
997 void RecomputeErrors(TH1* h)
999 for (int i=0; i<=h->GetXaxis()->GetNbins(); i++)
1001 cout<<h->GetBinContent(i)<<" "<<h->GetBinError(i)<<" error "<<TMath::Sqrt(h->GetBinContent(i))<<endl;
1002 h->SetBinError(i,TMath::Sqrt(h->GetBinContent(i)));
1006 void SetBintoOne(TH1* h)
1008 for (int i=0;i<=h->GetXaxis()->GetNbins(); i++)
1010 h->SetBinContent(i,1);
1011 h->SetBinError(i,0);
1016 void GetCorrectedSpectra(TH1F* corr,TH1F* raw,TH1F* eff, TH1F* con)
1018 for (int i=0;i<=corr->GetXaxis()->GetNbins(); i++)
1020 corr->SetBinContent(i,1);
1021 corr->SetBinError(i,0);
1028 corr->Multiply(raw);
1031 void GetCorrectedSpectraLeonardo(TH1F* spectra,TH1F* correction, TH1F* hprimaryData,TH1F* hprimaryMC)
1034 spectra->Multiply(correction);
1036 hprimaryData->Sumw2();
1037 spectra->Multiply(hprimaryData);
1038 hprimaryMC->Sumw2();
1039 spectra->Divide(hprimaryMC);
1042 void GFCorrection(TH1F **Spectra,Float_t tofpt,UInt_t options)
1044 if (options&kgeantflukaKaon)
1046 //Geant/Fluka Correction
1047 Printf("\nGF correction for Kaons");
1048 //Getting GF For Kaons in TPC
1049 TGraph *gGFCorrectionKaonPlus=new TGraph();
1050 gGFCorrectionKaonPlus->SetName("gGFCorrectionKaonPlus");
1051 gGFCorrectionKaonPlus->SetTitle("gGFCorrectionKaonPlus");
1052 TGraph *gGFCorrectionKaonMinus=new TGraph();
1053 gGFCorrectionKaonMinus->SetName("gGFCorrectionKaonMinus");
1054 gGFCorrectionKaonMinus->SetTitle("gGFCorrectionKaonMinus");
1055 TString fnameGeanFlukaK="GFCorrection/correctionForCrossSection.321.root";
1056 TFile *fGeanFlukaK= TFile::Open(fnameGeanFlukaK.Data());
1059 fnameGeanFlukaK="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/correctionForCrossSection.321.root";
1060 fGeanFlukaK= TFile::Open(fnameGeanFlukaK.Data());
1064 TH1F *hGeantFlukaKPos=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles");
1065 TH1F *hGeantFlukaKNeg=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles");
1066 //getting GF func for Kaons with TOF
1067 TF1 *fGFKPosTracking;
1068 fGFKPosTracking = TrackingEff_geantflukaCorrection(3,kPositive);
1069 TF1 *fGFKNegTracking;
1070 fGFKNegTracking = TrackingEff_geantflukaCorrection(3,kNegative);
1071 TF1 *fGFKPosMatching;
1072 fGFKPosMatching = TOFmatchMC_geantflukaCorrection(3,kPositive);
1073 TF1 *fGFKNegMatching;
1074 fGFKNegMatching = TOFmatchMC_geantflukaCorrection(3,kNegative);
1075 if(tcutsdata->GetUseTypeDependedTOFCut())
1076 tofpt=tcutsdata->GetPtTOFMatchingKaon();
1078 for(Int_t binK=0;binK<=Spectra[1]->GetNbinsX();binK++)
1080 if(Spectra[1]->GetBinCenter(binK)<tofpt)
1081 {//use TPC GeantFlukaCorrection
1082 Float_t FlukaCorrKPos=hGeantFlukaKPos->GetBinContent(hGeantFlukaKPos->FindBin(Spectra[1]->GetBinCenter(binK)));
1083 Float_t FlukaCorrKNeg=hGeantFlukaKNeg->GetBinContent(hGeantFlukaKNeg->FindBin(Spectra[4]->GetBinCenter(binK)));
1084 // Printf("TPC Geant/Fluka: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPos,FlukaCorrKNeg);
1085 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPos);
1086 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNeg);
1087 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPos);
1088 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNeg);
1089 gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),FlukaCorrKPos);
1090 gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),FlukaCorrKNeg);
1094 gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),0);
1095 gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),0);
1096 Float_t FlukaCorrKPosTracking=fGFKPosTracking->Eval(Spectra[1]->GetBinCenter(binK));
1097 Float_t FlukaCorrKNegTracking=fGFKNegTracking->Eval(Spectra[1]->GetBinCenter(binK));
1098 // Printf("TPC/TOF Geant/Fluka Tracking: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosTracking,FlukaCorrKNegTracking);
1099 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosTracking);
1100 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegTracking);
1101 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosTracking);
1102 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegTracking);
1103 Float_t FlukaCorrKPosMatching=fGFKPosMatching->Eval(Spectra[1]->GetBinCenter(binK));
1104 Float_t FlukaCorrKNegMatching=fGFKNegMatching->Eval(Spectra[1]->GetBinCenter(binK));
1105 // Printf("TPC/TOF Geant/Fluka Matching: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosMatching,FlukaCorrKNegMatching);
1106 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosMatching);
1107 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegMatching);
1108 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosMatching);
1109 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegMatching);
1113 if(!(options&kgeantflukaProton))
1115 //Geant Fluka for P in TPC
1116 Printf("\nGF correction for Protons");
1117 const Int_t kNCharge=2;
1120 TString fnameGFProtons= "GFCorrection/correctionForCrossSection.root";
1121 TFile* fGFProtons = TFile::Open(fnameGFProtons.Data());
1124 fnameGFProtons="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/correctionForCrossSection.root";
1125 fGFProtons = TFile::Open(fnameGFProtons.Data());
1132 TH2D * hCorrFluka[kNCharge];
1133 TH2D * hCorrFluka[2];
1134 hCorrFluka[kPos] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionProtons");
1135 hCorrFluka[kNeg] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionAntiProtons");
1136 //getting GF func for Kaons with TPCTOF
1137 TF1 *fGFpPosTracking;
1138 fGFpPosTracking = TrackingEff_geantflukaCorrection(4,kPositive);
1139 TF1 *fGFpNegTracking;
1140 fGFpNegTracking = TrackingEff_geantflukaCorrection(4,kNegative);
1141 TF1 *fGFpPosMatching;
1142 fGFpPosMatching = TOFmatchMC_geantflukaCorrection(4,kPositive);
1143 TF1 *fGFpNegMatching;
1144 fGFpNegMatching = TOFmatchMC_geantflukaCorrection(4,kNegative);
1147 Int_t nbins = Spectra[2]->GetNbinsX();
1148 if(tcutsdata->GetUseTypeDependedTOFCut())
1149 tofpt=tcutsdata->GetPtTOFMatchingProton();
1151 for(Int_t ibin = 0; ibin < nbins; ibin++)
1153 if(Spectra[2]->GetBinCenter(ibin)<tofpt)
1154 {//use TPC GeantFlukaCorrection
1155 for(Int_t icharge = 0; icharge < kNCharge; icharge++)
1157 Int_t nbinsy=hCorrFluka[icharge]->GetNbinsY();
1158 Float_t pt = Spectra[2]->GetBinCenter(ibin);
1159 Float_t minPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(1);
1160 Float_t maxPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1);
1161 if (pt < minPtCorrection)
1162 pt = minPtCorrection+0.0001;
1163 if (pt > maxPtCorrection)
1164 pt = maxPtCorrection;
1165 Float_t correction = hCorrFluka[icharge]->GetBinContent(1,hCorrFluka[icharge]->GetYaxis()->FindBin(pt));
1167 if (correction > 0.0)
1168 {// If the bin is empty this is a 0
1169 // cout<<icharge<<" "<<ibin<<" "<<correction<<endl;
1170 Spectra[icharge*3+2]->SetBinContent(ibin,Spectra[icharge*3+2]->GetBinContent(ibin)*correction);
1171 Spectra[icharge*3+2]->SetBinError(ibin,Spectra[icharge*3+2]->GetBinError(ibin)*correction);
1173 else if (Spectra[icharge*3+2]->GetBinContent(ibin) > 0.0)
1175 // If we are skipping a non-empty bin, we notify the user
1176 cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for " <<"protons Pt:"<< pt<< endl;
1177 cout << " Bin content: " << Spectra[icharge*3+2]->GetBinContent(ibin) << endl;
1183 Float_t FlukaCorrpPosTracking=fGFpPosTracking->Eval(Spectra[2]->GetBinCenter(ibin));
1184 Float_t FlukaCorrpNegTracking=fGFpNegTracking->Eval(Spectra[2]->GetBinCenter(ibin));
1185 // Printf("TPC/TOF Geant/Fluka Tracking: pt:%f Pos:%f Neg:%f",Spectra[2]->GetBinCenter(ibin),FlukaCorrpPosTracking,FlukaCorrpNegTracking);
1186 Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*FlukaCorrpPosTracking);
1187 Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*FlukaCorrpNegTracking);
1188 Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError(ibin)*FlukaCorrpPosTracking);
1189 Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError(ibin)*FlukaCorrpNegTracking);
1190 Float_t FlukaCorrpPosMatching=fGFpPosMatching->Eval(Spectra[2]->GetBinCenter(ibin));
1191 Float_t FlukaCorrpNegMatching=fGFpNegMatching->Eval(Spectra[2]->GetBinCenter(ibin));
1192 // Printf("TPC/TOF Geant/Fluka Matching: pt:%f Pos:%f Neg:%f",Spectra[2]->GetBinCenter(ibin),FlukaCorrpPosMatching,FlukaCorrpNegMatching);
1193 Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*FlukaCorrpPosMatching);
1194 Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*FlukaCorrpNegMatching);
1195 Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError(ibin)*FlukaCorrpPosMatching);
1196 Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError(ibin)*FlukaCorrpNegMatching);
1199 fGeanFlukaK->Close();
1206 TrackingEff_geantflukaCorrection(Int_t ipart, Int_t icharge)
1209 if (ipart == 3 && icharge == kNegative) {
1210 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.);
1213 else if (ipart == 4 && icharge == kNegative) {
1214 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.);
1217 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionNull(x)", 0., 5.);
1223 TrackingPtGeantFlukaCorrectionNull(Double_t pTmc)
1229 TrackingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
1231 return (1 - 0.129758 *TMath::Exp(-pTmc*0.679612));
1235 TrackingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
1237 return TMath::Min((0.972865 + 0.0117093*pTmc), 1.);
1239 ///////////////////////////////////////////
1241 TOFmatchMC_geantflukaCorrection(Int_t ipart, Int_t icharge)
1244 if (ipart == 3 && icharge == kNegative) {
1245 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.);
1248 else if (ipart == 4 && icharge == kNegative) {
1249 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.);
1252 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionNull(x)", 0., 5.);
1259 MatchingPtGeantFlukaCorrectionNull(Double_t pTmc)
1265 MatchingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
1267 Float_t ptTPCoutP =pTmc*(1-6.81059e-01*TMath::Exp(-pTmc*4.20094));
1268 return (TMath::Power(1 - 0.129758*TMath::Exp(-ptTPCoutP*0.679612),0.07162/0.03471));
1272 MatchingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
1274 Float_t ptTPCoutK=pTmc*(1- 3.37297e-03/pTmc/pTmc - 3.26544e-03/pTmc);
1275 return TMath::Min((TMath::Power(0.972865 + 0.0117093*ptTPCoutK,0.07162/0.03471)), 1.);
1277 void MatchingTOFEff(TH1F** Spectra, TList* list=0x0)
1279 if(TOFMatchingScalling[0]<0.0&&TOFMatchingScalling[1]<0.0)
1281 TH1F *hMatcEffPos_data=(TH1F*)tcutsdata->GetHistoNMatchedPos();
1282 hMatcEffPos_data->Sumw2();
1283 //hMatcEffPos_data->Divide((TH1F*)tcutsdata->GetHistoNSelectedPos());
1284 hMatcEffPos_data->Divide(hMatcEffPos_data,(TH1F*)tcutsdata->GetHistoNSelectedPos(),1,1,"B");
1285 hMatcEffPos_data->SetTitle("Matching Eff Pos - data");
1286 TH1F *hMatcEffNeg_data=(TH1F*)tcutsdata->GetHistoNMatchedNeg();
1287 hMatcEffNeg_data->Sumw2();
1288 //hMatcEffNeg_data->Divide((TH1F*)tcutsdata->GetHistoNSelectedNeg());
1289 hMatcEffNeg_data->Divide(hMatcEffNeg_data,(TH1F*)tcutsdata->GetHistoNSelectedNeg(),1,1,"B");
1290 hMatcEffNeg_data->SetTitle("Matching Eff Neg - data");
1291 TH1F *hMatcEffPos_mc=(TH1F*)tcutsmc->GetHistoNMatchedPos();
1292 hMatcEffPos_mc->Sumw2();
1293 //hMatcEffPos_mc->Divide((TH1F*)tcutsmc->GetHistoNSelectedPos());
1294 hMatcEffPos_mc->Divide(hMatcEffPos_mc,(TH1F*)tcutsmc->GetHistoNSelectedPos(),1,1,"B");
1295 hMatcEffPos_mc->SetTitle("Matching Eff Pos - mc");
1296 TH1F *hMatcEffNeg_mc=(TH1F*)tcutsmc->GetHistoNMatchedNeg();
1297 hMatcEffNeg_mc->Sumw2();
1298 //hMatcEffNeg_mc->Divide((TH1F*)tcutsmc->GetHistoNSelectedNeg());
1299 hMatcEffNeg_mc->Divide(hMatcEffNeg_mc,(TH1F*)tcutsmc->GetHistoNSelectedNeg(),1,1,"B");
1300 hMatcEffNeg_mc->SetTitle("Matching Eff Neg - mc");
1303 hMatcEffPos_data->Divide(hMatcEffPos_mc);
1304 hMatcEffNeg_data->Divide(hMatcEffNeg_mc);
1305 hMatcEffPos_data->SetName("MatchingTOFPos");
1306 hMatcEffNeg_data->SetName("MatchingTOFNeg");
1309 TF1 *pol0MatchPos_data=new TF1("pol0MatchPos_data","pol0",.6,5);
1310 hMatcEffPos_data->Fit("pol0MatchPos_data","MNR");
1311 TF1 *pol0MatchNeg_data=new TF1("pol0MatchNeg_data","pol0",.6,5);
1312 hMatcEffNeg_data->Fit("pol0MatchNeg_data","MNR");
1314 list->Add(hMatcEffPos_data);
1315 list->Add(hMatcEffNeg_data);
1318 TOFMatchingScalling[0]=pol0MatchPos_data->GetParameter(0);
1319 TOFMatchingScalling[1]=pol0MatchNeg_data->GetParameter(0);
1322 //Correction spectra for matching efficiency
1323 //For the moment I'm using the inclusive correction
1324 for(Int_t ipart=0;ipart<3;ipart++)
1327 for(Int_t ibin=1;ibin<Spectra[ipart]->GetNbinsX();ibin++)
1329 Float_t ptspectra=Spectra[ipart]->GetBinCenter(ibin);
1330 if(ptspectra<tcutsdata->GetPtTOFMatching())
1332 //Spectra[ipart]->SetBinContent(ibin,( Spectra[ipart]->GetBinContent(ibin)/hMatcEffPos_data->GetBinContent(hMatcEffPos_data->FindBin(ptspectra))));
1333 //Spectra[ipart+3]->SetBinContent(ibin,( Spectra[ipart+3]->GetBinContent(ibin)/hMatcEffNeg_data->GetBinContent(hMatcEffNeg_data->FindBin(ptspectra))));
1334 Spectra[ipart]->SetBinContent(ibin,( Spectra[ipart]->GetBinContent(ibin)/TOFMatchingScalling[0]));
1335 Spectra[ipart+3]->SetBinContent(ibin,( Spectra[ipart+3]->GetBinContent(ibin)/TOFMatchingScalling[1]));
1339 Double_t eta2y(Double_t pt, Double_t mass, Double_t eta)
1341 Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
1342 return TMath::ASinH(pt / mt * TMath::SinH(eta));
1345 TH1* GetSumAllCh(TH1F** spectra, Double_t* mass,Double_t etacut)
1347 TH1F* allch=(((TH1F*))spectra[0]->Clone("allCh"));
1349 for (int i=0;i<6;i++)
1351 Double_t masstmp=mass[i%3];
1352 for (int j=1;j<spectra[i]->GetXaxis()->GetNbins();j++)
1354 Float_t value=spectra[i]->GetBinContent(j);
1355 Float_t error=spectra[i]->GetBinError(j);
1358 Float_t pt=spectra[i]->GetXaxis()->GetBinCenter(j);
1359 Float_t mt2=masstmp*masstmp+pt*pt;
1360 Float_t mt=TMath::Sqrt(mt2);
1361 Float_t maxy=eta2y(pt,masstmp,etacut/2.0);
1362 Float_t conver=maxy*(TMath::Sqrt(1-masstmp*masstmp/(mt2*TMath::CosH(maxy)*TMath::CosH(maxy)))+TMath::Sqrt(1-masstmp*masstmp/(mt2*TMath::CosH(0.0)*TMath::CosH(0.0))));
1363 conver=conver/etacut;
1364 //cout<<maxy<<" "<<conver<<" "<<masstmp<<""<<spectra[i]->GetName()<<endl;
1365 Float_t bincontent=allch->GetBinContent(j);
1366 Float_t binerror=allch->GetBinError(j);
1367 bincontent+=conver*value;
1368 binerror=TMath::Sqrt(binerror*binerror+conver*conver*error*error);
1369 allch->SetBinContent(j,bincontent);
1370 allch->SetBinError(j,binerror);
1375 Divideby2pipt(allch);
1376 allch->SetTitle("N_{ch};p_{T} (GeV/c);1/(2 #pi p_{T})dN/p_{T} |#eta|<0.8");
1380 void Divideby2pipt(TH1* hist)
1383 for (int i=1;i<hist->GetXaxis()->GetNbins();i++)
1385 Float_t value=hist->GetBinContent(i);
1386 Float_t error=hist->GetBinError(i);
1387 Float_t pt=hist->GetXaxis()->GetBinCenter(i);
1388 hist->SetBinContent(i,value/(2.0*TMath::Pi()*pt));
1389 hist->SetBinError(i,error/(2.0*TMath::Pi()*pt));
1394 Short_t DCAfitsettings (Float_t pt, Int_t type)
1397 if (pt<maxptformaterial[type]&&pt>minptformaterial[type])
1399 if (pt<maxptforWD[type]&&pt>minptforWD[type])
1405 Float_t Normaliztionwithbin0integrals(UInt_t options)
1408 TH1F* bin0mcRec=(TH1F*)ecutsmc->GetHistoVtxGenerated()->Clone("Bin0_rec");
1409 TH1F* bin0mcMC=(TH1F*)ecutsmc->GetHistoVtxGenerated()->Clone("Bin0_MC");
1411 TH1F* vertexmc=ecutsmc->GetHistoVtxAftSelwithoutZvertexCut();
1412 TH1F* vertexmcMCz=ecutsmc->GetHistoVtxAftSelwithoutZvertexCutusingMCz();
1413 TH1F* vertexdata=ecutsdata->GetHistoVtxAftSelwithoutZvertexCut();
1415 TH1I* histodata=ecutsdata->GetHistoCuts();
1416 TH1I* histomc=ecutsmc->GetHistoCuts();
1418 Float_t dataevents=(Float_t)histodata->GetBinContent(3);
1419 //cout<<histodata->GetBinContent(2)<<endl;
1420 Float_t databin0events=((Float_t)histodata->GetBinContent(2))-((Float_t)histodata->GetBinContent(4));
1425 bin0mcRec->Add(vertexmc,-1);
1426 bin0mcMC->Add(vertexmcMCz,-1);
1428 bin0mcRec->Divide(vertexmc);
1429 bin0mcMC->Divide(vertexmcMCz);
1431 bin0mcRec->Multiply(vertexdata);
1432 bin0mcMC->Multiply(vertexdata);
1434 Float_t bin0mcRecN=0.0;
1435 Float_t bin0mcMCN=0.0;
1437 for (int i=0;i<=bin0mcRec->GetXaxis()->GetNbins();i++)
1439 bin0mcRecN+=bin0mcRec->GetBinContent(i);
1440 bin0mcMCN+=bin0mcMC->GetBinContent(i);
1443 bin0mcRec->Scale(databin0events/bin0mcRecN);
1444 bin0mcMC->Scale(databin0events/bin0mcMCN);
1446 Int_t binmin=bin0mcRec->GetXaxis()->FindBin(-10);
1447 Int_t binmax=bin0mcRec->GetXaxis()->FindBin(10)-1;
1448 cout<< bin0mcRec->GetXaxis()->GetBinLowEdge(binmin)<<" "<<bin0mcRec->GetXaxis()->GetBinUpEdge(binmax)<<endl;
1449 cout<<bin0mcRecN<<" "<<bin0mcMCN<<" "<<databin0events<<endl;
1450 cout<<dataevents<<" normalization "<<dataevents+bin0mcRec->Integral(binmin,binmax)<<" "<<dataevents+bin0mcMC->Integral(binmin,binmax)<<endl;
1451 cout<<histodata->GetBinContent(2)<<" "<<histodata->GetBinContent(4)<<endl;
1452 if ((options&knormalizationwithbin0integralsdata)==knormalizationwithbin0integralsdata)
1453 return dataevents+bin0mcRec->Integral(binmin,binmax);
1454 else if ((options&kknormalizationwithbin0integralsMC)==knormalizationwithbin0integralsMC)
1455 return dataevents+bin0mcMC->Integral(binmin,binmax) ;
1461 Bool_t ReadConfigFile(TString configfile)
1463 ifstream infile(configfile.Data());
1464 if(infile.is_open()==false)
1466 TString namesofSetting[]={"CutRangeMin","CutRangeMax","FitRangeMin","FitRangeMax","MinMatPionPlus","MaxMatPionPlus","MinMatKaonPlus","MaxMatKaonPlus","MinMatProtonPlus","MaxMatProtonPlus","MinMatPionMinus","MaxMatPionMinus","MinMatKaonMinus","MaxMatKaonMinus","MinMatProtonMinus","MaxMatProtonMinus","MinWDPionPlus","MaxWDPionPlus","MinWDKaonPlus","MaxWDKaonPlus","MinWDProtonPlus","MaxWDProtonPlus","MinWDPionMinus","MaxWDPionMinus","MinWDKaonMinus","MaxWDKaonMinus","MinWDProtonMinus","MaxWDProtonMinus","MaxContaminationPIDMC","MinPions","MaxPions","MinKaons","MaxKaons","MinProtons","MaxProtons","TOFPIDsignalmatchPion","TOFPIDsignalmatchKaon","TOFPIDsignalmatchProton","NamefileEff","NameFilePIDcon","NameFileSeccon"};
1469 while (infile.eof()==false)
1472 while (buffer[0]=='#'&&infile.eof()==false)
1473 infile.getline(buffer,256);
1474 TString tmpstring(buffer);
1476 if(tmpstring.Contains(namesofSetting[0]))
1477 CutRange[0]=(tmpstring.Remove(0,namesofSetting[0].Length()+1)).Atof();
1478 else if (tmpstring.Contains(namesofSetting[1]))
1479 CutRange[1]=(tmpstring.Remove(0,namesofSetting[1].Length()+1)).Atof();
1480 else if (tmpstring.Contains(namesofSetting[2]))
1481 FitRange[0]=(tmpstring.Remove(0,namesofSetting[2].Length()+1)).Atof();
1482 else if (tmpstring.Contains(namesofSetting[3]))
1483 FitRange[1]=(tmpstring.Remove(0,namesofSetting[3].Length()+1)).Atof();
1484 else if (tmpstring.Contains(namesofSetting[4]))
1485 minptformaterial[0]=(tmpstring.Remove(0,namesofSetting[4].Length()+1)).Atof();
1486 else if (tmpstring.Contains(namesofSetting[5]))
1487 maxptformaterial[0]=(tmpstring.Remove(0,namesofSetting[5].Length()+1)).Atof();
1488 else if (tmpstring.Contains(namesofSetting[6]))
1489 minptformaterial[1]=(tmpstring.Remove(0,namesofSetting[6].Length()+1)).Atof();
1490 else if (tmpstring.Contains(namesofSetting[7]))
1491 maxptformaterial[1]=(tmpstring.Remove(0,namesofSetting[7].Length()+1)).Atof();
1492 else if (tmpstring.Contains(namesofSetting[8]))
1493 minptformaterial[2]=(tmpstring.Remove(0,namesofSetting[8].Length()+1)).Atof();
1494 else if (tmpstring.Contains(namesofSetting[9]))
1495 maxptformaterial[2]=(tmpstring.Remove(0,namesofSetting[9].Length()+1)).Atof();
1496 else if (tmpstring.Contains(namesofSetting[10]))
1497 minptformaterial[3]=(tmpstring.Remove(0,namesofSetting[10].Length()+1)).Atof();
1498 else if (tmpstring.Contains(namesofSetting[11]))
1499 maxptformaterial[3]=(tmpstring.Remove(0,namesofSetting[11].Length()+1)).Atof();
1500 else if (tmpstring.Contains(namesofSetting[12]))
1501 minptformaterial[4]=(tmpstring.Remove(0,namesofSetting[12].Length()+1)).Atof();
1502 else if (tmpstring.Contains(namesofSetting[13]))
1503 maxptformaterial[4]=(tmpstring.Remove(0,namesofSetting[13].Length()+1)).Atof();
1504 else if (tmpstring.Contains(namesofSetting[14]))
1505 minptformaterial[5]=(tmpstring.Remove(0,namesofSetting[14].Length()+1)).Atof();
1506 else if (tmpstring.Contains(namesofSetting[15]))
1507 maxptformaterial[5]=(tmpstring.Remove(0,namesofSetting[15].Length()+1)).Atof();
1508 else if (tmpstring.Contains(namesofSetting[16]))
1509 minptforWD[0]=(tmpstring.Remove(0,namesofSetting[16].Length()+1)).Atof();
1510 else if (tmpstring.Contains(namesofSetting[17]))
1511 maxptforWD[0]=(tmpstring.Remove(0,namesofSetting[17].Length()+1)).Atof();
1512 else if (tmpstring.Contains(namesofSetting[18]))
1513 minptforWD[1]=(tmpstring.Remove(0,namesofSetting[18].Length()+1)).Atof();
1514 else if (tmpstring.Contains(namesofSetting[19]))
1515 maxptforWD[1]=(tmpstring.Remove(0,namesofSetting[19].Length()+1)).Atof();
1516 else if (tmpstring.Contains(namesofSetting[20]))
1517 minptforWD[2]=(tmpstring.Remove(0,namesofSetting[20].Length()+1)).Atof();
1518 else if (tmpstring.Contains(namesofSetting[21]))
1519 maxptforWD[2]=(tmpstring.Remove(0,namesofSetting[21].Length()+1)).Atof();
1520 else if (tmpstring.Contains(namesofSetting[22]))
1521 minptforWD[3]=(tmpstring.Remove(0,namesofSetting[22].Length()+1)).Atof();
1522 else if (tmpstring.Contains(namesofSetting[23]))
1523 maxptforWD[3]=(tmpstring.Remove(0,namesofSetting[23].Length()+1)).Atof();
1524 else if (tmpstring.Contains(namesofSetting[24]))
1525 minptforWD[4]=(tmpstring.Remove(0,namesofSetting[24].Length()+1)).Atof();
1526 else if (tmpstring.Contains(namesofSetting[25]))
1527 maxptforWD[4]=(tmpstring.Remove(0,namesofSetting[25].Length()+1)).Atof();
1528 else if (tmpstring.Contains(namesofSetting[26]))
1529 minptforWD[5]=(tmpstring.Remove(0,namesofSetting[26].Length()+1)).Atof();
1530 else if (tmpstring.Contains(namesofSetting[27]))
1531 maxptforWD[5]=(tmpstring.Remove(0,namesofSetting[27].Length()+1)).Atof();
1532 else if (tmpstring.Contains(namesofSetting[28]))
1533 fMaxContaminationPIDMC=(tmpstring.Remove(0,namesofSetting[28].Length()+1)).Atof();
1534 else if (tmpstring.Contains(namesofSetting[29]))
1535 minRanges[0]=(tmpstring.Remove(0,namesofSetting[29].Length()+1)).Atof();
1536 else if (tmpstring.Contains(namesofSetting[30]))
1537 maxRanges[0]=(tmpstring.Remove(0,namesofSetting[30].Length()+1)).Atof();
1538 else if (tmpstring.Contains(namesofSetting[31]))
1539 minRanges[1]=(tmpstring.Remove(0,namesofSetting[31].Length()+1)).Atof();
1540 else if (tmpstring.Contains(namesofSetting[32]))
1541 maxRanges[1]=(tmpstring.Remove(0,namesofSetting[32].Length()+1)).Atof();
1542 else if (tmpstring.Contains(namesofSetting[33]))
1543 minRanges[2]=(tmpstring.Remove(0,namesofSetting[33].Length()+1)).Atof();
1544 else if (tmpstring.Contains(namesofSetting[34]))
1545 maxRanges[2]=(tmpstring.Remove(0,namesofSetting[34].Length()+1)).Atof();
1546 else if (tmpstring.Contains(namesofSetting[35]))
1547 TOFPIDsignalmatching[0]=(tmpstring.Remove(0,namesofSetting[35].Length()+1)).Atof();
1548 else if (tmpstring.Contains(namesofSetting[36]))
1549 TOFPIDsignalmatching[1]=(tmpstring.Remove(0,namesofSetting[36].Length()+1)).Atof();
1550 else if (tmpstring.Contains(namesofSetting[37]))
1551 TOFPIDsignalmatching[2]=(tmpstring.Remove(0,namesofSetting[37].Length()+1)).Atof();
1552 else if (tmpstring.Contains(namesofSetting[38]))
1553 filenames[0]=tmpstring.Remove(0,namesofSetting[38].Length()+1);
1554 else if (tmpstring.Contains(namesofSetting[39]))
1555 filenames[1]=tmpstring.Remove(0,namesofSetting[39].Length()+1);
1556 else if (tmpstring.Contains(namesofSetting[40]))
1557 filenames[2]=tmpstring.Remove(0,namesofSetting[40].Length()+1);
1562 // Double_t minRanges[3]={0.3,0.3,0.45};
1563 //Double_t maxRanges[3]={1.5,1.2,2.2};
1564 //Double_t fMaxContaminationPIDMC=0.2;
1569 for(int i=0;i<6;i++)
1570 cout<<minptformaterial[i]<<" "<<maxptformaterial[i]<<" "<<minptforWD[i]<<" "<<maxptforWD[i]<<endl;
1571 cout<<FitRange[0]<<" "<<FitRange[1]<<" "<<CutRange[0]<<CutRange[1]<<endl;
1572 if(FitRange[0]>=FitRange[1])
1577 if(CutRange[0]>=CutRange[1])
1582 for(int i=0;i<6;i++)
1584 if((minptformaterial[i]>maxptformaterial[i]&&minptformaterial[i]>0.0)||minptformaterial[i]<0.0||maxptformaterial[i]<0.0)
1589 if((minptforWD[i]>maxptforWD[i]&&minptforWD[i]>0.0)||minptforWD[i]<0.0||maxptforWD[i]<0.0)
1595 for(int i=0;i<3;i++)
1596 if(minRanges[i]>maxRanges[i])
1603 void SubHistWithFullCorr(TH1F* h1, TH1F* h2, Float_t factor1=1.0, Float_t factor2=1.0)
1605 if(h1->GetNbinsX()!=h2->GetNbinsX())
1607 for (int i=0;i<=h1->GetNbinsX();i++)
1609 Float_t tmpvalue=factor1*h1->GetBinContent(i)-factor2*h2->GetBinContent(i);
1610 Float_t tmperror=TMath::Abs(factor1*factor1*h1->GetBinError(i)*h1->GetBinError(i)-factor2*factor2*h2->GetBinError(i)*h2->GetBinError(i));
1611 h1->SetBinContent(i,tmpvalue);
1612 h1->SetBinError(i,TMath::Sqrt(tmperror));
1617 void TOFMatchingForNch(TH1* h)
1619 if(TOFMatchingScalling[0]>0.0&&TOFMatchingScalling[1]>0.0)
1621 Float_t factor=0.5*TOFMatchingScalling[0]+0.5*TOFMatchingScalling[1];
1622 for(Int_t ibin=1;ibin<h->GetNbinsX();ibin++)
1624 Float_t ptspectra=h->GetBinCenter(ibin);
1625 if(ptspectra<tcutsdata->GetPtTOFMatching())
1627 h->SetBinContent(ibin,(h->GetBinContent(ibin)/factor));
1636 void TOFPIDsignalmatchingApply(TH1* h, Float_t factor)
1640 for(Int_t ibin=1;ibin<h->GetNbinsX();ibin++)
1642 Float_t ptspectra=h->GetBinCenter(ibin);
1643 if(ptspectra<tcutsdata->GetPtTOFMatching())
1645 h->SetBinContent(ibin,(h->GetBinContent(ibin)*factor));
1650 void CalculateDoubleCounts(TH1* doubleconunts,TH1** rawspectra,Int_t ipar, Bool_t dataflag)
1654 tmphist=(TH2F*)managerdata->GetGenMulvsRawMulHistogram("hHistDoubleCounts");
1656 tmphist=(TH2F*)managermc->GetGenMulvsRawMulHistogram("hHistDoubleCounts");
1660 TH1F* tmphist1=(TH1F*)rawspectra[ipar]->Clone("test");
1661 tmphist1->Add(rawspectra[ipar+3]);
1662 doubleconunts->Add(tmphist->ProjectionX(doubleconunts->GetName(),1,1));
1664 doubleconunts->Add(tmphist->ProjectionX("pi+k",2,2));
1666 doubleconunts->Add(tmphist->ProjectionX("pi+p",3,3));
1668 doubleconunts->Add(tmphist->ProjectionX("k+p",4,4));
1669 doubleconunts->Divide(doubleconunts,tmphist1,1,1,"B");
1676 void CopyCorrectionFromFile(TString filename,TString correctionname,TH1F** corrtab)
1678 TFile* ftmp=TFile::Open(filename.Data());
1680 TList* ltmp=(TList*)ftmp->Get("output");
1682 for(Int_t icharge=0;icharge<2;icharge++)
1684 for(Int_t ipart=0;ipart<3;ipart++)
1686 Int_t index=ipart+3*icharge;
1687 tmp.Form("%s%s%s",correctionname.Data(),Particle[ipart].Data(),Sign[icharge].Data());
1688 TH1* histtmp=(TH1*)ltmp->FindObject(tmp.Data());
1691 for(Int_t ibin=1;ibin<=histtmp->GetNbinsX();ibin++)
1693 Float_t content=histtmp->GetBinContent(ibin);
1694 Float_t error=histtmp->GetBinError(ibin);
1695 Int_t bin=corrtab[index]->FindBin(histtmp->GetBinCenter(ibin));
1696 corrtab[index]->SetBinContent(bin,content);
1697 corrtab[index]->SetBinError(bin,error);