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;
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
58 Bool_t OpenFile(TString dirname, TString outputname, Bool_t mcflag,Bool_t mcasdata=false);
59 void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TString outnamemc="",TString configfile="" )
61 gStyle->SetOptStat(0);
62 TH1::AddDirectory(kFALSE);
63 gSystem->Load("libCore.so");
64 gSystem->Load("libPhysics.so");
65 gSystem->Load("libTree");
66 gSystem->Load("libMatrix");
67 gSystem->Load("libSTEERBase");
68 gSystem->Load("libESD");
69 gSystem->Load("libAOD");
70 gSystem->Load("libANALYSIS");
71 gSystem->Load("libOADB");
72 gSystem->Load("libANALYSISalice");
73 gSystem->Load("libTENDER");
74 gSystem->Load("libCORRFW");
75 gSystem->Load("libPWGTools");
76 gSystem->Load("libPWGLFspectra");
78 gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/QAPlotsBoth.C");
80 mass[0] = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
81 mass[1] = TDatabasePDG::Instance()->GetParticle("K+")->Mass();
82 mass[2] = TDatabasePDG::Instance()->GetParticle("proton")->Mass();
84 TFormula* dcacutxy=0x0;
85 if(!(options&kdonotusedcacuts))
88 AliESDtrackCuts* esdtrackcuts= AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE);
89 TString formulastring(esdtrackcuts->GetMaxDCAToVertexXYPtDep());
90 formulastring.ReplaceAll("pt","x");
91 dcacutxy=new TFormula("dcacutxy",formulastring.Data());
93 if(options&kuserangeonfigfile)
94 if(!ReadConfigFile(configfile))
96 TList* lout=new TList();
99 TString indirname=Form("/output/train%s",outdate.Data());
100 //TString indirname("/output/train24072012");
101 if(outnamemc.Length()==0)
102 outnamemc=outnamedata;
103 cout<<indirname.Data()<<" "<<outnamemc.Data()<<endl;
107 OpenFile(indirname,outnamemc,true);
108 OpenFile(indirname,outnamedata,false,((Bool_t)(options&kmcisusedasdata)));
109 if(!managermc||!managerdata)
111 cout<<managermc<<" "<<managerdata<<endl;
114 TH1F* rawspectradata[6];
115 TH1F* rawspectramc[6];
128 TH1F* contMatfit[12];
129 TH1F* primaryfit[12];
134 TH1F* spectraLeonardo[6];
136 TH1F* corrLeonardo[6];
137 TH1F* doubleconuntsdata[3];
138 TH1F* doubleconuntsMC[3];
139 //GetSpectra(managerdata,rawspectradata,true);
140 //GetSpectra(managermc,rawspectramc,true,true);
142 GetPtHistFromPtDCAhisto("hHistPtRecSigma","SpectraMC",managermc,rawspectramc,dcacutxy);
143 GetPtHistFromPtDCAhisto("hHistPtRecSigma","SpectraDATA",managerdata,rawspectradata,dcacutxy);
144 GetPtHistFromPtDCAhisto("hHistPtRecTruePrimary","eff",managermc,eff,dcacutxy);
145 GetPtHistFromPtDCAhisto("hHistPtRecTrue","conPID",managermc,contPID,dcacutxy);
146 GetPtHistFromPtDCAhisto("hHistPtRecSigmaSecondaryWeakDecay","conWD",managermc,contWD,dcacutxy);
147 GetPtHistFromPtDCAhisto("hHistPtRecSigmaSecondaryMaterial","conMat",managermc,contMat,dcacutxy);
148 GetPtHistFromPtDCAhisto("hHistPtRecSigmaPrimary","conPIDprimary",managermc,contPIDpri,dcacutxy);
151 Double_t neventsmcall = 1 ; //if loop over MC is done after or befor events cuts this will be changed
152 Double_t neventsdata = 1;
153 Double_t neventsmc = 1;
155 //Normaliztion of MCtruth depends if the loop was done after of before ESD event cuts.
156 //In currect code this cannot be check on the level of macro.
157 //If the loop was done before MC should be done to all processed events (NumberOfProcessedEvents())
158 //If loop was done after MC should be normalized to all accepted events (NumberOfEvents())
159 // The option one will be alaways use.
161 neventsmcall= ecutsmc->NumberOfProcessedEvents();
162 if(options&knormalizationtoeventspassingPhySel)
164 //neventsmcall= ecutsmc->NumberOfProcessedEvents();
165 neventsdata=ecutsdata->NumberOfPhysSelEvents();
166 neventsmc=ecutsmc->NumberOfPhysSelEvents();
168 else if ((options&knormalizationwithbin0integralsdata)||(options&knormalizationwithbin0integralsMC))
170 neventsdata=Normaliztionwithbin0integrals(options);
171 neventsmc=ecutsmc->NumberOfPhysSelEvents();
175 neventsdata=ecutsdata->NumberOfEvents(); //number of accepted events
176 neventsmc=ecutsmc->NumberOfEvents();
177 neventsmcall= ecutsmc->NumberOfEvents();
181 cout<<neventsdata<<" Events"<<endl;
184 TH1F* allgen=((TH1F*)managermc->GetPtHistogram1D("hHistPtGen",1,1))->Clone();
185 allgen->SetName("AllGen");
186 TH1F* allrecMC=GetOneHistFromPtDCAhisto("hHistPtRec","rawallMC",managermc,dcacutxy);
187 TH1F* alleff=GetOneHistFromPtDCAhisto("hHistPtRecPrimary","effall",managermc,dcacutxy);
188 TH1F* allrecdata=GetOneHistFromPtDCAhisto("hHistPtRec","rawalldata",managerdata,dcacutxy);
193 TH1F* spectraall=(TH1F*)allrecdata->Clone("recNch");
194 spectraall->SetTitle("recNch");
195 TH1F* contall=(TH1F*)allrecMC->Clone("contall");
196 contall->SetTitle("contall");
197 //contall->Add(alleff,-1);
198 SubHistWithFullCorr(contall,alleff);
199 alleff->Divide(alleff,allgen,1,1,"B");
200 contall->Divide(contall,allrecMC,1,1,"B");
202 GetCorrectedSpectra(spectraall,allrecdata,alleff,contall);
203 Divideby2pipt(spectraall);
205 allrecdata->Scale(1./neventsdata,"width");
206 allgen->Scale(1./neventsmcall,"width");
207 allrecMC->Scale(1./neventsmc,"width");
208 spectraall->Scale(1./neventsdata,"width");
214 lout->Add(allrecdata);
215 lout->Add(spectraall);
218 for (int i=0;i<6;i++)
222 TString tmpname(rawspectramc[i]->GetTitle());
223 tmpname.ReplaceAll("SpectraMC","%s");
224 contallMC[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contallMC"));
225 contallMC[i]->SetTitle(Form(tmpname.Data(),"contallMC"));
226 contfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contfit"));
227 contfit[i]->SetTitle(Form(tmpname.Data(),"contfit"));
228 contWDfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contWDfit"));
229 contWDfit[i]->SetTitle(Form(tmpname.Data(),"contWDfit"));
230 contMatfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contMatfit"));
231 contMatfit[i]->SetTitle(Form(tmpname.Data(),"contMatfit"));
232 primaryfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"primaryfit"));
233 primaryfit[i]->SetTitle(Form(tmpname.Data(),"primaryfit"));
234 contfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contfitonMC"));
235 contfit[i+6]->SetTitle(Form(tmpname.Data(),"contfitonMC"));
236 contWDfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contWDfitonMC"));
237 contWDfit[i+6]->SetTitle(Form(tmpname.Data(),"contWDfitonMC"));
238 contMatfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contMatfitonMC"));
239 contMatfit[i+6]->SetTitle(Form(tmpname.Data(),"contMatfitonMC"));
240 primaryfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"primaryfitMC"));
241 primaryfit[i+6]->SetTitle(Form(tmpname.Data(),"primaryfitMC"));
247 contWDfit[i]->Reset();
248 contMatfit[i]->Reset();
249 primaryfit[i]->Reset();
252 contfit[i+6]->Reset();
253 contWDfit[i+6]->Reset();
254 contMatfit[i+6]->Reset();
255 primaryfit[i+6]->Reset();
257 SetBintoOne(primaryfit[i]);
258 SetBintoOne(primaryfit[i+6]);
259 spectra[i]=(TH1F*)rawspectradata[i]->Clone(Form(tmpname.Data(),"SpectraFinal"));
260 spectra[i]->SetTitle(Form(tmpname.Data(),"SpectraFinal"));
261 spectraLeonardo[i]=(TH1F*)rawspectradata[i]->Clone(Form(tmpname.Data(),"SpectraFinalLeonardo"));
262 spectraLeonardo[i]->SetTitle(Form(tmpname.Data(),"SpectraFinalLeonardo"));
264 corrLeonardo[i]=(TH1F*)MCTruth[i]->Clone(Form(tmpname.Data(),"CorrFactLeonardo"));
265 corrLeonardo[i]->SetTitle(Form(tmpname.Data(),"CorrFactLeonardo"));
266 corrLeonardo[i]->Divide(corrLeonardo[i],rawspectramc[i],1,1,"B");
270 //contallMC[i]->Add(eff[i],-1.0);
271 SubHistWithFullCorr(contallMC[i],eff[i]);
272 //RecomputeErrors(contallMC[i]);
273 contallMC[i]->Sumw2();
274 contallMC[i]->Divide(contallMC[i],rawspectramc[i],1,1,"B");
276 // contamintaion from PID but only primaries
277 //contPIDpri[i]->Add(eff[i],-1.0);
278 SubHistWithFullCorr(contPIDpri[i],eff[i]);
280 //RecomputeErrors(contPIDpri[i]);
281 contPIDpri[i]->Divide(contPIDpri[i],rawspectramc[i],1,1,"B");
283 eff[i]->Divide(eff[i],MCTruth[i],1,1,"B");
287 rawspectramc[i]->Sumw2();
288 //contPID[i]->Add(contPID[i],rawspectramc[i],-1,1);
289 SubHistWithFullCorr(contPID[i],rawspectramc[i]);
290 contPID[i]->Scale(-1.0);
292 //RecomputeErrors(contPID[i]);
293 contPID[i]->ResetStats();
295 contPID[i]->Divide(contPID[i],rawspectramc[i],1,1,"B");
297 if(options&kuseprimaryPIDcont)
298 confinal[i]=(TH1F*)contPIDpri[i]->Clone(Form(tmpname.Data(),"confinal"));
300 confinal[i]=(TH1F*)contPID[i]->Clone(Form(tmpname.Data(),"confinal"));
301 confinal[i]->SetTitle(Form(tmpname.Data(),"confinal"));
303 contSecMC[i]=(TH1F*)contWD[i]->Clone(Form(tmpname.Data(),"conSecMC"));
304 contSecMC[i]->Add(contMat[i]);
305 contWD[i]->Divide(contWD[i],rawspectramc[i],1,1,"B");
306 contMat[i]->Divide(contMat[i],rawspectramc[i],1,1,"B");
307 contSecMC[i]->Divide(contSecMC[i],rawspectramc[i],1,1,"B");
310 doubleconuntsdata[i-3]=(TH1F*)rawspectradata[i]->Clone(Form("DoublecountsDATA%s",Particle[i-3].Data()));
311 doubleconuntsdata[i-3]->Reset();
312 doubleconuntsMC[i-3]=(TH1F*)rawspectramc[i]->Clone(Form("DoublecountsMC%s",Particle[i-3].Data()));
313 doubleconuntsMC[i-3]->Reset();
315 CalculateDoubleCounts(doubleconuntsdata[i-3],rawspectradata,i-3,kTRUE);
316 CalculateDoubleCounts(doubleconuntsMC[i-3],rawspectramc,i-3,kFALSE);
324 rawspectradata[i]->Scale(1./neventsdata,"width");
325 rawspectramc[i]->Scale(1./neventsmc,"width");
326 MCTruth[i]->Scale(1./neventsmcall,"width");
327 spectraLeonardo[i]->Scale(1./neventsdata,"width");
331 lout->Add(rawspectradata[i]);
332 lout->Add(rawspectramc[i]);
333 lout->Add(MCTruth[i]);
335 lout->Add(contallMC[i]);
336 lout->Add(contPID[i]);
337 lout->Add(contWD[i]);
338 lout->Add(contMat[i]);
339 lout->Add(contfit[i]);
340 lout->Add(contWDfit[i]);
341 lout->Add(contMatfit[i]);
342 lout->Add(primaryfit[i]);
343 lout->Add(contfit[i+6]);
344 lout->Add(contWDfit[i+6]);
345 lout->Add(contMatfit[i+6]);
346 lout->Add(primaryfit[i+6]);
347 lout->Add(spectra[i]);
348 lout->Add(spectraLeonardo[i]);
349 lout->Add(confinal[i]);
350 lout->Add(contPIDpri[i]);
351 lout->Add(contSecMC[i]);
354 lout->Add(doubleconuntsdata[i-3]);
355 lout->Add(doubleconuntsMC[i-3]);
358 outdate.ReplaceAll("/","_");
359 configfile.ReplaceAll(".","_");
361 if(configfile.Length()>0&&(options&kuserangeonfigfile))
362 fout=new TFile(Form("./results/ResMY_%s_%s_%#X_%s.root",outnamemc.Data(),outdate.Data(),options,configfile.Data()),"RECREATE");
364 fout=new TFile(Form("./results/ResMY_%s_%s_%#X.root",outnamemc.Data(),outdate.Data(),options),"RECREATE");
366 DCACorrectionMarek(managerdata,managermc,dcacutxy,fout,contfit,contWDfit,contMatfit,primaryfit);
367 for (int i=0;i<6;i++)
371 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
372 confinal[i]->Add(contfit[i]);
373 GetCorrectedSpectra(spectra[i],rawspectradata[i],eff[i],confinal[i]);
377 GetCorrectedSpectra(spectra[i],rawspectradata[i],eff[i],contallMC[i]);
379 GetCorrectedSpectraLeonardo(spectraLeonardo[i],corrLeonardo[i],primaryfit[i],primaryfit[i+6]);
380 if(options&kskipconcutonspectra)
382 if(options&kuseprimaryPIDcont)
384 CleanHisto(spectra[i],-1,100,contPIDpri[i]);
385 CleanHisto(spectraLeonardo[i],-1,100,contPIDpri[i]);
389 CleanHisto(spectra[i],-1,100,contPID[i]);
390 CleanHisto(spectraLeonardo[i],-1,100,contPID[i]);
392 // Apply correction for wrongly simulated TOF signal in MC
393 if(options&kuseTOFcorrforPIDsignalmatching)
394 TOFPIDsignalmatchingApply(spectra[i],TOFPIDsignalmatching[i%3]);
398 GFCorrection(spectra,tcutsdata->GetPtTOFMatching(),options);
399 GFCorrection(spectraLeonardo,tcutsdata->GetPtTOFMatching(),options);
401 Double_t ycut=tcutsdata->GetY();
402 Double_t etacut=tcutsdata->GetEtaMax()-tcutsdata->GetEtaMin();
405 cout<<"using eta window 1.6"<<endl;
410 TH1F* allch=GetSumAllCh(spectra,mass,etacut);
412 if(options&kuseTOFmatchingcorrection)
414 MatchingTOFEff(spectra,lout);
415 MatchingTOFEff(spectraLeonardo);
416 TOFMatchingForNch(spectraall);
420 TList* listqa=new TList();
421 TList* canvaslist=new TList();
422 Float_t vertexcorrection=1.0;
423 Float_t corrbadchunksvtx=1.0;
424 if ((options&knormalizationwithbin0integralsdata)||(options&knormalizationwithbin0integralsMC))
425 corrbadchunksvtx=QAPlotsBoth(managerdata,managermc,ecutsdata,ecutsmc,tcutsdata,tcutsmc,listqa,canvaslist,0);
427 corrbadchunksvtx=QAPlotsBoth(managerdata,managermc,ecutsdata,ecutsmc,tcutsdata,tcutsmc,listqa,canvaslist,1);
428 if (options&kveretxcorrectionandbadchunkscorr)
429 vertexcorrection=corrbadchunksvtx;
430 cout<<" VTX corr="<<vertexcorrection<<endl;
431 if(TMath::Abs(ycut)>0.0)
432 vertexcorrection=vertexcorrection/(2.0*ycut);
433 for (int i=0;i<6;i++)
435 spectra[i]->Scale(vertexcorrection);
436 spectraLeonardo[i]->Scale(vertexcorrection);
437 if(TMath::Abs(ycut)>0.0)
439 rawspectradata[i]->Scale(1.0/(2.0*ycut));
440 rawspectramc[i]->Scale(1.0/(2.0*ycut));
441 MCTruth[i]->Scale(1.0/(2.0*ycut));
444 allch->Scale(vertexcorrection);
445 spectraall->Scale(vertexcorrection/etacut);
447 //spectraall->Scale(1.0/1.6);
448 lout->Write("output",TObject::kSingleKey);
449 listqa->Write("outputQA",TObject::kSingleKey);
450 canvaslist->Write("outputcanvas",TObject::kSingleKey);
453 //Normaliztionwithbin0integrals();
457 Bool_t OpenFile(TString dirname,TString outputname, Bool_t mcflag, Bool_t mcasdata)
461 TString nameFile = Form("./%s/AnalysisResults%s.root",dirname.Data(),(mcflag?"MC":"DATA"));
462 TFile *file = TFile::Open(nameFile.Data());
465 cout<<"no file"<<endl;
468 TString sname=Form("OutputBothSpectraTask_%s_%s",(mcflag?"MC":"Data"),outputname.Data());
471 cout<<"using MC as data "<<endl;
472 sname=Form("OutputBothSpectraTask_%s_%s","MC",outputname.Data());
475 TDirectoryFile *dir=(TDirectoryFile*)file->Get(sname.Data());
478 // cout<<"no dir "<<sname.Data()<<endl;
481 cout<<"using MC as data "<<endl;
482 sname=Form("OutputAODSpectraTask_%s_%s","MC",outputname.Data());
485 sname=Form("OutputAODSpectraTask_%s_%s",(mcflag?"MC":"Data"),outputname.Data());
486 // cout<<"trying "<<sname.Data()<<endl;
487 dir=(TDirectoryFile*)file->Get(sname.Data());
490 cout<<"no dir "<<sname.Data()<<endl;
494 cout << " -- Info about " <<(mcflag?"MC":"DATA") <<" -- "<< endl;
497 managermc= (AliSpectraBothHistoManager*) dir->Get("SpectraHistos");
498 ecutsmc = (AliSpectraBothEventCuts*) dir->Get("Event Cuts");
499 tcutsmc = (AliSpectraBothTrackCuts*) dir->Get("Track Cuts");
500 ecutsmc->PrintCuts();
501 tcutsmc->PrintCuts();
502 if(!managermc||!ecutsmc||!tcutsmc)
504 if(managermc->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()!=ecutsmc->GetHistoCuts()->GetBinContent(3))
505 cout<<"Please check MC file possible problem with merging"<<" "<<managermc->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()<<" "<<ecutsmc->GetHistoCuts()->GetBinContent(3)<<endl;
509 managerdata= (AliSpectraBothHistoManager*) dir->Get("SpectraHistos");
510 ecutsdata = (AliSpectraBothEventCuts*) dir->Get("Event Cuts");
511 tcutsdata = (AliSpectraBothTrackCuts*) dir->Get("Track Cuts");
512 ecutsdata->PrintCuts();
513 tcutsdata->PrintCuts();
514 if(!managerdata||!ecutsdata||!tcutsdata)
516 if(managerdata->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()!=ecutsdata->GetHistoCuts()->GetBinContent(3))
517 cout<<"Please check DATA file possible problem with merging"<<" "<<managerdata->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()<<" "<<ecutsdata->GetHistoCuts()->GetBinContent(3)<<endl;
523 void GetMCTruth(TH1F** MCTruth)
525 for(Int_t icharge=0;icharge<2;icharge++)
527 for(Int_t ipart=0;ipart<3;ipart++)
529 Int_t index=ipart+3*icharge;
530 TString hname=Form("hHistPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
531 MCTruth[index]=(TH1F*)((TH1F*)managermc->GetPtHistogram1D(hname.Data(),1,1))->Clone();
532 MCTruth[index]->SetName(Form("MCTruth_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
533 MCTruth[index]->SetTitle(Form("MCTruth_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
534 MCTruth[index]->Sumw2();
539 TH1F* GetOneHistFromPtDCAhisto(TString name,TString hnameout,AliSpectraBothHistoManager* hman,TFormula* dcacutxy)
541 histo =(TH1F*)((TH1F*) hman->GetPtHistogram1D(name.Data(),-1,-1))->Clone();
542 histo->SetName(hnameout.Data());
543 histo->SetTitle(hnameout.Data());
547 for(int ibin=1;ibin<histo->GetNbinsX();ibin++)
549 Double_t lowedge=histo->GetBinLowEdge(ibin);
550 Float_t cut=dcacutxy->Eval(lowedge);
551 TH1F* dcahist=(TH1F*)hman->GetDCAHistogram1D(name.Data(),lowedge,lowedge);
552 // Float_t inyield=dcahist->Integral(dcahist->GetXaxis()->FindBin(-1.0*cut),dcahist->GetXaxis()->FindBin(cut));
553 Float_t testyield=0.0;
554 Float_t testerror=0.0;
555 for (int itest=dcahist->GetXaxis()->FindBin(-1.0*cut);itest<=dcahist->GetXaxis()->FindBin(cut);itest++)
557 testyield+=dcahist->GetBinContent(itest);
558 testerror+=dcahist->GetBinError(itest)*dcahist->GetBinError(itest);
560 // cout<<"corr data "<<histo->GetBinContent(ibin)<<" "<<inyield<<" "<<dcahist->Integral()<<" "<<hnameout.Data()<<endl;
561 // 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;
563 // cout<<testyield<<" "<<TMath::Sqrt(testerror)<<" error2 "<<inyield<<" "<<TMath::Sqrt(inyield)<<endl;
565 //histo->SetBinContent(ibin,inyield);
566 //histo->SetBinError(ibin,TMath::Sqrt(inyield));
567 histo->SetBinContent(ibin,testyield);
568 histo->SetBinError(ibin,TMath::Sqrt(testerror));
579 void GetPtHistFromPtDCAhisto(TString hnamein, TString hnameout, AliSpectraBothHistoManager* hman,TH1F** histo,TFormula* dcacutxy)
581 //Float_t min[3]={0.3,0.3,0.45};
582 //Float_t max[3]={1.5,1.2,2.2};
583 for(Int_t icharge=0;icharge<2;icharge++)
585 for(Int_t ipart=0;ipart<3;ipart++)
587 Int_t index=ipart+3*icharge;
588 Printf("Getting %s",hnamein.Data());
589 TString nameinfinal=Form("%s%s%s",hnamein.Data(),Particle[ipart].Data(),Sign[icharge].Data());
590 TString nameoutfinal=Form("%s%s%s",hnameout.Data(),Particle[ipart].Data(),Sign[icharge].Data());
593 histo[index]=GetOneHistFromPtDCAhisto(nameinfinal,nameoutfinal,hman,dcacutxy);
594 CleanHisto(histo[index],minRanges[ipart],maxRanges[ipart]);
598 void CleanHisto(TH1F* h, Float_t minV, Float_t maxV,TH1* contpid=0x0)
600 for (int i=0;i<=h->GetNbinsX();i++)
602 if(h->GetXaxis()->GetBinCenter(i)<minV||h->GetXaxis()->GetBinCenter(i)>maxV)
604 h->SetBinContent(i,0);
609 if(contpid->GetBinContent(i)>fMaxContaminationPIDMC)
611 h->SetBinContent(i,0);
619 void DCACorrectionMarek(AliSpectraBothHistoManager* hman_data, AliSpectraBothHistoManager* hman_mc,TFormula* fun,TFile *fout,TH1F** hcon,TH1F** hconWD,TH1F** hconMat,TH1F** hprimary)
621 printf("\n\n-> DCA Correction");
625 TString sample[2]={"data","mc"};
626 ofstream debug("debugDCA.txt");
627 TList* listofdcafits=new TList();
628 for(Int_t icharge=0;icharge<2;icharge++)
630 for(Int_t ipart=0;ipart<3;ipart++)
632 Int_t index=ipart+3*icharge;
633 for(Int_t isample=0;isample<2;isample++)
636 for(Int_t ibin_data=7;ibin_data<40;ibin_data++)
639 Double_t lowedge=hcon[index]->GetBinLowEdge(ibin_data);
640 Double_t binwidth=hcon[index]->GetBinWidth(ibin_data);
641 debug<<"NEW "<<Particle[ipart].Data()<<" "<<Sign[icharge].Data()<<" "<<lowedge<<endl;
644 CutRange[1]=fun->Eval(lowedge);
645 CutRange[0]=-1.0*CutRange[1];
647 debug<<"cut "<<CutRange[1]<<" "<<CutRange[0]<<endl;
648 Short_t fitsettings=DCAfitsettings(lowedge+0.5*binwidth,index);
649 debug<<"settings "<< fitsettings<<endl;
653 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);
654 cDCA->SetMargin(0.1,0.02,0.1,0.02);
655 TLegend* Leg1 = new TLegend(0.6,0.6,0.85,0.85,"","NDC");
656 Leg1->SetFillStyle(kFALSE);
657 Leg1->SetLineColor(kWhite);
658 Leg1->SetBorderSize(0);
662 TH1F *hToFit =(TH1F*) ((TH1F*)hman_data->GetDCAHistogram1D(Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
664 TH1F *hToFit =(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
665 debug<<Particle[ipart].Data()<<" "<<Sign[icharge].Data()<<" "<<lowedge<<endl;
666 TH1F *hmc1=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
667 TH1F *hmc2=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaSecondaryWeakDecay%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
668 TH1F *hmc3=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaSecondaryMaterial%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
669 Double_t minentries=1;
670 debug<<hToFit->GetEntries()<<" "<<hmc1->GetEntries()<<" "<<hmc2->GetEntries()<<" "<<hmc3->GetEntries()<<endl;
671 debug<<((fitsettings&0x1)&&hmc2->GetEntries()<=minentries)<<" "<<((fitsettings&0x2)&&hmc3->GetEntries()<=minentries)<<endl;
672 if(hToFit->GetEntries()<=minentries || hmc1->GetEntries()<=minentries || ((fitsettings&0x1)&&hmc2->GetEntries()<=minentries) || ((fitsettings&0x2)&&hmc3->GetEntries()<=minentries))
679 Float_t corrforrebinning[4]={1.0,1.0,1.0,1.0};
680 Int_t binCutRange[]={hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1])};
683 if(hmc3->GetNbinsX()>300)
686 corrforrebinning[0]=hToFit->Integral(binCutRange[0],binCutRange[1]);
687 corrforrebinning[1]=hmc1->Integral(binCutRange[0],binCutRange[1]);
688 corrforrebinning[2]=hmc2->Integral(binCutRange[0],binCutRange[1]);
689 corrforrebinning[3]=hmc3->Integral(binCutRange[0],binCutRange[1]);
696 binCutRange[0]=hmc1->GetXaxis()->FindBin(CutRange[0]);
697 binCutRange[1]=hmc1->GetXaxis()->FindBin(CutRange[1]);
700 //after rebbing we lose resolution of the dca this correction also us to do obtain inside used dca
702 if(hToFit->Integral(binCutRange[0],binCutRange[1])>0.0)
703 corrforrebinning[0]=corrforrebinning[0]/hToFit->Integral(binCutRange[0],binCutRange[1]);
705 corrforrebinning[0]=1.0;
706 if(hmc1->Integral(binCutRange[0],binCutRange[1])>0.0)
707 corrforrebinning[1]=corrforrebinning[1]/hmc1->Integral(binCutRange[0],binCutRange[1]);
709 corrforrebinning[1]=1.0;
710 if(hmc2->Integral(binCutRange[0],binCutRange[1])>0.0)
711 corrforrebinning[2]=corrforrebinning[2]/hmc2->Integral(binCutRange[0],binCutRange[1]);
713 corrforrebinning[2]=1.0;
714 if(hmc3->Integral(binCutRange[0],binCutRange[1])>0.0)
715 corrforrebinning[3]=corrforrebinning[3]/hmc3->Integral(binCutRange[0],binCutRange[1]);
717 corrforrebinning[3]=1.0;
731 mc = new TObjArray(3); // MC histograms are put in this array
734 mc = new TObjArray(2);
743 TFractionFitter* fit = new TFractionFitter(hToFit,mc); // initialise
744 fit->Constrain(0,0.0,1.0); // constrain fraction 1 to be between 0 and 1
745 Double_t defaultStep = 0.01;
746 for (int iparm = 0; iparm < Npar; ++iparm)
748 TString name("frac"); name += iparm;
750 fit->GetFitter()->SetParameter(iparm,name.Data(),0.85,0.01,0.0,1.0);
752 fit->GetFitter()->SetParameter(iparm,name.Data(),0.15,0.01,0.0,1.0);
754 fit->GetFitter()->SetParameter(iparm,name.Data(),0.075,0.01,0.0,1.0);
758 fit->Constrain(1,0.0,1.0); // constrain fraction 1 to be between 0 and 1
760 fit->Constrain(1+(fitsettings&0x1),0.0,1.0); // constrain fraction 1 to be between 0 and 1
762 Int_t binFitRange[]={hmc1->GetXaxis()->FindBin(FitRange[0]),hmc1->GetXaxis()->FindBin(FitRange[1])};
763 fit->SetRangeX(binFitRange[0],binFitRange[1]);
764 hToFit->GetXaxis()->SetRange(binFitRange[0],binFitRange[1]);
765 // hToFit->SetTitle(Form("DCA distr - %s %s %s %lf",sample[isample].Data(),Particle[ipart].Data(),Sign[icharge].Data(),lowedge));
766 hToFit->SetTitle("");
767 Int_t status = fit->Fit(); // perform the fit
768 cout << "fit status: " << status << endl;
769 debug<<"fit status: " << status << endl;
773 Double_t v1=0.0,v2=0.0,v3=0.0;
774 Double_t ev1=0.0,ev2=0.0,ev3=0.0;
777 // check on fit status
778 TH1F* result = (TH1F*) fit->GetPlot();
779 TH1F* PrimMCPred=(TH1F*)fit->GetMCPrediction(0);
781 TH1F* secStMCPred=0X0;
785 //first method, use directly the fit result
786 fit->GetResult(0,v1,ev1);
790 fit->GetResult(1,v2,ev2);
791 secStMCPred=(TH1F*)fit->GetMCPrediction(1);
796 fit->GetResult(1+(fitsettings&0x1),v3,ev3);
797 secMCPred=(TH1F*)fit->GetMCPrediction(1+(fitsettings&0x1));
799 cov=fit->GetFitter()->GetCovarianceMatrixElement(1,2);
805 debug<<v1<<" "<<ev1<<" "<<v2<<" "<<ev2<<" "<<v3<<" "<<ev3<<" "<<" "<<cov<<endl;
807 // becuase dca cut range is not a fit range the results from TFractionFitter should be rescale
808 // This can be done in two ways or use input histograms or output histograms
809 // The difference between those two methods should be on the level of statistical error
810 // I use output histograms
812 // Method 1 input histo
814 Float_t normalizationdata=hToFit->Integral(hToFit->GetXaxis()->FindBin(CutRange[0]),hToFit->GetXaxis()->FindBin(CutRange[1]))/hToFit->Integral(binFitRange[0],binFitRange[1]);
815 normalizationdata*=corrforrebinning[0];
817 Float_t normalizationmc1=(hmc1->Integral(binCutRange[0],binCutRange[1])/hmc1->Integral(binFitRange[0],binFitRange[1]))/normalizationdata;
818 Float_t normalizationmc2=0.0;
820 normalizationmc2=(hmc2->Integral(binCutRange[0],binCutRange[1])/hmc2->Integral(binFitRange[0],binFitRange[1]))/normalizationdata;
821 Float_t normalizationmc3=0.0;
823 normalizationmc3=(hmc3->Integral(binCutRange[0],binCutRange[1])/hmc3->Integral(binFitRange[0],binFitRange[1]))/normalizationdata;
825 normalizationmc1*=corrforrebinning[1];
826 normalizationmc2*=corrforrebinning[2];
827 normalizationmc3*=corrforrebinning[3];
829 debug<<"After Nor"<<endl;
830 debug<<v1*normalizationmc1<<" "<<ev1*normalizationmc1<<" "<<v2*normalizationmc2<<" "<<ev2*normalizationmc2<<" "<<v3*normalizationmc3<<" "<<ev3*normalizationmc3<<" "<<endl;
831 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;
832 debug<<"addtional info"<<endl;
835 //Method 2 output histo
837 Float_t normalizationdata1=result->Integral(binCutRange[0],binCutRange[1])/result->Integral(binFitRange[0],binFitRange[1]);
840 // if the cut range is bigger the fit range we should calculate the normalization factor for data using the data histogram
841 // because result histogram has entries only in fits range
842 if(FitRange[0]>CutRange[0]||FitRange[1]<CutRange[1])
843 normalizationdata1=normalizationdata;
845 normalizationdata1*=corrforrebinning[0];
848 Float_t normalizationmc11=(PrimMCPred->Integral(binCutRange[0],binCutRange[1])/PrimMCPred->Integral(binFitRange[0],binFitRange[1]))/normalizationdata1;
849 Float_t normalizationmc21=0.0;
851 normalizationmc21=(secStMCPred->Integral(binCutRange[0],binCutRange[1])/secStMCPred->Integral(binFitRange[0],binFitRange[1]))/normalizationdata1;
852 Float_t normalizationmc31=0.0;
854 normalizationmc31=(secMCPred->Integral(binCutRange[0],binCutRange[1])/secMCPred->Integral(binFitRange[0],binFitRange[1]))/normalizationdata1;
856 normalizationmc11*=corrforrebinning[1];
857 normalizationmc21*=corrforrebinning[2];
858 normalizationmc31*=corrforrebinning[3];
860 debug<<"After Nor 2"<<endl;
861 debug<<v1*normalizationmc11<<" "<<ev1*normalizationmc11<<" "<<v2*normalizationmc21<<" "<<ev2*normalizationmc21<<" "<<v3*normalizationmc31<<" "<<ev3*normalizationmc31<<endl;
863 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;
866 hconWD[index+6*isample]->SetBinContent(ibin_data,v2*normalizationmc21);
867 hconWD[index+6*isample]->SetBinError(ibin_data,ev2*normalizationmc21);
868 hconMat[index+6*isample]->SetBinContent(ibin_data,v3*normalizationmc31);
869 hconMat[index+6*isample]->SetBinError(ibin_data,ev3*normalizationmc31);
870 hprimary[index+6*isample]->SetBinContent(ibin_data,v1*normalizationmc11);
871 hprimary[index+6*isample]->SetBinError(ibin_data,ev1*normalizationmc11);
872 hcon[index+6*isample]->SetBinContent(ibin_data,v2*normalizationmc21+v3*normalizationmc31);
873 hcon[index+6*isample]->SetBinError(ibin_data,TMath::Sqrt(ev2*ev2*normalizationmc21*normalizationmc21+ev3*ev3*normalizationmc31*normalizationmc31+cov*normalizationmc31*normalizationmc21));
878 result->Scale(1.0/result->Integral(result->GetXaxis()->FindBin(FitRange[0]),result->GetXaxis()->FindBin(FitRange[1])));
879 hToFit->Scale(1.0/hToFit->Integral(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1])));
880 PrimMCPred->Scale(v1/PrimMCPred->Integral(PrimMCPred->GetXaxis()->FindBin(FitRange[0]),PrimMCPred->GetXaxis()->FindBin(FitRange[1])));
882 hToFit->SetMinimum(0.0001);
883 hToFit->GetXaxis()->SetTitle("DCA_{xy} (cm)");
884 hToFit->GetYaxis()->SetTitle("N_{counts}/N_{counts}(-3cm;3cm)");
885 hToFit->GetYaxis()->SetTitleOffset(1.3);
886 hToFit->DrawClone("E1x0");
887 Leg1->AddEntry(hToFit,"data","p");
888 result->SetTitle("Fit result");
889 result->SetLineColor(kBlack);
890 Leg1->AddEntry(result,"fit result","l");
891 result->DrawClone("histsame");
893 PrimMCPred->SetLineColor(kGreen+2);
894 PrimMCPred->SetLineStyle(2);
895 PrimMCPred->SetLineWidth(3.0);
896 Leg1->AddEntry(PrimMCPred,"primaries","l");
897 PrimMCPred->DrawClone("histsame");
901 secStMCPred->Scale(v2/secStMCPred->Integral(secStMCPred->GetXaxis()->FindBin(FitRange[0]),secStMCPred->GetXaxis()->FindBin(FitRange[1])));
902 secStMCPred->SetLineColor(kRed);
903 secStMCPred->SetLineWidth(3.0);
905 secStMCPred->SetLineStyle(3);
906 Leg1->AddEntry(secStMCPred,"weak decays","l");
907 secStMCPred->DrawClone("histsame");
913 secMCPred->Scale(v3/secMCPred->Integral(secMCPred->GetXaxis()->FindBin(FitRange[0]),secMCPred->GetXaxis()->FindBin(FitRange[1])));
914 secMCPred->SetLineColor(kBlue);
915 secMCPred->SetLineWidth(3.0);
917 secMCPred->SetLineStyle(4);
918 Leg1->AddEntry(secMCPred,"material","l");
919 secMCPred->DrawClone("histsame");
925 hconWD[index+6*isample]->SetBinContent(ibin_data,0.0);
926 hconWD[index+6*isample]->SetBinError(ibin_data,0.0);
927 hconMat[index+6*isample]->SetBinContent(ibin_data,0.0);
928 hconMat[index+6*isample]->SetBinError(ibin_data,0.0);
929 hcon[index+6*isample]->SetBinContent(ibin_data,0.0);
930 hcon[index+6*isample]->SetBinError(ibin_data,0.0);
931 hprimary[index+6*isample]->SetBinContent(ibin_data,1.0);
932 hprimary[index+6*isample]->SetBinError(ibin_data,0.0);
935 TLatex* texttitle=new TLatex();
937 texttitle->SetTextSize(0.04);
938 texttitle->DrawLatex(0.12,0.92,Form("%s %.2f<#it{p}_{T} < %.2f (GeV/#it{c})",symboles[index].Data(),lowedge,binwidth+lowedge));
939 listofdcafits->Add(cDCA);
951 listofdcafits->Write("DCAfits",TObject::kSingleKey);
954 void RecomputeErrors(TH1* h)
956 for (int i=0; i<=h->GetXaxis()->GetNbins(); i++)
958 cout<<h->GetBinContent(i)<<" "<<h->GetBinError(i)<<" error "<<TMath::Sqrt(h->GetBinContent(i))<<endl;
959 h->SetBinError(i,TMath::Sqrt(h->GetBinContent(i)));
963 void SetBintoOne(TH1* h)
965 for (int i=0;i<=h->GetXaxis()->GetNbins(); i++)
967 h->SetBinContent(i,1);
973 void GetCorrectedSpectra(TH1F* corr,TH1F* raw,TH1F* eff, TH1F* con)
975 for (int i=0;i<=corr->GetXaxis()->GetNbins(); i++)
977 corr->SetBinContent(i,1);
978 corr->SetBinError(i,0);
988 void GetCorrectedSpectraLeonardo(TH1F* spectra,TH1F* correction, TH1F* hprimaryData,TH1F* hprimaryMC)
991 spectra->Multiply(correction);
993 hprimaryData->Sumw2();
994 spectra->Multiply(hprimaryData);
996 spectra->Divide(hprimaryMC);
999 void GFCorrection(TH1F **Spectra,Float_t tofpt,UInt_t options)
1001 if (options&kgeantflukaKaon)
1003 //Geant/Fluka Correction
1004 Printf("\nGF correction for Kaons");
1005 //Getting GF For Kaons in TPC
1006 TGraph *gGFCorrectionKaonPlus=new TGraph();
1007 gGFCorrectionKaonPlus->SetName("gGFCorrectionKaonPlus");
1008 gGFCorrectionKaonPlus->SetTitle("gGFCorrectionKaonPlus");
1009 TGraph *gGFCorrectionKaonMinus=new TGraph();
1010 gGFCorrectionKaonMinus->SetName("gGFCorrectionKaonMinus");
1011 gGFCorrectionKaonMinus->SetTitle("gGFCorrectionKaonMinus");
1012 TString fnameGeanFlukaK="GFCorrection/correctionForCrossSection.321.root";
1013 TFile *fGeanFlukaK= new TFile(fnameGeanFlukaK.Data());
1016 fnameGeanFlukaK="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/correctionForCrossSection.321.root";
1017 fGeanFlukaK= new TFile(fnameGeanFlukaK.Data());
1021 TH1F *hGeantFlukaKPos=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles");
1022 TH1F *hGeantFlukaKNeg=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles");
1023 //getting GF func for Kaons with TOF
1024 TF1 *fGFKPosTracking;
1025 fGFKPosTracking = TrackingEff_geantflukaCorrection(3,kPositive);
1026 TF1 *fGFKNegTracking;
1027 fGFKNegTracking = TrackingEff_geantflukaCorrection(3,kNegative);
1028 TF1 *fGFKPosMatching;
1029 fGFKPosMatching = TOFmatchMC_geantflukaCorrection(3,kPositive);
1030 TF1 *fGFKNegMatching;
1031 fGFKNegMatching = TOFmatchMC_geantflukaCorrection(3,kNegative);
1032 for(Int_t binK=0;binK<=Spectra[1]->GetNbinsX();binK++)
1034 if(Spectra[1]->GetBinCenter(binK)<tofpt)
1035 {//use TPC GeantFlukaCorrection
1036 Float_t FlukaCorrKPos=hGeantFlukaKPos->GetBinContent(hGeantFlukaKPos->FindBin(Spectra[1]->GetBinCenter(binK)));
1037 Float_t FlukaCorrKNeg=hGeantFlukaKNeg->GetBinContent(hGeantFlukaKNeg->FindBin(Spectra[4]->GetBinCenter(binK)));
1038 // Printf("TPC Geant/Fluka: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPos,FlukaCorrKNeg);
1039 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPos);
1040 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNeg);
1041 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPos);
1042 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNeg);
1043 gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),FlukaCorrKPos);
1044 gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),FlukaCorrKNeg);
1048 gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),0);
1049 gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),0);
1050 Float_t FlukaCorrKPosTracking=fGFKPosTracking->Eval(Spectra[1]->GetBinCenter(binK));
1051 Float_t FlukaCorrKNegTracking=fGFKNegTracking->Eval(Spectra[1]->GetBinCenter(binK));
1052 // Printf("TPC/TOF Geant/Fluka Tracking: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosTracking,FlukaCorrKNegTracking);
1053 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosTracking);
1054 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegTracking);
1055 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosTracking);
1056 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegTracking);
1057 Float_t FlukaCorrKPosMatching=fGFKPosMatching->Eval(Spectra[1]->GetBinCenter(binK));
1058 Float_t FlukaCorrKNegMatching=fGFKNegMatching->Eval(Spectra[1]->GetBinCenter(binK));
1059 // Printf("TPC/TOF Geant/Fluka Matching: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosMatching,FlukaCorrKNegMatching);
1060 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosMatching);
1061 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegMatching);
1062 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosMatching);
1063 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegMatching);
1067 if(!(options&kgeantflukaProton))
1069 //Geant Fluka for P in TPC
1070 Printf("\nGF correction for Protons");
1071 const Int_t kNCharge=2;
1074 TString fnameGFProtons= "GFCorrection/correctionForCrossSection.root";
1075 TFile* fGFProtons = new TFile (fnameGFProtons.Data());
1078 fnameGFProtons=="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/correctionForCrossSection.root";
1079 TFile* fGFProtons = new TFile (fnameGFProtons.Data());
1086 TH2D * hCorrFluka[kNCharge];
1087 TH2D * hCorrFluka[2];
1088 hCorrFluka[kPos] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionProtons");
1089 hCorrFluka[kNeg] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionAntiProtons");
1090 //getting GF func for Kaons with TPCTOF
1091 TF1 *fGFpPosTracking;
1092 fGFpPosTracking = TrackingEff_geantflukaCorrection(4,kPositive);
1093 TF1 *fGFpNegTracking;
1094 fGFpNegTracking = TrackingEff_geantflukaCorrection(4,kNegative);
1095 TF1 *fGFpPosMatching;
1096 fGFpPosMatching = TOFmatchMC_geantflukaCorrection(4,kPositive);
1097 TF1 *fGFpNegMatching;
1098 fGFpNegMatching = TOFmatchMC_geantflukaCorrection(4,kNegative);
1101 Int_t nbins = Spectra[2]->GetNbinsX();
1103 for(Int_t ibin = 0; ibin < nbins; ibin++)
1105 if(Spectra[2]->GetBinCenter(ibin)<tofpt)
1106 {//use TPC GeantFlukaCorrection
1107 for(Int_t icharge = 0; icharge < kNCharge; icharge++)
1109 Int_t nbinsy=hCorrFluka[icharge]->GetNbinsY();
1110 Float_t pt = Spectra[2]->GetBinCenter(ibin);
1111 Float_t minPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(1);
1112 Float_t maxPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1);
1113 if (pt < minPtCorrection)
1114 pt = minPtCorrection+0.0001;
1115 if (pt > maxPtCorrection)
1116 pt = maxPtCorrection;
1117 Float_t correction = hCorrFluka[icharge]->GetBinContent(1,hCorrFluka[icharge]->GetYaxis()->FindBin(pt));
1119 if (correction > 0.0)
1120 {// If the bin is empty this is a 0
1121 // cout<<icharge<<" "<<ibin<<" "<<correction<<endl;
1122 Spectra[icharge*3+2]->SetBinContent(ibin,Spectra[icharge*3+2]->GetBinContent(ibin)*correction);
1123 Spectra[icharge*3+2]->SetBinError(ibin,Spectra[icharge*3+2]->GetBinError(ibin)*correction);
1125 else if (Spectra[icharge*3+2]->GetBinContent(ibin) > 0.0)
1127 // If we are skipping a non-empty bin, we notify the user
1128 cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for " <<"protons Pt:"<< pt<< endl;
1129 cout << " Bin content: " << Spectra[icharge*3+2]->GetBinContent(ibin) << endl;
1135 Float_t FlukaCorrpPosTracking=fGFpPosTracking->Eval(Spectra[2]->GetBinCenter(ibin));
1136 Float_t FlukaCorrpNegTracking=fGFpNegTracking->Eval(Spectra[2]->GetBinCenter(ibin));
1137 // Printf("TPC/TOF Geant/Fluka Tracking: pt:%f Pos:%f Neg:%f",Spectra[2]->GetBinCenter(ibin),FlukaCorrpPosTracking,FlukaCorrpNegTracking);
1138 Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*FlukaCorrpPosTracking);
1139 Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*FlukaCorrpNegTracking);
1140 Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError(ibin)*FlukaCorrpPosTracking);
1141 Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError(ibin)*FlukaCorrpNegTracking);
1142 Float_t FlukaCorrpPosMatching=fGFpPosMatching->Eval(Spectra[2]->GetBinCenter(ibin));
1143 Float_t FlukaCorrpNegMatching=fGFpNegMatching->Eval(Spectra[2]->GetBinCenter(ibin));
1144 // Printf("TPC/TOF Geant/Fluka Matching: pt:%f Pos:%f Neg:%f",Spectra[2]->GetBinCenter(ibin),FlukaCorrpPosMatching,FlukaCorrpNegMatching);
1145 Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*FlukaCorrpPosMatching);
1146 Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*FlukaCorrpNegMatching);
1147 Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError(ibin)*FlukaCorrpPosMatching);
1148 Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError(ibin)*FlukaCorrpNegMatching);
1151 fGeanFlukaK->Close();
1158 TrackingEff_geantflukaCorrection(Int_t ipart, Int_t icharge)
1161 if (ipart == 3 && icharge == kNegative) {
1162 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.);
1165 else if (ipart == 4 && icharge == kNegative) {
1166 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.);
1169 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionNull(x)", 0., 5.);
1175 TrackingPtGeantFlukaCorrectionNull(Double_t pTmc)
1181 TrackingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
1183 return (1 - 0.129758 *TMath::Exp(-pTmc*0.679612));
1187 TrackingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
1189 return TMath::Min((0.972865 + 0.0117093*pTmc), 1.);
1191 ///////////////////////////////////////////
1193 TOFmatchMC_geantflukaCorrection(Int_t ipart, Int_t icharge)
1196 if (ipart == 3 && icharge == kNegative) {
1197 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.);
1200 else if (ipart == 4 && icharge == kNegative) {
1201 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.);
1204 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionNull(x)", 0., 5.);
1211 MatchingPtGeantFlukaCorrectionNull(Double_t pTmc)
1217 MatchingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
1219 Float_t ptTPCoutP =pTmc*(1-6.81059e-01*TMath::Exp(-pTmc*4.20094));
1220 return (TMath::Power(1 - 0.129758*TMath::Exp(-ptTPCoutP*0.679612),0.07162/0.03471));
1224 MatchingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
1226 Float_t ptTPCoutK=pTmc*(1- 3.37297e-03/pTmc/pTmc - 3.26544e-03/pTmc);
1227 return TMath::Min((TMath::Power(0.972865 + 0.0117093*ptTPCoutK,0.07162/0.03471)), 1.);
1229 void MatchingTOFEff(TH1F** Spectra, TList* list=0x0)
1231 if(TOFMatchingScalling[0]<0.0&&TOFMatchingScalling[1]<0.0)
1233 TH1F *hMatcEffPos_data=(TH1F*)tcutsdata->GetHistoNMatchedPos();
1234 hMatcEffPos_data->Sumw2();
1235 //hMatcEffPos_data->Divide((TH1F*)tcutsdata->GetHistoNSelectedPos());
1236 hMatcEffPos_data->Divide(hMatcEffPos_data,(TH1F*)tcutsdata->GetHistoNSelectedPos(),1,1,"B");
1237 hMatcEffPos_data->SetTitle("Matching Eff Pos - data");
1238 TH1F *hMatcEffNeg_data=(TH1F*)tcutsdata->GetHistoNMatchedNeg();
1239 hMatcEffNeg_data->Sumw2();
1240 //hMatcEffNeg_data->Divide((TH1F*)tcutsdata->GetHistoNSelectedNeg());
1241 hMatcEffNeg_data->Divide(hMatcEffNeg_data,(TH1F*)tcutsdata->GetHistoNSelectedNeg(),1,1,"B");
1242 hMatcEffNeg_data->SetTitle("Matching Eff Neg - data");
1243 TH1F *hMatcEffPos_mc=(TH1F*)tcutsmc->GetHistoNMatchedPos();
1244 hMatcEffPos_mc->Sumw2();
1245 //hMatcEffPos_mc->Divide((TH1F*)tcutsmc->GetHistoNSelectedPos());
1246 hMatcEffPos_mc->Divide(hMatcEffPos_mc,(TH1F*)tcutsmc->GetHistoNSelectedPos(),1,1,"B");
1247 hMatcEffPos_mc->SetTitle("Matching Eff Pos - mc");
1248 TH1F *hMatcEffNeg_mc=(TH1F*)tcutsmc->GetHistoNMatchedNeg();
1249 hMatcEffNeg_mc->Sumw2();
1250 //hMatcEffNeg_mc->Divide((TH1F*)tcutsmc->GetHistoNSelectedNeg());
1251 hMatcEffNeg_mc->Divide(hMatcEffNeg_mc,(TH1F*)tcutsmc->GetHistoNSelectedNeg(),1,1,"B");
1252 hMatcEffNeg_mc->SetTitle("Matching Eff Neg - mc");
1255 hMatcEffPos_data->Divide(hMatcEffPos_mc);
1256 hMatcEffNeg_data->Divide(hMatcEffNeg_mc);
1257 hMatcEffPos_data->SetName("MatchingTOFPos");
1258 hMatcEffNeg_data->SetName("MatchingTOFNeg");
1261 TF1 *pol0MatchPos_data=new TF1("pol0MatchPos_data","pol0",.6,5);
1262 hMatcEffPos_data->Fit("pol0MatchPos_data","MNR");
1263 TF1 *pol0MatchNeg_data=new TF1("pol0MatchNeg_data","pol0",.6,5);
1264 hMatcEffNeg_data->Fit("pol0MatchNeg_data","MNR");
1266 list->Add(hMatcEffPos_data);
1267 list->Add(hMatcEffNeg_data);
1270 TOFMatchingScalling[0]=pol0MatchPos_data->GetParameter(0);
1271 TOFMatchingScalling[1]=pol0MatchNeg_data->GetParameter(0);
1274 //Correction spectra for matching efficiency
1275 //For the moment I'm using the inclusive correction
1276 for(Int_t ipart=0;ipart<3;ipart++)
1279 for(Int_t ibin=1;ibin<Spectra[ipart]->GetNbinsX();ibin++)
1281 Float_t ptspectra=Spectra[ipart]->GetBinCenter(ibin);
1282 if(ptspectra<tcutsdata->GetPtTOFMatching())
1284 //Spectra[ipart]->SetBinContent(ibin,( Spectra[ipart]->GetBinContent(ibin)/hMatcEffPos_data->GetBinContent(hMatcEffPos_data->FindBin(ptspectra))));
1285 //Spectra[ipart+3]->SetBinContent(ibin,( Spectra[ipart+3]->GetBinContent(ibin)/hMatcEffNeg_data->GetBinContent(hMatcEffNeg_data->FindBin(ptspectra))));
1286 Spectra[ipart]->SetBinContent(ibin,( Spectra[ipart]->GetBinContent(ibin)/TOFMatchingScalling[0]));
1287 Spectra[ipart+3]->SetBinContent(ibin,( Spectra[ipart+3]->GetBinContent(ibin)/TOFMatchingScalling[1]));
1291 Double_t eta2y(Double_t pt, Double_t mass, Double_t eta)
1293 Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
1294 return TMath::ASinH(pt / mt * TMath::SinH(eta));
1297 TH1* GetSumAllCh(TH1F** spectra, Double_t* mass,Double_t etacut)
1299 TH1F* allch=(((TH1F*))spectra[0]->Clone("allCh"));
1301 for (int i=0;i<6;i++)
1303 Double_t masstmp=mass[i%3];
1304 for (int j=1;j<spectra[i]->GetXaxis()->GetNbins();j++)
1306 Float_t value=spectra[i]->GetBinContent(j);
1307 Float_t error=spectra[i]->GetBinError(j);
1310 Float_t pt=spectra[i]->GetXaxis()->GetBinCenter(j);
1311 Float_t mt2=masstmp*masstmp+pt*pt;
1312 Float_t mt=TMath::Sqrt(mt2);
1313 Float_t maxy=eta2y(pt,masstmp,etacut/2.0);
1314 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))));
1315 conver=conver/etacut;
1316 //cout<<maxy<<" "<<conver<<" "<<masstmp<<""<<spectra[i]->GetName()<<endl;
1317 Float_t bincontent=allch->GetBinContent(j);
1318 Float_t binerror=allch->GetBinError(j);
1319 bincontent+=conver*value;
1320 binerror=TMath::Sqrt(binerror*binerror+conver*conver*error*error);
1321 allch->SetBinContent(j,bincontent);
1322 allch->SetBinError(j,binerror);
1327 Divideby2pipt(allch);
1328 allch->SetTitle("N_{ch};p_{T} (GeV/c);1/(2 #pi p_{T})dN/p_{T} |#eta|<0.8");
1332 void Divideby2pipt(TH1* hist)
1335 for (int i=1;i<hist->GetXaxis()->GetNbins();i++)
1337 Float_t value=hist->GetBinContent(i);
1338 Float_t error=hist->GetBinError(i);
1339 Float_t pt=hist->GetXaxis()->GetBinCenter(i);
1340 hist->SetBinContent(i,value/(2.0*TMath::Pi()*pt));
1341 hist->SetBinError(i,error/(2.0*TMath::Pi()*pt));
1346 Short_t DCAfitsettings (Float_t pt, Int_t type)
1349 if (pt<maxptformaterial[type]&&pt>minptformaterial[type])
1351 if (pt<maxptforWD[type]&&pt>minptforWD[type])
1357 Float_t Normaliztionwithbin0integrals(UInt_t options)
1360 TH1F* bin0mcRec=(TH1F*)ecutsmc->GetHistoVtxGenerated()->Clone("Bin0_rec");
1361 TH1F* bin0mcMC=(TH1F*)ecutsmc->GetHistoVtxGenerated()->Clone("Bin0_MC");
1363 TH1F* vertexmc=ecutsmc->GetHistoVtxAftSelwithoutZvertexCut();
1364 TH1F* vertexmcMCz=ecutsmc->GetHistoVtxAftSelwithoutZvertexCutusingMCz();
1365 TH1F* vertexdata=ecutsdata->GetHistoVtxAftSelwithoutZvertexCut();
1367 TH1I* histodata=ecutsdata->GetHistoCuts();
1368 TH1I* histomc=ecutsmc->GetHistoCuts();
1370 Float_t dataevents=(Float_t)histodata->GetBinContent(3);
1371 //cout<<histodata->GetBinContent(2)<<endl;
1372 Float_t databin0events=((Float_t)histodata->GetBinContent(2))-((Float_t)histodata->GetBinContent(4));
1377 bin0mcRec->Add(vertexmc,-1);
1378 bin0mcMC->Add(vertexmcMCz,-1);
1380 bin0mcRec->Divide(vertexmc);
1381 bin0mcMC->Divide(vertexmcMCz);
1383 bin0mcRec->Multiply(vertexdata);
1384 bin0mcMC->Multiply(vertexdata);
1386 Float_t bin0mcRecN=0.0;
1387 Float_t bin0mcMCN=0.0;
1389 for (int i=0;i<=bin0mcRec->GetXaxis()->GetNbins();i++)
1391 bin0mcRecN+=bin0mcRec->GetBinContent(i);
1392 bin0mcMCN+=bin0mcMC->GetBinContent(i);
1395 bin0mcRec->Scale(databin0events/bin0mcRecN);
1396 bin0mcMC->Scale(databin0events/bin0mcMCN);
1398 Int_t binmin=bin0mcRec->GetXaxis()->FindBin(-10);
1399 Int_t binmax=bin0mcRec->GetXaxis()->FindBin(10)-1;
1400 cout<< bin0mcRec->GetXaxis()->GetBinLowEdge(binmin)<<" "<<bin0mcRec->GetXaxis()->GetBinUpEdge(binmax)<<endl;
1401 cout<<bin0mcRecN<<" "<<bin0mcMCN<<" "<<databin0events<<endl;
1402 cout<<dataevents<<" normalization "<<dataevents+bin0mcRec->Integral(binmin,binmax)<<" "<<dataevents+bin0mcMC->Integral(binmin,binmax)<<endl;
1403 cout<<histodata->GetBinContent(2)<<" "<<histodata->GetBinContent(4)<<endl;
1404 if ((options&knormalizationwithbin0integralsdata)==knormalizationwithbin0integralsdata)
1405 return dataevents+bin0mcRec->Integral(binmin,binmax);
1406 else if ((options&kknormalizationwithbin0integralsMC)==knormalizationwithbin0integralsMC)
1407 return dataevents+bin0mcMC->Integral(binmin,binmax) ;
1413 Bool_t ReadConfigFile(TString configfile)
1415 ifstream infile(configfile.Data());
1416 if(infile.is_open()==false)
1418 TString namesofSetting[38]={"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"};
1421 while (infile.eof()==false)
1424 while (buffer[0]=='#'&&infile.eof()==false)
1425 infile.getline(buffer,256);
1426 TString tmpstring(buffer);
1428 if(tmpstring.Contains(namesofSetting[0]))
1429 CutRange[0]=(tmpstring.Remove(0,namesofSetting[0].Length()+1)).Atof();
1430 else if (tmpstring.Contains(namesofSetting[1]))
1431 CutRange[1]=(tmpstring.Remove(0,namesofSetting[1].Length()+1)).Atof();
1432 else if (tmpstring.Contains(namesofSetting[2]))
1433 FitRange[0]=(tmpstring.Remove(0,namesofSetting[2].Length()+1)).Atof();
1434 else if (tmpstring.Contains(namesofSetting[3]))
1435 FitRange[1]=(tmpstring.Remove(0,namesofSetting[3].Length()+1)).Atof();
1436 else if (tmpstring.Contains(namesofSetting[4]))
1437 minptformaterial[0]=(tmpstring.Remove(0,namesofSetting[4].Length()+1)).Atof();
1438 else if (tmpstring.Contains(namesofSetting[5]))
1439 maxptformaterial[0]=(tmpstring.Remove(0,namesofSetting[5].Length()+1)).Atof();
1440 else if (tmpstring.Contains(namesofSetting[6]))
1441 minptformaterial[1]=(tmpstring.Remove(0,namesofSetting[6].Length()+1)).Atof();
1442 else if (tmpstring.Contains(namesofSetting[7]))
1443 maxptformaterial[1]=(tmpstring.Remove(0,namesofSetting[7].Length()+1)).Atof();
1444 else if (tmpstring.Contains(namesofSetting[8]))
1445 minptformaterial[2]=(tmpstring.Remove(0,namesofSetting[8].Length()+1)).Atof();
1446 else if (tmpstring.Contains(namesofSetting[9]))
1447 maxptformaterial[2]=(tmpstring.Remove(0,namesofSetting[9].Length()+1)).Atof();
1448 else if (tmpstring.Contains(namesofSetting[10]))
1449 minptformaterial[3]=(tmpstring.Remove(0,namesofSetting[10].Length()+1)).Atof();
1450 else if (tmpstring.Contains(namesofSetting[11]))
1451 maxptformaterial[3]=(tmpstring.Remove(0,namesofSetting[11].Length()+1)).Atof();
1452 else if (tmpstring.Contains(namesofSetting[12]))
1453 minptformaterial[4]=(tmpstring.Remove(0,namesofSetting[12].Length()+1)).Atof();
1454 else if (tmpstring.Contains(namesofSetting[13]))
1455 maxptformaterial[4]=(tmpstring.Remove(0,namesofSetting[13].Length()+1)).Atof();
1456 else if (tmpstring.Contains(namesofSetting[14]))
1457 minptformaterial[5]=(tmpstring.Remove(0,namesofSetting[14].Length()+1)).Atof();
1458 else if (tmpstring.Contains(namesofSetting[15]))
1459 maxptformaterial[5]=(tmpstring.Remove(0,namesofSetting[15].Length()+1)).Atof();
1460 else if (tmpstring.Contains(namesofSetting[16]))
1461 minptforWD[0]=(tmpstring.Remove(0,namesofSetting[16].Length()+1)).Atof();
1462 else if (tmpstring.Contains(namesofSetting[17]))
1463 maxptforWD[0]=(tmpstring.Remove(0,namesofSetting[17].Length()+1)).Atof();
1464 else if (tmpstring.Contains(namesofSetting[18]))
1465 minptforWD[1]=(tmpstring.Remove(0,namesofSetting[18].Length()+1)).Atof();
1466 else if (tmpstring.Contains(namesofSetting[19]))
1467 maxptforWD[1]=(tmpstring.Remove(0,namesofSetting[19].Length()+1)).Atof();
1468 else if (tmpstring.Contains(namesofSetting[20]))
1469 minptforWD[2]=(tmpstring.Remove(0,namesofSetting[20].Length()+1)).Atof();
1470 else if (tmpstring.Contains(namesofSetting[21]))
1471 maxptforWD[2]=(tmpstring.Remove(0,namesofSetting[21].Length()+1)).Atof();
1472 else if (tmpstring.Contains(namesofSetting[22]))
1473 minptforWD[3]=(tmpstring.Remove(0,namesofSetting[22].Length()+1)).Atof();
1474 else if (tmpstring.Contains(namesofSetting[23]))
1475 maxptforWD[3]=(tmpstring.Remove(0,namesofSetting[23].Length()+1)).Atof();
1476 else if (tmpstring.Contains(namesofSetting[24]))
1477 minptforWD[4]=(tmpstring.Remove(0,namesofSetting[24].Length()+1)).Atof();
1478 else if (tmpstring.Contains(namesofSetting[25]))
1479 maxptforWD[4]=(tmpstring.Remove(0,namesofSetting[25].Length()+1)).Atof();
1480 else if (tmpstring.Contains(namesofSetting[26]))
1481 minptforWD[5]=(tmpstring.Remove(0,namesofSetting[26].Length()+1)).Atof();
1482 else if (tmpstring.Contains(namesofSetting[27]))
1483 maxptforWD[5]=(tmpstring.Remove(0,namesofSetting[27].Length()+1)).Atof();
1484 else if (tmpstring.Contains(namesofSetting[28]))
1485 fMaxContaminationPIDMC=(tmpstring.Remove(0,namesofSetting[28].Length()+1)).Atof();
1486 else if (tmpstring.Contains(namesofSetting[29]))
1487 minRanges[0]=(tmpstring.Remove(0,namesofSetting[29].Length()+1)).Atof();
1488 else if (tmpstring.Contains(namesofSetting[30]))
1489 maxRanges[0]=(tmpstring.Remove(0,namesofSetting[30].Length()+1)).Atof();
1490 else if (tmpstring.Contains(namesofSetting[31]))
1491 minRanges[1]=(tmpstring.Remove(0,namesofSetting[31].Length()+1)).Atof();
1492 else if (tmpstring.Contains(namesofSetting[32]))
1493 maxRanges[1]=(tmpstring.Remove(0,namesofSetting[32].Length()+1)).Atof();
1494 else if (tmpstring.Contains(namesofSetting[33]))
1495 minRanges[2]=(tmpstring.Remove(0,namesofSetting[33].Length()+1)).Atof();
1496 else if (tmpstring.Contains(namesofSetting[34]))
1497 maxRanges[2]=(tmpstring.Remove(0,namesofSetting[34].Length()+1)).Atof();
1498 else if (tmpstring.Contains(namesofSetting[35]))
1499 TOFPIDsignalmatching[0]=(tmpstring.Remove(0,namesofSetting[35].Length()+1)).Atof();
1500 else if (tmpstring.Contains(namesofSetting[36]))
1501 TOFPIDsignalmatching[1]=(tmpstring.Remove(0,namesofSetting[36].Length()+1)).Atof();
1502 else if (tmpstring.Contains(namesofSetting[37]))
1503 TOFPIDsignalmatching[2]=(tmpstring.Remove(0,namesofSetting[37].Length()+1)).Atof();
1508 // Double_t minRanges[3]={0.3,0.3,0.45};
1509 //Double_t maxRanges[3]={1.5,1.2,2.2};
1510 //Double_t fMaxContaminationPIDMC=0.2;
1515 for(int i=0;i<6;i++)
1516 cout<<minptformaterial[i]<<" "<<maxptformaterial[i]<<" "<<minptforWD[i]<<" "<<maxptforWD[i]<<endl;
1517 cout<<FitRange[0]<<" "<<FitRange[1]<<" "<<CutRange[0]<<CutRange[1]<<endl;
1518 if(FitRange[0]>=FitRange[1])
1523 if(CutRange[0]>=CutRange[1])
1528 for(int i=0;i<6;i++)
1530 if((minptformaterial[i]>maxptformaterial[i]&&minptformaterial[i]>0.0)||minptformaterial[i]<0.0||maxptformaterial[i]<0.0)
1535 if((minptforWD[i]>maxptforWD[i]&&minptforWD[i]>0.0)||minptforWD[i]<0.0||maxptforWD[i]<0.0)
1541 for(int i=0;i<3;i++)
1542 if(minRanges[i]>maxRanges[i])
1549 void SubHistWithFullCorr(TH1F* h1, TH1F* h2, Float_t factor1=1.0, Float_t factor2=1.0)
1551 if(h1->GetNbinsX()!=h2->GetNbinsX())
1553 for (int i=0;i<=h1->GetNbinsX();i++)
1555 Float_t tmpvalue=factor1*h1->GetBinContent(i)-factor2*h2->GetBinContent(i);
1556 Float_t tmperror=TMath::Abs(factor1*factor1*h1->GetBinError(i)*h1->GetBinError(i)-factor2*factor2*h2->GetBinError(i)*h2->GetBinError(i));
1557 h1->SetBinContent(i,tmpvalue);
1558 h1->SetBinError(i,TMath::Sqrt(tmperror));
1563 void TOFMatchingForNch(TH1* h)
1565 if(TOFMatchingScalling[0]>0.0&&TOFMatchingScalling[1]>0.0)
1567 Float_t factor=0.5*TOFMatchingScalling[0]+0.5*TOFMatchingScalling[1];
1568 for(Int_t ibin=1;ibin<h->GetNbinsX();ibin++)
1570 Float_t ptspectra=h->GetBinCenter(ibin);
1571 if(ptspectra<tcutsdata->GetPtTOFMatching())
1573 h->SetBinContent(ibin,(h->GetBinContent(ibin)/factor));
1582 void TOFPIDsignalmatchingApply(TH1* h, Float_t factor)
1586 for(Int_t ibin=1;ibin<h->GetNbinsX();ibin++)
1588 Float_t ptspectra=h->GetBinCenter(ibin);
1589 if(ptspectra<tcutsdata->GetPtTOFMatching())
1591 h->SetBinContent(ibin,(h->GetBinContent(ibin)*factor));
1596 void CalculateDoubleCounts(TH1* doubleconunts,TH1** rawspectra,Int_t ipar, Bool_t dataflag)
1600 tmphist=(TH2F*)managerdata->GetGenMulvsRawMulHistogram("hHistDoubleCounts");
1602 tmphist=(TH2F*)managermc->GetGenMulvsRawMulHistogram("hHistDoubleCounts");
1606 TH1F* tmphist1=(TH1F*)rawspectra[ipar]->Clone("test");
1607 tmphist1->Add(rawspectra[ipar+3]);
1608 doubleconunts->Add(tmphist->ProjectionX(doubleconunts->GetName(),1,1));
1610 doubleconunts->Add(tmphist->ProjectionX("pi+k",2,2));
1612 doubleconunts->Add(tmphist->ProjectionX("pi+p",3,3));
1614 doubleconunts->Add(tmphist->ProjectionX("k+p",4,4));
1615 doubleconunts->Divide(doubleconunts,tmphist1,1,1,"B");