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