small changes in analysing macros
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / AnalysisBoth.C
CommitLineData
88d02ae4 1class AliSpectraBothHistoManager;
2class AliSpectraBothEventCuts;
3class AliSpectraBothTrackCuts;
4TString Charge[]={"Pos","Neg"};
5TString Sign[]={"Plus","Minus"};
6TString Particle[]={"Pion","Kaon","Proton"};
7AliSpectraBothHistoManager* managerdata=0x0;
8AliSpectraBothEventCuts* ecutsdata=0x0;
9AliSpectraBothTrackCuts* tcutsdata=0x0;
10
11AliSpectraBothHistoManager* managermc=0x0;
12AliSpectraBothEventCuts* ecutsmc=0x0;
13AliSpectraBothTrackCuts* tcutsmc=0x0;
14
15Float_t TOFMatchingScalling[2]={-1,-1};
16Int_t Color[3]={1,2,4};
17Int_t Marker[6]={20,21,22,24,25,26};
18Double_t Range[3]={0.3,0.3,0.5}; // LowPt range for pi k p
19
20enum ECharge_t {
21 kPositive,
22 kNegative,
23 kNCharges
24};
25
26Bool_t OpenFile(TString dirname, TString outputname, Bool_t mcflag);
27void AnalysisBoth (Bool_t dca=kTRUE,TString outdate, TString outnamedata, TString outnamemc="" )
28{
29 TH1::AddDirectory(kFALSE);
30 gSystem->Load("libCore.so");
31 gSystem->Load("libPhysics.so");
32 gSystem->Load("libTree");
33 gSystem->Load("libMatrix");
34 gSystem->Load("libSTEERBase");
35 gSystem->Load("libESD");
36 gSystem->Load("libAOD");
37 gSystem->Load("libANALYSIS");
38 gSystem->Load("libOADB");
39 gSystem->Load("libANALYSISalice");
40 gSystem->Load("libTENDER");
41 gSystem->Load("libCORRFW");
42 gSystem->Load("libPWGTools");
43 gSystem->Load("libPWGLFspectra");
44
45 gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/QAPlotsBoth.C");
46 Double_t mass[3];
47 mass[0] = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
48 mass[1] = TDatabasePDG::Instance()->GetParticle("K+")->Mass();
49 mass[2] = TDatabasePDG::Instance()->GetParticle("proton")->Mass();
50
51 AliESDtrackCuts* esdtrackcuts= AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE);
52 TString formulastring(esdtrackcuts->GetMaxDCAToVertexXYPtDep());
53 formulastring.ReplaceAll("pt","x");
54 TFormula* dcacutxy=new TFormula("dcacutxy",formulastring.Data());
55
56 TList* lout=new TList();
57
58
59 TString indirname=Form("/output/train%s",outdate.Data());
60 //TString indirname("/output/train24072012");
61 if(outnamemc.Length()==0)
62 outnamemc=outnamedata;
63 cout<<indirname.Data()<<" "<<outnamemc.Data()<<endl;
64 // Do the job
65
66
67 OpenFile(indirname,outnamemc,true);
68 OpenFile(indirname,outnamedata,false);
69 if(!managermc||!managerdata)
70 {
71 cout<<managermc<<" "<<managerdata<<endl;
72 return;
73 }
74 TH1F* rawspectradata[6];
75 TH1F* rawspectramc[6];
76 TH1F* MCTruth[6];
77 TH1F* eff[6];
78 TH1F* contallMC[6];
79 TH1F* contPID[6];
80 TH1F* contWD[6];
81 TH1F* contMat[6];
82 TH1F* confinal[6];
83
84 TH1F* contfit[12];
85 TH1F* contWDfit[12];
86 TH1F* contMatfit[12];
87 TH1F* primaryfit[12];
88
89
90
91 TH1F* spectra[6];
92 TH1F* spectraLeonardo[6];
93
94 TH1F* corrLeonardo[6];
95 //GetSpectra(managerdata,rawspectradata,true);
96 //GetSpectra(managermc,rawspectramc,true,true);
97
98 GetPtHistFromPtDCAhisto("hHistPtRecSigma","SpectraMC",managermc,rawspectramc,dcacutxy);
99 GetPtHistFromPtDCAhisto("hHistPtRecSigma","SpectraDATA",managerdata,rawspectradata,dcacutxy);
100 GetPtHistFromPtDCAhisto("hHistPtRecTruePrimary","eff",managermc,eff,dcacutxy);
101 GetPtHistFromPtDCAhisto("hHistPtRecTrue","conPID",managermc,contPID,dcacutxy);
102 GetPtHistFromPtDCAhisto("hHistPtRecSigmaSecondaryWeakDecay","conWD",managermc,contWD,dcacutxy);
103 GetPtHistFromPtDCAhisto("hHistPtRecSigmaSecondaryMaterial","conMat",managermc,contMat,dcacutxy);
104
105
106// Double_t neventsdata = ecutsdata->NumberOfEvents();
107 Double_t neventsmcall = ecutsmc->NumberOfEvents();
108 Double_t neventsdata = ecutsdata->NumberOfPhysSelEvents();
109 Double_t neventsmc = ecutsmc->NumberOfPhysSelEvents();
110
111 GetMCTruth(MCTruth);
112
113
114
115 TH1F* allgen=((TH1F*)managermc->GetPtHistogram1D("hHistPtGen",1,1))->Clone();
116 allgen->SetName("AllGen");
117 TH1F* allrecMC=GetOneHistFromPtDCAhisto("hHistPtRec","rawallMC",managermc,dcacutxy);
118 TH1F* alleff=GetOneHistFromPtDCAhisto("hHistPtRecPrimary","effall",managermc,dcacutxy);
119 TH1F* allrecdata=GetOneHistFromPtDCAhisto("hHistPtRec","rawalldata",managerdata,dcacutxy);
120
121
122
123
124 TH1F* spectraall=(TH1F*)allrecdata->Clone("recNch");
125 TH1F* contall=(TH1F*)allrecMC->Clone("contall");
126 contall->Add(alleff,-1);
127 alleff->Divide(alleff,allgen,1,1,"B");
128 contall->Divide(contall,allgen,1,1,"B");
129
130 GetCorrectedSpectra(spectraall,allrecdata,alleff,contall);
131 Divideby2pipt(spectraall);
132
133 allrecdata->Scale(1./neventsdata,"width");
134 allgen->Scale(1./neventsmcall,"width");
135 allrecMC->Scale(1./neventsmc,"width");
136 spectraall->Scale(1./neventsdata,"width");
137
138
139 lout->Add(allgen);
140 lout->Add(allrecMC);
141 lout->Add(alleff);
142 lout->Add(allrecdata);
143 lout->Add(spectraall);
144 lout->Add(contall);
145 for (int i=0;i<6;i++)
146 {
147
148
149 TString tmpname(rawspectramc[i]->GetTitle());
150 tmpname.ReplaceAll("SpectraMC","%s");
151 contallMC[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contallMC"));
152 contfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contfit"));
153 contWDfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contWDfit"));
154 contMatfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contMatfit"));
155 primaryfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"primaryfit"));
156
157 contfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contfitonMC"));
158 contWDfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contWDfitonMC"));
159 contMatfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contMatfitonMC"));
160 primaryfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"primaryfitMC"));
161
162 contfit[i]->Reset();
163 contWDfit[i]->Reset();
164 contMatfit[i]->Reset();
165 primaryfit[i]->Reset();
166
167
168 contfit[i+6]->Reset();
169 contWDfit[i+6]->Reset();
170 contMatfit[i+6]->Reset();
171 primaryfit[i+6]->Reset();
172
173 SetBintoOne(primaryfit[i]);
174 SetBintoOne(primaryfit[i+6]);
175 spectra[i]=(TH1F*)rawspectradata[i]->Clone(Form(tmpname.Data(),"SpectraFinal"));
176
177 spectraLeonardo[i]=(TH1F*)rawspectradata[i]->Clone(Form(tmpname.Data(),"SpectraFinalLeonardo"));
178 corrLeonardo[i]=(TH1F*)MCTruth[i]->Clone(Form(tmpname.Data(),"CorrFactLeonardo"));
179
180 corrLeonardo[i]->Divide(corrLeonardo[i],rawspectramc[i],1,1,"B");
181
182
183
184 contallMC[i]->Add(eff[i],-1.0);
185 RecomputeErrors(contallMC[i]);
186 contallMC[i]->Sumw2();
187 contallMC[i]->Divide(contallMC[i],rawspectramc[i],1,1,"B");
188
189 eff[i]->Divide(eff[i],MCTruth[i],1,1,"B");
190
191
192 contPID[i]->Sumw2();
193 rawspectramc[i]->Sumw2();
194 contPID[i]->Add(contPID[i],rawspectramc[i],-1,1);
195 RecomputeErrors(contPID[i]);
196 contPID[i]->ResetStats();
197 contPID[i]->Sumw2();
198 contPID[i]->Divide(contPID[i],rawspectramc[i],1,1,"B");
199
200 confinal[i]=(TH1F*)contPID[i]->Clone(Form(tmpname.Data(),"confinal"));
201
202
203 contWD[i]->Divide(contWD[i],rawspectramc[i],1,1,"B");
204 contMat[i]->Divide(contMat[i],rawspectramc[i],1,1,"B");
205
206
207
208 rawspectradata[i]->Scale(1./neventsdata,"width");
209 rawspectramc[i]->Scale(1./neventsmc,"width");
210 MCTruth[i]->Scale(1./neventsmcall,"width");
211 spectraLeonardo[i]->Scale(1./neventsdata,"width");
212
213
214
215 lout->Add(rawspectradata[i]);
216 lout->Add(rawspectramc[i]);
217 lout->Add(MCTruth[i]);
218 lout->Add(eff[i]);
219 lout->Add(contallMC[i]);
220 lout->Add(contPID[i]);
221 lout->Add(contWD[i]);
222 lout->Add(contMat[i]);
223 lout->Add(contfit[i]);
224 lout->Add(contWDfit[i]);
225 lout->Add(contMatfit[i]);
226 lout->Add(contfit[i+6]);
227 lout->Add(contWDfit[i+6]);
228 lout->Add(contMatfit[i+6]);
229 lout->Add(spectra[i]);
230 lout->Add(spectraLeonardo[i]);
231 lout->Add(confinal[i]);
232 }
233 TFile* fout=new TFile(Form("./results/ResMY_%s_%s.root",outnamemc.Data(),outdate.Data()),"RECREATE");
234 if (dca)
235 DCACorrectionMarek(managerdata,managermc,dcacutxy,fout,contfit,contWDfit,contMatfit,primaryfit);
236 for (int i=0;i<6;i++)
237 {
238 if(dca)
239 {
240 confinal[i]->Add(contfit[i]);
241 GetCorrectedSpectra(spectra[i],rawspectradata[i],eff[i],confinal[i]);
242 }
243 else
244 {
245 GetCorrectedSpectra(spectra[i],rawspectradata[i],eff[i],contallMC[i]);
246 }
247 GetCorrectedSpectraLeonardo(spectraLeonardo[i],corrLeonardo[i],primaryfit[i],primaryfit[i+6]);
248 CleanHisto(spectra[i],-1,100,contPID[i]);
249 CleanHisto(spectraLeonardo[i],-1,100,contPID[i]);
250 }
251
252 GFCorrection(spectra,tcutsdata->GetPtTOFMatching());
253 GFCorrection(spectraLeonardo,tcutsdata->GetPtTOFMatching());
254 MatchingTOFEff(spectra,lout);
255 MatchingTOFEff(spectraLeonardo);
256 TH1F* allch=GetSumAllCh(spectra,mass);
257 lout->Add(allch);
258
259// lout->ls();
260 fout->cd();
261 TList* listqa=new TList();
33d374d5 262 TList* canvaslist=new TList();
263 Float_t vertexcorrection=QAPlotsBoth(managerdata,managermc,ecutsdata,ecutsmc,tcutsdata,tcutsmc,listqa,canvaslist);
88d02ae4 264 cout<<" VTX corr="<<vertexcorrection<<endl;
265 for (int i=0;i<6;i++)
266 {
267 spectra[i]->Scale(vertexcorrection);
268 spectraLeonardo[i]->Scale(vertexcorrection);
269 }
270 allch->Scale(vertexcorrection);
271 spectraall->Scale(vertexcorrection/1.6);
272
273 //spectraall->Scale(1.0/1.6);
274 lout->Write("output",TObject::kSingleKey);
275 listqa->Write("outputQA",TObject::kSingleKey);
33d374d5 276 canvaslist->Write("outputcanvas",TObject::kSingleKey);
277
88d02ae4 278 fout->Close();
279
280}
281
282Bool_t OpenFile(TString dirname,TString outputname, Bool_t mcflag)
283{
284
285 TString nameFile = Form("./%s/AnalysisResults%s.root",dirname.Data(),(mcflag?"MC":"DATA"));
286 TFile *file = TFile::Open(nameFile.Data());
287 if(!file)
288 {
289 cout<<"no file"<<endl;
290 return false;
291 }
292 TString sname=Form("OutputBothSpectraTask_%s_%s",(mcflag?"MC":"Data"),outputname.Data());
293 file->ls();
294 TDirectoryFile *dir=(TDirectoryFile*)file->Get(sname.Data());
295 if(!dir)
296 {
297 // cout<<"no dir "<<sname.Data()<<endl;
298 sname=Form("OutputAODSpectraTask_%s_%s",(mcflag?"MC":"Data"),outputname.Data());
299 // cout<<"trying "<<sname.Data()<<endl;
300 dir=(TDirectoryFile*)file->Get(sname.Data());
301 if(!dir)
302 {
303 cout<<"no dir "<<sname.Data()<<endl;
304 return false;
305 }
306 }
307 cout << " -- Info about " <<(mcflag?"MC":"DATA") <<" -- "<< endl;
308 if(mcflag)
309 {
310 managermc= (AliSpectraBothHistoManager*) dir->Get("SpectraHistos");
311 ecutsmc = (AliSpectraBothEventCuts*) dir->Get("Event Cuts");
312 tcutsmc = (AliSpectraBothTrackCuts*) dir->Get("Track Cuts");
313 ecutsmc->PrintCuts();
314 tcutsmc->PrintCuts();
315 if(!managermc||!ecutsmc||!tcutsmc)
316 return false;
317 }
318 else
319 {
320 managerdata= (AliSpectraBothHistoManager*) dir->Get("SpectraHistos");
321 ecutsdata = (AliSpectraBothEventCuts*) dir->Get("Event Cuts");
322 tcutsdata = (AliSpectraBothTrackCuts*) dir->Get("Track Cuts");
323 ecutsdata->PrintCuts();
324 tcutsdata->PrintCuts();
325 if(!managerdata||!ecutsdata||!tcutsdata)
326 return false;
327 }
328 return true;
329}
330
331 void GetMCTruth(TH1F** MCTruth)
332 {
333 for(Int_t icharge=0;icharge<2;icharge++)
334 {
335 for(Int_t ipart=0;ipart<3;ipart++)
336 {
337 Int_t index=ipart+3*icharge;
338 TString hname=Form("hHistPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
339 MCTruth[index]=(TH1F*)((TH1F*)managermc->GetPtHistogram1D(hname.Data(),1,1))->Clone();
340 MCTruth[index]->SetName(Form("MCTruth_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
341 MCTruth[index]->SetTitle(Form("MCTruth_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
342 MCTruth[index]->Sumw2();
343 }
344 }
345}
346
347TH1F* GetOneHistFromPtDCAhisto(TString name,TString hnameout,AliSpectraBothHistoManager* hman,TFormula* dcacutxy)
348{
349 histo =(TH1F*)((TH1F*) hman->GetPtHistogram1D(name.Data(),-1,-1))->Clone();
350 histo->SetName(hnameout.Data());
351 histo->SetTitle(hnameout.Data());
352
353 if(dcacutxy)
354 {
355 for(int ibin=1;ibin<histo->GetNbinsX();ibin++)
356 {
357 Double_t lowedge=histo->GetBinLowEdge(ibin);
358 Float_t cut=dcacutxy->Eval(lowedge);
359 TH1F* dcahist=(TH1F*)hman->GetDCAHistogram1D(name.Data(),lowedge,lowedge));
360 Float_t inyield=dcahist->Integral(dcahist->GetXaxis()->FindBin(-1.0*cut),dcahist->GetXaxis()->FindBin(cut));
361 cout<<"corr data "<<histo->GetBinContent(ibin)<<" "<<inyield<<" "<<dcahist->Integral()<<" "<<hnameout.Data()<<endl;
362 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;
363 histo->SetBinContent(ibin,inyield);
364 histo->SetBinError(ibin,TMath::Sqrt(inyield));
365 }
366 }
367 histo->Sumw2();
368 return histo;
369}
370
371
372
373
374void GetPtHistFromPtDCAhisto(TString hnamein, TString hnameout, AliSpectraBothHistoManager* hman,TH1F** histo,TFormula* dcacutxy)
375{
376 Float_t min[3]={0.3,0.3,0.4};
377 Float_t max[3]={1.5,1.2,2.2};
378 for(Int_t icharge=0;icharge<2;icharge++)
379 {
380 for(Int_t ipart=0;ipart<3;ipart++)
381 {
382 Int_t index=ipart+3*icharge;
383 //TString hname=Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data());
384 printf("Getting %s",hnamein.Data());
385 TString nameinfinal=Form("%s%s%s",hnamein.Data(),Particle[ipart].Data(),Sign[icharge].Data());
386 TString nameoutfinal=Form("%s%s%s",hnameout.Data(),Particle[ipart].Data(),Sign[icharge].Data());
387
388
389 histo[index]=GetOneHistFromPtDCAhisto(nameinfinal,nameoutfinal,hman,dcacutxy);
390 /*histo[index] =(TH1F*)((TH1F*) hman->GetPtHistogram1D(Form("%s%s%s",hnamein.Data(),Particle[ipart].Data(),Sign[icharge].Data()),-1,-1))->Clone();
391 histo[index]->SetName(Form("%s%s%s",hnameout.Data(),Particle[ipart].Data(),Sign[icharge].Data()));
392 histo[index]->SetTitle(Form("%s%s%s",hnameout.Data(),Particle[ipart].Data(),Sign[icharge].Data()));
393
394 if(dcacutxy)
395 {
396 for(int ibin=1;ibin<histo[index]->GetNbinsX();ibin++)
397 {
398 Double_t lowedge=histo[index]->GetBinLowEdge(ibin);
399 Float_t cut=dcacutxy->Eval(lowedge);
400 TH1F* dcahist=(TH1F*)hman->GetDCAHistogram1D(Form("%s%s%s",hnamein.Data(),Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge));
401 Float_t inyield=dcahist->Integral(dcahist->GetXaxis()->FindBin(-1.0*cut),dcahist->GetXaxis()->FindBin(cut));
402 cout<<"corr data "<<histo[index]->GetBinContent(ibin)<<" "<<inyield<<" "<<dcahist->Integral()<<endl;
403 histo[index]->SetBinContent(ibin,inyield);
404 histo[index]->SetBinError(ibin,TMath::Sqrt(inyield));
405 }
406 }*/
407 CleanHisto(histo[index],min[ipart],max[ipart]);
408 // histo[index]->Sumw2();
409 }
410 }
411}
412void CleanHisto(TH1F* h, Float_t minV, Float_t maxV,TH1* contpid=0x0)
413{
414 for (int i=0;i<=h->GetNbinsX();i++)
415 {
416 if(h->GetXaxis()->GetBinCenter(i)<minV||h->GetXaxis()->GetBinCenter(i)>maxV)
417 {
418 h->SetBinContent(i,0);
419 h->SetBinError(i,0);
420 }
421 if(contpid)
422 {
423 if(contpid->GetBinContent(i)>0.2)
424 {
425 h->SetBinContent(i,0);
426 h->SetBinError(i,0);
427 }
428 }
429 }
430}
431
432
433void DCACorrectionMarek(AliSpectraBothHistoManager* hman_data, AliSpectraBothHistoManager* hman_mc,TFormula* fun,TFile *fout,TH1F** hcon,TH1F** hconWD,TH1F** hconMat,TH1F** hprimary)
434{
435 printf("\n\n-> DCA Correction");
436 Double_t FitRange[2]={-3.,3.};
437 Double_t CutRange[2]={-1.0,1.0};
438 Double_t minptformaterial[6]={0.0,100.0,0.0,0.0,100.0,0.0};
439 Double_t maxptformaterial[6]={0.0,-100.0,1.3,0.0,-100.0,0.0};
440 Bool_t usefit[3]={true,false,true};
441 Double_t maxbinforfit[6]={1.5,0,2.0,1.5,0,2.0};
442 Printf("\DCACorr");
443 // TH1F *hcorrection[2];
444 /* TCanvas *ccorrection=new TCanvas("DCAcorrection","DCAcorrection",700,500);
445 TCanvas *cRatiocorrection=new TCanvas("DCARatiocorrection","DCARatiocorrection",700,500);
446 cRatiocorrection->Divide(2,1);
447 ccorrection->Divide(2,1);*/
448 TString sample[2]={"data","mc"};
449 ofstream debug("debugDCA.txt");
450 TList* listofdcafits=new TList();
451 for(Int_t icharge=0;icharge<2;icharge++)
452 {
453 for(Int_t ipart=0;ipart<3;ipart++)
454 {
455 Int_t index=ipart+3*icharge;
456 for(Int_t isample=0;isample<2;isample++)
457 {
458
459 /*hcorrection[isample]=(TH1F*)Spectra[index]->Clone();
460 hcorrection[isample]->Reset("all");*/
461 for(Int_t ibin_data=7;ibin_data<40;ibin_data++)
462 {
463
464 Double_t lowedge=hcon[index]->GetBinLowEdge(ibin_data);
465 Double_t binwidth=hcon[index]->GetBinWidth(ibin_data);
466 debug<<"NEW "<<Particle[ipart].Data()<<" "<<Sign[icharge].Data()<<" "<<lowedge<<endl;
467 CutRange[1]=fun->Eval(lowedge);
468 CutRange[0]=-1.0*CutRange[1];
469 debug<<"cut "<<CutRange[1]<<" "<<CutRange[0]<<endl;
470 Bool_t useMaterial=kFALSE;
471 cout<<"try to fit "<<lowedge<<" "<<maxbinforfit[index]<<endl;
472 if(lowedge>maxbinforfit[index])
473 continue;
474 if(lowedge>minptformaterial[index]&&lowedge<maxptformaterial[index])
475 useMaterial=kTRUE;
476
477 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);
478 if(isample==0)
479 TH1F *hToFit =(TH1F*) ((TH1F*)hman_data->GetDCAHistogram1D(Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
480 if(isample==1)
481 TH1F *hToFit =(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
482 debug<<Particle[ipart].Data()<<" "<<Sign[icharge].Data()<<" "<<lowedge<<endl;
483 TH1F *hmc1=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
484 TH1F *hmc2=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaSecondaryWeakDecay%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
485 TH1F *hmc3=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaSecondaryMaterial%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
486 Double_t minentries=1;
487 debug<<" Entries "<<isample<<" "<<hToFit->GetEntries()<<" "<<hmc1->GetEntries()<<" "<<hmc2->GetEntries()<<" "<<hmc3->GetEntries()<<endl;
488 debug<<"2 Entries "<<isample<<" "<<hToFit->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))<<" "<<hmc1->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))<<" "<<hmc2->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))<<" "<<hmc3->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))<<endl;
489 debug<< CutRange[0]<<" "<<CutRange[1]<<" "<<lowedge<<endl;
490 if(hToFit->GetEntries()<=minentries || hmc1->GetEntries()<=minentries || hmc2->GetEntries()<=minentries || hmc3->GetEntries()<=minentries)
491 continue;
492 //Data and MC can have different stat
493 hToFit->Sumw2();
494 hmc1->Sumw2();
495 hmc2->Sumw2();
496 hmc3->Sumw2();
497
498 Float_t corrforrebinning[4]={1.0,1.0,1.0,1.0};
499
500
501 if(hmc3->GetNbinsX()>300)
502 {
503
504 corrforrebinning[0]=hToFit->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]));
505 corrforrebinning[1]=hmc1->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]));
506 corrforrebinning[2]=hmc2->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]));
507 corrforrebinning[3]=hmc3->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]));
508
509 hToFit->Rebin(30);
510 hmc1->Rebin(30);
511 hmc2->Rebin(30);
512 hmc3->Rebin(30);
513 if(hToFit->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))>0.0)
514 corrforrebinning[0]=corrforrebinning[0]/hToFit->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]));
515 else
516 corrforrebinning[0]=1.0;
517 if(hmc1->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))>0.0)
518 corrforrebinning[1]=corrforrebinning[1]/hmc1->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]));
519 else
520 corrforrebinning[1]=1.0;
521 if(hmc2->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))>0.0)
522 corrforrebinning[2]=corrforrebinning[2]/hmc2->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]));
523 else
524 corrforrebinning[2]=1.0;
525 if(hmc3->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))>0.0)
526 corrforrebinning[3]=corrforrebinning[3]/hmc3->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]));
527 else
528 corrforrebinning[3]=1.0;
529
530 debug<<" cor bin "<<corrforrebinning[0]<<" "<<corrforrebinning[1]<<" "<<corrforrebinning[2]<<" "<<corrforrebinning[3]<<endl;
531
532
533 }
534
535 cDCA->cd();
536 gPad->SetGridy();
537 gPad->SetGridx();
538 gPad->SetLogy();
539
540 TObjArray *mc=0x0;
541 if(useMaterial)
542 mc = new TObjArray(3); // MC histograms are put in this array
543 else
544 mc = new TObjArray(2);
545 mc->Add(hmc1);
546 mc->Add(hmc2);
547 if(useMaterial)
548 mc->Add(hmc3);
549 TFractionFitter* fit = new TFractionFitter(hToFit,mc); // initialise
550 fit->Constrain(0,0.0,1.0); // constrain fraction 1 to be between 0 and 1
551 fit->Constrain(1,0.0,1.0); // constrain fraction 1 to be between 0 and 1
552 if(useMaterial)
553 fit->Constrain(2,0.0,1.0); // constrain fraction 1 to be between 0 and 1
554 fit->SetRangeX(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1]));
555 hToFit->GetXaxis()->SetRange(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1]));
556 hToFit->SetTitle(Form("DCA distr - %s %s %s %lf",sample[isample].Data(),Particle[ipart].Data(),Sign[icharge].Data(),lowedge));
557 Int_t status = fit->Fit(); // perform the fit
558 cout << "fit status: " << status << endl;
559 debug<<"fit status: " << status << endl;
560
561 if (status == 0 && usefit[ipart])
562 { // check on fit status
563 TH1F* result = (TH1F*) fit->GetPlot();
564 TH1F* PrimMCPred=(TH1F*)fit->GetMCPrediction(0);
565 TH1F* secStMCPred=(TH1F*)fit->GetMCPrediction(1);
566
567 TH1F* secMCPred=0x0;
568 if(useMaterial)
569 secMCPred=(TH1F*)fit->GetMCPrediction(2);
570
571 Double_t v1=0,v2=0,v3=0;
572 Double_t ev1=0,ev2=0,ev3=0;
573 Double_t cov=0.0;
574 //first method, use directly the fit result
575 fit->GetResult(0,v1,ev1);
576 fit->GetResult(1,v2,ev2);
577 if(useMaterial)
578 {
579 fit->GetResult(2,v3,ev3);
580 fit->GetFitter()->GetCovarianceMatrixElement(1,2);
581 }
582 debug<<v1<<" "<<ev1<<" "<<v2<<" "<<ev2<<" "<<v3<<" "<<ev3<<" "<<endl;
583
584
585
586 Float_t normalizationdata=hToFit->Integral(hToFit->GetXaxis()->FindBin(CutRange[0]),hToFit->GetXaxis()->FindBin(CutRange[1]))/hToFit->Integral(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1]));
587
588 Float_t normalizationmc1=hmc1->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))/hmc1->Integral(hmc1->GetXaxis()->FindBin(FitRange[0]),hmc1->GetXaxis()->FindBin(FitRange[1]))/normalizationdata;
589 Float_t normalizationmc2=hmc2->Integral(hmc2->GetXaxis()->FindBin(CutRange[0]),hmc2->GetXaxis()->FindBin(CutRange[1]))/hmc2->Integral(hmc1->GetXaxis()->FindBin(FitRange[0]),hmc2->GetXaxis()->FindBin(FitRange[1]))/normalizationdata;
590 Float_t normalizationmc3=hmc3->Integral(hmc3->GetXaxis()->FindBin(CutRange[0]),hmc3->GetXaxis()->FindBin(CutRange[1]))/hmc3->Integral(hmc1->GetXaxis()->FindBin(FitRange[0]),hmc3->GetXaxis()->FindBin(FitRange[1]))/normalizationdata;
591 debug<<"After Nor"<<endl;
592 debug<<v1*normalizationmc1<<" "<<ev1*normalizationmc1<<" "<<v2*normalizationmc2<<" "<<ev2*normalizationmc2<<" "<<v3*normalizationmc3<<" "<<ev3*normalizationmc3<<" "<<endl;
593 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;
594 Float_t normalizationdata1=result->Integral(result->GetXaxis()->FindBin(CutRange[0]),result->GetXaxis()->FindBin(CutRange[1]))/result->Integral(result->GetXaxis()->FindBin(FitRange[0]),result->GetXaxis()->FindBin(FitRange[1]));
595
596
597 normalizationdata1*=corrforrebinning[0];
598
599
600 Float_t normalizationmc11=PrimMCPred->Integral(PrimMCPred->GetXaxis()->FindBin(CutRange[0]),PrimMCPred->GetXaxis()->FindBin(CutRange[1]))/PrimMCPred->Integral(PrimMCPred->GetXaxis()->FindBin(FitRange[0]),PrimMCPred->GetXaxis()->FindBin(FitRange[1]))/normalizationdata1;
601 Float_t normalizationmc21=secStMCPred->Integral(secStMCPred->GetXaxis()->FindBin(CutRange[0]),secStMCPred->GetXaxis()->FindBin(CutRange[1]))/secStMCPred->Integral(secStMCPred->GetXaxis()->FindBin(FitRange[0]),secStMCPred->GetXaxis()->FindBin(FitRange[1]))/normalizationdata1;
602 Float_t normalizationmc31=0;
603 if(useMaterial)
604 normalizationmc31=secMCPred->Integral(secMCPred->GetXaxis()->FindBin(CutRange[0]),secMCPred->GetXaxis()->FindBin(CutRange[1]))/secMCPred->Integral(secMCPred->GetXaxis()->FindBin(FitRange[0]),secMCPred->GetXaxis()->FindBin(FitRange[1]))/normalizationdata1;
605
606 normalizationmc11*=corrforrebinning[1];
607 normalizationmc21*=corrforrebinning[2];
608 normalizationmc31*=corrforrebinning[3];
609
610 debug<<"After Nor 2"<<endl;
611 debug<<v1*normalizationmc11<<" "<<ev1*normalizationmc11<<" "<<v2*normalizationmc21<<" "<<ev2*normalizationmc21<<" "<<v3*normalizationmc31<<" "<<ev3*normalizationmc31<<endl;
612
613 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;
614
615 debug<<CutRange[0]<<" "<<CutRange[1]<<endl;
616 debug<<" Entries "<<isample<<" "<<hToFit->GetEntries()<<" "<<hmc1->GetEntries()<<" "<<hmc2->GetEntries()<<" "<<hmc3->GetEntries()<<endl;
617 debug<<"2 Entries "<<isample<<" "<<hToFit->Integral(result->GetXaxis()->FindBin(CutRange[0]),result->GetXaxis()->FindBin(CutRange[1]))<<" "<<hmc1->Integral(result->GetXaxis()->FindBin(CutRange[0]),result->GetXaxis()->FindBin(CutRange[1]))<<" "<<hmc2->Integral(result->GetXaxis()->FindBin(CutRange[0]),result->GetXaxis()->FindBin(CutRange[1]))<<" "<<hmc3->Integral(result->GetXaxis()->FindBin(CutRange[0]),result->GetXaxis()->FindBin(CutRange[1]))<<endl;
618
619
620 hconWD[index+6*isample]->SetBinContent(ibin_data,v2*normalizationmc21);
621 hconWD[index+6*isample]->SetBinError(ibin_data,ev2*normalizationmc21);
622 hconMat[index+6*isample]->SetBinContent(ibin_data,v3*normalizationmc31);
623 hconMat[index+6*isample]->SetBinError(ibin_data,ev3*normalizationmc31);
624 hprimary[index+6*isample]->SetBinContent(ibin_data,v1*normalizationmc11);
625 hprimary[index+6*isample]->SetBinContent(ibin_data,v1*normalizationmc11);
626 if(useMaterial)
627 {
628 hcon[index+6*isample]->SetBinContent(ibin_data,v2*normalizationmc21+v3*normalizationmc31);
629 hcon[index+6*isample]->SetBinError(ibin_data,TMath::Sqrt(ev2*ev2*normalizationmc21*normalizationmc21+ev3*ev3*normalizationmc31*normalizationmc31+cov*normalizationmc31*normalizationmc21));
630 }
631 else
632 {
633 hcon[index+6*isample]->SetBinContent(ibin_data,v2*normalizationmc21);
634 hcon[index+6*isample]->SetBinError(ibin_data,ev2*normalizationmc21);
635 }
636
637
638
639 //Drawing section
640 result->Scale(1.0/result->Integral(result->GetXaxis()->FindBin(FitRange[0]),result->GetXaxis()->FindBin(FitRange[1])));
641 hToFit->Scale(1.0/hToFit->Integral(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1])));
642 PrimMCPred->Scale(v1/PrimMCPred->Integral(PrimMCPred->GetXaxis()->FindBin(FitRange[0]),PrimMCPred->GetXaxis()->FindBin(FitRange[1])));
643 secStMCPred->Scale(v2/secStMCPred->Integral(secStMCPred->GetXaxis()->FindBin(FitRange[0]),secStMCPred->GetXaxis()->FindBin(FitRange[1])));
644 if(useMaterial)
645 secMCPred->Scale(v3/secMCPred->Integral(secMCPred->GetXaxis()->FindBin(FitRange[0]),secMCPred->GetXaxis()->FindBin(FitRange[1])));
646
647 result->SetLineColor(kBlack);
648 PrimMCPred->SetLineColor(kGreen+2);
649 secStMCPred->SetLineColor(kRed);
650
651 hToFit->SetMinimum(0.0001);
652 hToFit->DrawClone("E1x0");
653 result->SetTitle("Fit result");
654 result->DrawClone("lhistsame");
655 PrimMCPred->DrawClone("lhistsame");
656 secStMCPred->DrawClone("lhistsame");
657 if(useMaterial)
658 {
659 secMCPred->SetLineColor(kBlue);
660 secMCPred->DrawClone("lhistsame");
661 }
662 }
663 else
664 {
665 hconWD[index+6*isample]->SetBinContent(ibin_data,0.0);
666 hconWD[index+6*isample]->SetBinError(ibin_data,0.0);
667 hconMat[index+6*isample]->SetBinContent(ibin_data,0.0);
668 hconMat[index+6*isample]->SetBinError(ibin_data,0.0);
669 hcon[index+6*isample]->SetBinContent(ibin_data,0.0);
670 hcon[index+6*isample]->SetBinError(ibin_data,0.0);
671 hprimary[index+6*isample]->SetBinContent(ibin_data,1.0);
672 hprimary[index+6*isample]->SetBinError(ibin_data,0.0);
673 }
674 listofdcafits->Add(cDCA);
675
676 //cDCA->Write();
677 delete hToFit;
678 }
679
680
681 }
682
683 }
684 }
685 fout->cd();
686 listofdcafits->Write("DCAfits",TObject::kSingleKey);
687}
688
689void RecomputeErrors(TH1* h)
690{
691 for (int i=0; i<=h->GetXaxis()->GetNbins(); i++)
692 h->SetBinError(i,TMath::Sqrt(h->GetBinContent(i)));
693 h->Sumw2();
694}
695void SetBintoOne(TH1* h)
696{
697 for (int i=0;i<=h->GetXaxis()->GetNbins(); i++)
698 {
699 h->SetBinContent(i,1);
700 h->SetBinError(i,0);
701 }
702}
703
704
705void GetCorrectedSpectra(TH1F* corr,TH1F* raw,TH1F* eff, TH1F* con)
706{
707 for (int i=0;i<=corr->GetXaxis()->GetNbins(); i++)
708 {
709 corr->SetBinContent(i,1);
710 corr->SetBinError(i,0);
711 }
712 corr->Sumw2();
713 corr->Add(con,-1);
714 corr->Sumw2();
715 corr->Divide(eff);
716 corr->Sumw2();
717 corr->Multiply(raw);
718 corr->Sumw2();
719}
720void GetCorrectedSpectraLeonardo(TH1F* spectra,TH1F* correction, TH1F* hprimaryData,TH1F* hprimaryMC)
721{
722 spectra->Sumw2();
723 spectra->Multiply(correction);
724 spectra->Sumw2();
725 hprimaryData->Sumw2();
726 spectra->Multiply(hprimaryData);
727 hprimaryMC->Sumw2();
728 spectra->Divide(hprimaryMC);
729}
730
731void GFCorrection(TH1F **Spectra,Float_t tofpt)
732{
733 //Geant/Fluka Correction
734 Printf("\nGF correction for Kaons");
735 //Getting GF For Kaons in TPC
736 TGraph *gGFCorrectionKaonPlus=new TGraph();
737 gGFCorrectionKaonPlus->SetName("gGFCorrectionKaonPlus");
738 gGFCorrectionKaonPlus->SetTitle("gGFCorrectionKaonPlus");
739 TGraph *gGFCorrectionKaonMinus=new TGraph();
740 gGFCorrectionKaonMinus->SetName("gGFCorrectionKaonMinus");
741 gGFCorrectionKaonMinus->SetTitle("gGFCorrectionKaonMinus");
742 TString fnameGeanFlukaK="GFCorrection/correctionForCrossSection.321.root";
743 TFile *fGeanFlukaK= new TFile(fnameGeanFlukaK.Data());
744 if (!fGeanFlukaK)
745 return;
746 TH1F *hGeantFlukaKPos=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles");
747 TH1F *hGeantFlukaKNeg=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles");
748 //getting GF func for Kaons with TOF
749 TF1 *fGFKPosTracking;
750 fGFKPosTracking = TrackingEff_geantflukaCorrection(3,kPositive);
751 TF1 *fGFKNegTracking;
752 fGFKNegTracking = TrackingEff_geantflukaCorrection(3,kNegative);
753 TF1 *fGFKPosMatching;
754 fGFKPosMatching = TOFmatchMC_geantflukaCorrection(3,kPositive);
755 TF1 *fGFKNegMatching;
756 fGFKNegMatching = TOFmatchMC_geantflukaCorrection(3,kNegative);
757 for(Int_t binK=0;binK<=Spectra[1]->GetNbinsX();binK++)
758 {
759 if(Spectra[1]->GetBinCenter(binK)<tofpt)
760 {//use TPC GeantFlukaCorrection
761 Float_t FlukaCorrKPos=hGeantFlukaKPos->GetBinContent(hGeantFlukaKPos->FindBin(Spectra[1]->GetBinCenter(binK)));
762 Float_t FlukaCorrKNeg=hGeantFlukaKNeg->GetBinContent(hGeantFlukaKNeg->FindBin(Spectra[4]->GetBinCenter(binK)));
763 // Printf("TPC Geant/Fluka: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPos,FlukaCorrKNeg);
764 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPos);
765 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNeg);
766 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPos);
767 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNeg);
768 gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),FlukaCorrKPos);
769 gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),FlukaCorrKNeg);
770 }
771 else
772 {
773 gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),0);
774 gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),0);
775 Float_t FlukaCorrKPosTracking=fGFKPosTracking->Eval(Spectra[1]->GetBinCenter(binK));
776 Float_t FlukaCorrKNegTracking=fGFKNegTracking->Eval(Spectra[1]->GetBinCenter(binK));
777 // Printf("TPC/TOF Geant/Fluka Tracking: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosTracking,FlukaCorrKNegTracking);
778 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosTracking);
779 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegTracking);
780 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosTracking);
781 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegTracking);
782 Float_t FlukaCorrKPosMatching=fGFKPosMatching->Eval(Spectra[1]->GetBinCenter(binK));
783 Float_t FlukaCorrKNegMatching=fGFKNegMatching->Eval(Spectra[1]->GetBinCenter(binK));
784 // Printf("TPC/TOF Geant/Fluka Matching: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosMatching,FlukaCorrKNegMatching);
785 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosMatching);
786 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegMatching);
787 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosMatching);
788 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegMatching);
789 }
790 }
791
792 //Geant Fluka for P in TPC
793 Printf("\nGF correction for Protons");
794 const Int_t kNCharge=2;
795 Int_t kPos=0;
796 Int_t kNeg=1;
797 TFile* fGFProtons = new TFile ("GFCorrection/correctionForCrossSection.root");
798 TH2D * hCorrFluka[kNCharge];
799 TH2D * hCorrFluka[2];
800 hCorrFluka[kPos] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionProtons");
801 hCorrFluka[kNeg] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionAntiProtons");
802 //getting GF func for Kaons with TPCTOF
803 TF1 *fGFpPosTracking;
804 fGFpPosTracking = TrackingEff_geantflukaCorrection(4,kPositive);
805 TF1 *fGFpNegTracking;
806 fGFpNegTracking = TrackingEff_geantflukaCorrection(4,kNegative);
807 TF1 *fGFpPosMatching;
808 fGFpPosMatching = TOFmatchMC_geantflukaCorrection(4,kPositive);
809 TF1 *fGFpNegMatching;
810 fGFpNegMatching = TOFmatchMC_geantflukaCorrection(4,kNegative);
811
812
813 Int_t nbins = Spectra[2]->GetNbinsX();
814
815 for(Int_t ibin = 0; ibin < nbins; ibin++)
816 {
817 if(Spectra[2]->GetBinCenter(ibin)<tofpt)
818 {//use TPC GeantFlukaCorrection
819 for(Int_t icharge = 0; icharge < kNCharge; icharge++)
820 {
821 Int_t nbinsy=hCorrFluka[icharge]->GetNbinsY();
822 Float_t pt = Spectra[2]->GetBinCenter(ibin);
823 Float_t minPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(1);
824 Float_t maxPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1);
825 if (pt < minPtCorrection)
826 pt = minPtCorrection+0.0001;
827 if (pt > maxPtCorrection)
828 pt = maxPtCorrection;
829 Float_t correction = hCorrFluka[icharge]->GetBinContent(1,hCorrFluka[icharge]->GetYaxis()->FindBin(pt));
830
831 if (correction > 0.0)
832 {// If the bin is empty this is a 0
833 // cout<<icharge<<" "<<ibin<<" "<<correction<<endl;
834 Spectra[icharge*3+2]->SetBinContent(ibin,Spectra[icharge*3+2]->GetBinContent(ibin)*correction);
835 Spectra[icharge*3+2]->SetBinError(ibin,Spectra[icharge*3+2]->GetBinError(ibin)*correction);
836 }
837 else if (Spectra[icharge*3+2]->GetBinContent(ibin) > 0.0)
838 {
839 // If we are skipping a non-empty bin, we notify the user
840 cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for " <<"protons Pt:"<< pt<< endl;
841 cout << " Bin content: " << Spectra[icharge*3+2]->GetBinContent(ibin) << endl;
842 }
843 }
844 }
845 else
846 {
847 Float_t FlukaCorrpPosTracking=fGFpPosTracking->Eval(Spectra[2]->GetBinCenter(ibin));
848 Float_t FlukaCorrpNegTracking=fGFpNegTracking->Eval(Spectra[2]->GetBinCenter(ibin));
849 // Printf("TPC/TOF Geant/Fluka Tracking: pt:%f Pos:%f Neg:%f",Spectra[2]->GetBinCenter(ibin),FlukaCorrpPosTracking,FlukaCorrpNegTracking);
850 Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*FlukaCorrpPosTracking);
851 Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*FlukaCorrpNegTracking);
852 Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError(ibin)*FlukaCorrpPosTracking);
853 Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError(ibin)*FlukaCorrpNegTracking);
854 Float_t FlukaCorrpPosMatching=fGFpPosMatching->Eval(Spectra[2]->GetBinCenter(ibin));
855 Float_t FlukaCorrpNegMatching=fGFpNegMatching->Eval(Spectra[2]->GetBinCenter(ibin));
856 // Printf("TPC/TOF Geant/Fluka Matching: pt:%f Pos:%f Neg:%f",Spectra[2]->GetBinCenter(ibin),FlukaCorrpPosMatching,FlukaCorrpNegMatching);
857 Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*FlukaCorrpPosMatching);
858 Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*FlukaCorrpNegMatching);
859 Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError(ibin)*FlukaCorrpPosMatching);
860 Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError(ibin)*FlukaCorrpNegMatching);
861 }
862 }
863 fGeanFlukaK->Close();
864 delete fGeanFlukaK;
865}
866
867
868///////////
869TF1 *
870TrackingEff_geantflukaCorrection(Int_t ipart, Int_t icharge)
871{
872
873 if (ipart == 3 && icharge == kNegative) {
874 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.);
875 return f;
876 }
877 else if (ipart == 4 && icharge == kNegative) {
878 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.);
879 }
880 else
881 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionNull(x)", 0., 5.);
882
883 return f;
884}
885
886Double_t
887TrackingPtGeantFlukaCorrectionNull(Double_t pTmc)
888{
889 return 1.;
890}
891
892Double_t
893TrackingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
894{
895 return (1 - 0.129758 *TMath::Exp(-pTmc*0.679612));
896}
897
898Double_t
899TrackingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
900{
901 return TMath::Min((0.972865 + 0.0117093*pTmc), 1.);
902}
903///////////////////////////////////////////
904TF1 *
905TOFmatchMC_geantflukaCorrection(Int_t ipart, Int_t icharge)
906{
907
908 if (ipart == 3 && icharge == kNegative) {
909 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.);
910 return f;
911 }
912 else if (ipart == 4 && icharge == kNegative) {
913 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.);
914 }
915 else
916 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionNull(x)", 0., 5.);
917
918 return f;
919}
920
921
922Double_t
923MatchingPtGeantFlukaCorrectionNull(Double_t pTmc)
924{
925 return 1.;
926}
927
928Double_t
929MatchingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
930{
931 Float_t ptTPCoutP =pTmc*(1-6.81059e-01*TMath::Exp(-pTmc*4.20094));
932 return (TMath::Power(1 - 0.129758*TMath::Exp(-ptTPCoutP*0.679612),0.07162/0.03471));
933}
934
935Double_t
936MatchingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
937{
938 Float_t ptTPCoutK=pTmc*(1- 3.37297e-03/pTmc/pTmc - 3.26544e-03/pTmc);
939 return TMath::Min((TMath::Power(0.972865 + 0.0117093*ptTPCoutK,0.07162/0.03471)), 1.);
940}
941void MatchingTOFEff(TH1F** Spectra, TList* list=0x0)
942{
943 if(TOFMatchingScalling[0]<0.0&&TOFMatchingScalling[1]<0.0)
944 {
945 TH1F *hMatcEffPos_data=(TH1F*)tcutsdata->GetHistoNMatchedPos();
946 hMatcEffPos_data->Divide((TH1F*)tcutsdata->GetHistoNSelectedPos());
947 hMatcEffPos_data->SetTitle("Matching Eff Pos - data");
948 TH1F *hMatcEffNeg_data=(TH1F*)tcutsdata->GetHistoNMatchedNeg();
949 hMatcEffNeg_data->Divide((TH1F*)tcutsdata->GetHistoNSelectedNeg());
950 hMatcEffNeg_data->SetTitle("Matching Eff Neg - data");
951 TH1F *hMatcEffPos_mc=(TH1F*)tcutsmc->GetHistoNMatchedPos();
952 hMatcEffPos_mc->Divide((TH1F*)tcutsmc->GetHistoNSelectedPos());
953 hMatcEffPos_mc->SetTitle("Matching Eff Pos - mc");
954 TH1F *hMatcEffNeg_mc=(TH1F*)tcutsmc->GetHistoNMatchedNeg();
955 hMatcEffNeg_mc->Divide((TH1F*)tcutsmc->GetHistoNSelectedNeg());
956 hMatcEffNeg_mc->SetTitle("Matching Eff Neg - mc");
957
958
959 hMatcEffPos_data->Divide(hMatcEffPos_mc);
960 hMatcEffNeg_data->Divide(hMatcEffNeg_mc);
961 hMatcEffPos_data->SetName("MatchingTOFPos");
962 hMatcEffNeg_data->SetName("MatchingTOFNeg");
963
964
965 TF1 *pol0MatchPos_data=new TF1("pol0MatchPos_data","pol0",.6,5);
966 hMatcEffPos_data->Fit("pol0MatchPos_data","MNR");
967 TF1 *pol0MatchNeg_data=new TF1("pol0MatchNeg_data","pol0",.6,5);
968 hMatcEffNeg_data->Fit("pol0MatchNeg_data","MNR");
969
970 list->Add(hMatcEffPos_data);
971 list->Add(hMatcEffNeg_data);
972
973
974 TOFMatchingScalling[0]=pol0MatchPos_data->GetParameter(0);
975 TOFMatchingScalling[1]=pol0MatchNeg_data->GetParameter(0);
976 }
977 //Correction spectra for matching efficiency
978 //For the moment I'm using the inclusive correction
979 for(Int_t ipart=0;ipart<3;ipart++)
980 {
981
982 for(Int_t ibin=1;ibin<Spectra[ipart]->GetNbinsX();ibin++)
983 {
984 Float_t ptspectra=Spectra[ipart]->GetBinCenter(ibin);
985 if(ptspectra<tcutsdata->GetPtTOFMatching())
986 continue;
987 //Spectra[ipart]->SetBinContent(ibin,( Spectra[ipart]->GetBinContent(ibin)/hMatcEffPos_data->GetBinContent(hMatcEffPos_data->FindBin(ptspectra))));
988 //Spectra[ipart+3]->SetBinContent(ibin,( Spectra[ipart+3]->GetBinContent(ibin)/hMatcEffNeg_data->GetBinContent(hMatcEffNeg_data->FindBin(ptspectra))));
989 Spectra[ipart]->SetBinContent(ibin,( Spectra[ipart]->GetBinContent(ibin)/TOFMatchingScalling[0]));
990 Spectra[ipart+3]->SetBinContent(ibin,( Spectra[ipart+3]->GetBinContent(ibin)/TOFMatchingScalling[1]));
991 }
992 }
993}
994Double_t eta2y(Double_t pt, Double_t mass, Double_t eta)
995{
996 Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
997 return TMath::ASinH(pt / mt * TMath::SinH(eta));
998}
999
1000TH1* GetSumAllCh(TH1F** spectra, Double_t* mass )
1001{
1002 TH1F* allch=(((TH1F*))spectra[0]->Clone("allCh"));
1003 allch->Reset();
1004 for (int i=0;i<6;i++)
1005 {
1006 Double_t masstmp=mass[i%3];
1007 for (int j=1;j<spectra[i]->GetXaxis()->GetNbins();j++)
1008 {
1009 Float_t value=spectra[i]->GetBinContent(j);
1010 Float_t error=spectra[i]->GetBinError(j);
1011 if(value>0.0)
1012 {
1013 Float_t pt=spectra[i]->GetXaxis()->GetBinCenter(j);
1014 Float_t mt2=masstmp*masstmp+pt*pt;
1015 Float_t mt=TMath::Sqrt(mt2);
1016 Float_t maxy=eta2y(pt,masstmp,0.8);
1017 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))));
1018 conver=conver/1.6;
1019 cout<<maxy<<" "<<conver<<" "<<masstmp<<""<<spectra[i]->GetName()<<endl;
1020 Float_t bincontent=allch->GetBinContent(j);
1021 Float_t binerror=allch->GetBinError(j);
1022 bincontent+=conver*value;
1023 binerror=TMath::Sqrt(binerror*binerror+conver*conver*error*error);
1024 allch->SetBinContent(j,bincontent);
1025 allch->SetBinError(j,binerror);
1026 }
1027
1028 }
1029 }
1030 Divideby2pipt(allch);
1031 allch->SetTitle("N_{ch};p_{T} (GeV/c);1/(2 #pi p_{T})dN/p_{T} |#eta|<0.8");
1032 return allch;
1033}
1034
1035void Divideby2pipt(TH1* hist)
1036{
1037
1038 for (int i=1;i<hist->GetXaxis()->GetNbins();i++)
1039 {
1040 Float_t value=hist->GetBinContent(i);
1041 Float_t error=hist->GetBinError(i);
1042 Float_t pt=hist->GetXaxis()->GetBinCenter(i);
1043 hist->SetBinContent(i,value/(2.0*TMath::Pi()*pt));
1044 hist->SetBinError(i,error/(2.0*TMath::Pi()*pt));
1045
1046 }
1047}