small changes in analysing macros
[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
26 Bool_t OpenFile(TString dirname, TString outputname, Bool_t mcflag);
27 void 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();
262         TList* canvaslist=new TList();
263         Float_t vertexcorrection=QAPlotsBoth(managerdata,managermc,ecutsdata,ecutsmc,tcutsdata,tcutsmc,listqa,canvaslist);
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);
276         canvaslist->Write("outputcanvas",TObject::kSingleKey);
277
278         fout->Close();
279
280 }
281
282 Bool_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
347 TH1F* 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         
374 void 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 }
412 void 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
433 void 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
689 void 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 }
695 void 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
705 void 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 }
720 void 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
731 void 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 ///////////
869 TF1 *
870 TrackingEff_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
886 Double_t
887 TrackingPtGeantFlukaCorrectionNull(Double_t pTmc)
888 {
889   return 1.;
890 }
891
892 Double_t
893 TrackingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
894 {
895   return (1 - 0.129758 *TMath::Exp(-pTmc*0.679612));
896 }
897
898 Double_t
899 TrackingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
900 {
901   return TMath::Min((0.972865 + 0.0117093*pTmc), 1.);
902 }
903 ///////////////////////////////////////////
904 TF1 *
905 TOFmatchMC_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
922 Double_t
923 MatchingPtGeantFlukaCorrectionNull(Double_t pTmc)
924 {
925   return 1.;
926 }
927
928 Double_t 
929 MatchingPtGeantFlukaCorrectionPrMinus(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
935 Double_t
936 MatchingPtGeantFlukaCorrectionKaMinus(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 }
941 void 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 }
994 Double_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
1000 TH1* 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
1035 void 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 }