Update of analysis macro
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / AnalysisBoth.C
1
2 class AliSpectraBothHistoManager;
3 class AliSpectraBothEventCuts; 
4 class AliSpectraBothTrackCuts;
5 TString Charge[]={"Pos","Neg"};
6 TString Sign[]={"Plus","Minus"};
7 TString Particle[]={"Pion","Kaon","Proton"};
8 AliSpectraBothHistoManager* managerdata=0x0;
9 AliSpectraBothEventCuts* ecutsdata=0x0; 
10 AliSpectraBothTrackCuts* tcutsdata=0x0;
11         
12 AliSpectraBothHistoManager* managermc=0x0;
13 AliSpectraBothEventCuts* ecutsmc=0x0; 
14 AliSpectraBothTrackCuts* tcutsmc=0x0;
15
16 Float_t TOFMatchingScalling[2]={-1,-1};
17 Int_t Color[3]={1,2,4};
18 Int_t Marker[6]={20,21,22,24,25,26};
19 Double_t Range[3]={0.3,0.3,0.5}; // LowPt range for pi k p
20
21
22 Double_t FitRange[2]={-2.0,2.0};
23 Double_t CutRange[2]={-3.0,3.0};
24 Double_t minptformaterial[6]={0.0,0.2,0.0,0.0,0.2,0.0};
25 Double_t maxptformaterial[6]={0.0,0.6,1.3,0.0,0.6,0.0};
26 Double_t minptforWD[6]={0.2,100.0,0.3,0.2,100.0,0.3};
27 Double_t maxptforWD[6]={1.5,-100.0,2.0,1.5,-100.0,2.0};
28 Double_t minRanges[3]={0.3,0.3,0.45};
29 Double_t maxRanges[3]={1.5,1.2,2.2};
30 Double_t fMaxContaminationPIDMC=0.2;
31
32
33 enum ECharge_t {
34   kPositive,
35   kNegative,
36   kNCharges
37 };
38 enum {
39  kdodca=0x1, //dca fits are made 
40  kgeantflukaKaon=0x2,// geant fluka correction is used for kaons 
41  kgeantflukaProton=0x4, // geant fluka correction is used for protons and antiprotons
42  knormalizationtoeventspassingPhySel=0x8,// spectra are divided by number of events passing physic selection   
43  kveretxcorrectionandbadchunkscorr=0x10, // correction for difference in z vertex distribution in data and MC and correction for bad chunks is applied
44  kmcisusedasdata=0x20, // the result of the looping over MC is used as data input 
45  kdonotusedcacuts=0x40, // allows to use the constant dca cut for all pt bins not the pt dependet defined in stardrad track cuts 2011
46  kuseprimaryPIDcont=0x80, //pid contamination is calculated using only primiary particle in this case K should use dca fits 
47  knormalizationwithbin0integralsdata=0x100, // the normalization factor is calcualte using integral over z vertex distributions (in this case reconstructed vertex disitrbution uses z vertex for data) 
48  knormalizationwithbin0integralsMC=0x200, //in this case reconstructed vertex disitrbution uses z vertex for data, those to options will be use only if knormalizationtoeventspassingPhySel is not set
49  kuserangeonfigfile=0x400, // use of config file for dca fit settings
50  kskipconcutonspectra=0x800, //do not use conPID<02 cut  useful for syst. studies
51  kuseTOFmatchingcorrection=0x1000 // if set tof matching correction is applied. 
52                                                         
53 };      
54
55 Bool_t OpenFile(TString dirname, TString outputname, Bool_t mcflag,Bool_t mcasdata=false);
56 void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TString outnamemc="",TString configfile="" )
57 {
58         gStyle->SetOptStat(0);  
59         TH1::AddDirectory(kFALSE);
60         gSystem->Load("libCore.so");  
61         gSystem->Load("libPhysics.so");
62         gSystem->Load("libTree");
63         gSystem->Load("libMatrix");
64         gSystem->Load("libSTEERBase");
65         gSystem->Load("libESD");
66         gSystem->Load("libAOD");
67         gSystem->Load("libANALYSIS");
68         gSystem->Load("libOADB");
69         gSystem->Load("libANALYSISalice");
70         gSystem->Load("libTENDER");
71         gSystem->Load("libCORRFW");
72         gSystem->Load("libPWGTools");
73         gSystem->Load("libPWGLFspectra");
74         
75         gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/QAPlotsBoth.C");
76         Double_t mass[3];
77         mass[0]   = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
78         mass[1]   = TDatabasePDG::Instance()->GetParticle("K+")->Mass();
79         mass[2] = TDatabasePDG::Instance()->GetParticle("proton")->Mass();
80
81         TFormula* dcacutxy=0x0;
82         if(!(options&kdonotusedcacuts))
83         {
84         
85                 AliESDtrackCuts* esdtrackcuts= AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE);
86                 TString formulastring(esdtrackcuts->GetMaxDCAToVertexXYPtDep());        
87                 formulastring.ReplaceAll("pt","x");
88                 dcacutxy=new TFormula("dcacutxy",formulastring.Data());
89         }
90         if(options&kuserangeonfigfile)
91                 if(!ReadConfigFile(configfile))
92                         return;         
93         TList* lout=new TList();
94         
95
96         TString indirname=Form("/output/train%s",outdate.Data());
97         //TString indirname("/output/train24072012");
98         if(outnamemc.Length()==0)
99         outnamemc=outnamedata;
100         cout<<indirname.Data()<<" "<<outnamemc.Data()<<endl;
101         // Do the job 
102
103
104         OpenFile(indirname,outnamemc,true);
105         OpenFile(indirname,outnamedata,false,((Bool_t)(options&kmcisusedasdata)));
106         if(!managermc||!managerdata)
107         {
108                 cout<<managermc<<" "<<managerdata<<endl;
109                 return; 
110         }
111         TH1F* rawspectradata[6];
112         TH1F* rawspectramc[6];
113         TH1F* MCTruth[6];
114         TH1F* eff[6];
115         TH1F* contallMC[6];
116         TH1F* contPID[6];
117         TH1F* contWD[6];
118         TH1F* contMat[6];
119         TH1F* confinal[6];
120         TH1F* contPIDpri[6];
121         TH1F* contSecMC[6];
122         
123         TH1F* contfit[12];
124         TH1F* contWDfit[12];
125         TH1F* contMatfit[12];
126         TH1F* primaryfit[12];
127
128         
129         
130         TH1F* spectra[6];
131         TH1F* spectraLeonardo[6];
132         
133         TH1F* corrLeonardo[6]; 
134         //GetSpectra(managerdata,rawspectradata,true);
135         //GetSpectra(managermc,rawspectramc,true,true);
136         
137         GetPtHistFromPtDCAhisto("hHistPtRecSigma","SpectraMC",managermc,rawspectramc,dcacutxy);
138         GetPtHistFromPtDCAhisto("hHistPtRecSigma","SpectraDATA",managerdata,rawspectradata,dcacutxy);
139         GetPtHistFromPtDCAhisto("hHistPtRecTruePrimary","eff",managermc,eff,dcacutxy);
140         GetPtHistFromPtDCAhisto("hHistPtRecTrue","conPID",managermc,contPID,dcacutxy);
141         GetPtHistFromPtDCAhisto("hHistPtRecSigmaSecondaryWeakDecay","conWD",managermc,contWD,dcacutxy);
142         GetPtHistFromPtDCAhisto("hHistPtRecSigmaSecondaryMaterial","conMat",managermc,contMat,dcacutxy);
143         GetPtHistFromPtDCAhisto("hHistPtRecSigmaPrimary","conPIDprimary",managermc,contPIDpri,dcacutxy);
144
145         
146         Double_t neventsmcall = 1 ;  //if loop over MC is done after or befor events cuts this will be changed 
147         Double_t neventsdata =  1;
148         Double_t neventsmc =  1;
149
150         //Normaliztion of MCtruth depends if the loop was done after of before ESD event cuts.
151         //In currect code this cannot be check on the level of macro.
152         //If the loop was done before MC should be done to all processed events (NumberOfProcessedEvents())
153         //If loop was done after MC should be normalized to all accepted events (NumberOfEvents()) 
154         // The option one will be alaways use.
155         
156         neventsmcall= ecutsmc->NumberOfProcessedEvents();
157         if(options&knormalizationtoeventspassingPhySel)
158         {
159                 //neventsmcall= ecutsmc->NumberOfProcessedEvents();
160                  neventsdata=ecutsdata->NumberOfPhysSelEvents();
161                  neventsmc=ecutsmc->NumberOfPhysSelEvents();
162         }
163         else if ((options&knormalizationwithbin0integralsdata)||(options&knormalizationwithbin0integralsMC))
164         {
165                 neventsdata=Normaliztionwithbin0integrals(options);
166                 neventsmc=ecutsmc->NumberOfPhysSelEvents();
167         }
168         else
169         {
170                 neventsdata=ecutsdata->NumberOfEvents(); //number of accepted events
171                  neventsmc=ecutsmc->NumberOfEvents();
172                 neventsmcall= ecutsmc->NumberOfEvents();
173
174         }
175         GetMCTruth(MCTruth);
176         cout<<neventsdata<<" Events"<<endl;
177         
178         
179         TH1F* allgen=((TH1F*)managermc->GetPtHistogram1D("hHistPtGen",1,1))->Clone();
180         allgen->SetName("AllGen");
181         TH1F* allrecMC=GetOneHistFromPtDCAhisto("hHistPtRec","rawallMC",managermc,dcacutxy);
182         TH1F* alleff=GetOneHistFromPtDCAhisto("hHistPtRecPrimary","effall",managermc,dcacutxy);
183         TH1F* allrecdata=GetOneHistFromPtDCAhisto("hHistPtRec","rawalldata",managerdata,dcacutxy);
184         
185         
186         
187         
188         TH1F* spectraall=(TH1F*)allrecdata->Clone("recNch");
189         spectraall->SetTitle("recNch");
190         TH1F* contall=(TH1F*)allrecMC->Clone("contall");
191         contall->SetTitle("contall");
192         //contall->Add(alleff,-1);
193         SubHistWithFullCorr(contall,alleff);
194         alleff->Divide(alleff,allgen,1,1,"B");
195         contall->Divide(contall,allrecMC,1,1,"B");
196         
197         GetCorrectedSpectra(spectraall,allrecdata,alleff,contall);
198         Divideby2pipt(spectraall);
199
200         allrecdata->Scale(1./neventsdata,"width");
201         allgen->Scale(1./neventsmcall,"width");
202         allrecMC->Scale(1./neventsmc,"width");
203         spectraall->Scale(1./neventsdata,"width");
204
205
206         lout->Add(allgen);
207         lout->Add(allrecMC);
208         lout->Add(alleff);
209         lout->Add(allrecdata);
210         lout->Add(spectraall);
211         lout->Add(contall);
212         
213         for (int i=0;i<6;i++)
214         {
215         
216                 
217                 TString tmpname(rawspectramc[i]->GetTitle());
218                 tmpname.ReplaceAll("SpectraMC","%s");
219                 contallMC[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contallMC"));
220                 contallMC[i]->SetTitle(Form(tmpname.Data(),"contallMC"));
221                 contfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contfit"));
222                 contfit[i]->SetTitle(Form(tmpname.Data(),"contfit"));
223                 contWDfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contWDfit"));
224                 contWDfit[i]->SetTitle(Form(tmpname.Data(),"contWDfit"));
225                 contMatfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contMatfit"));
226                 contMatfit[i]->SetTitle(Form(tmpname.Data(),"contMatfit"));
227                 primaryfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"primaryfit"));
228                 primaryfit[i]->SetTitle(Form(tmpname.Data(),"primaryfit"));
229                 contfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contfitonMC"));
230                 contfit[i+6]->SetTitle(Form(tmpname.Data(),"contfitonMC"));
231                 contWDfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contWDfitonMC"));
232                 contWDfit[i+6]->SetTitle(Form(tmpname.Data(),"contWDfitonMC"));
233                 contMatfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contMatfitonMC"));
234                 contMatfit[i+6]->SetTitle(Form(tmpname.Data(),"contMatfitonMC"));
235                 primaryfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"primaryfitMC"));
236                 primaryfit[i+6]->SetTitle(Form(tmpname.Data(),"primaryfitMC"));
237                 
238
239
240                 contfit[i]->Reset();
241                 contWDfit[i]->Reset();
242                 contMatfit[i]->Reset();
243                 primaryfit[i]->Reset();
244                 
245
246                 contfit[i+6]->Reset();
247                 contWDfit[i+6]->Reset();
248                 contMatfit[i+6]->Reset();
249                 primaryfit[i+6]->Reset();
250                 
251                 SetBintoOne(primaryfit[i]);
252                 SetBintoOne(primaryfit[i+6]);           
253                 spectra[i]=(TH1F*)rawspectradata[i]->Clone(Form(tmpname.Data(),"SpectraFinal"));
254                 spectra[i]->SetTitle(Form(tmpname.Data(),"SpectraFinal"));
255                 spectraLeonardo[i]=(TH1F*)rawspectradata[i]->Clone(Form(tmpname.Data(),"SpectraFinalLeonardo"));
256                 spectraLeonardo[i]->SetTitle(Form(tmpname.Data(),"SpectraFinalLeonardo"));
257
258                 corrLeonardo[i]=(TH1F*)MCTruth[i]->Clone(Form(tmpname.Data(),"CorrFactLeonardo"));
259                 corrLeonardo[i]->SetTitle(Form(tmpname.Data(),"CorrFactLeonardo"));                     
260                 corrLeonardo[i]->Divide(corrLeonardo[i],rawspectramc[i],1,1,"B");
261                 
262                 
263                 
264                 //contallMC[i]->Add(eff[i],-1.0);
265                 SubHistWithFullCorr(contallMC[i],eff[i]);
266                 //RecomputeErrors(contallMC[i]);
267                 contallMC[i]->Sumw2(); 
268                 contallMC[i]->Divide(contallMC[i],rawspectramc[i],1,1,"B");
269         
270                 // contamintaion from PID but only primaries
271                 //contPIDpri[i]->Add(eff[i],-1.0);
272                 SubHistWithFullCorr(contPIDpri[i],eff[i]);
273
274                 //RecomputeErrors(contPIDpri[i]);
275                 contPIDpri[i]->Divide(contPIDpri[i],rawspectramc[i],1,1,"B");
276         
277                 eff[i]->Divide(eff[i],MCTruth[i],1,1,"B");
278                 
279                 
280                 contPID[i]->Sumw2();
281                 rawspectramc[i]->Sumw2();
282                 //contPID[i]->Add(contPID[i],rawspectramc[i],-1,1);
283                 SubHistWithFullCorr(contPID[i],rawspectramc[i]);
284                 contPID[i]->Scale(-1.0);
285
286                 //RecomputeErrors(contPID[i]);
287                 contPID[i]->ResetStats();
288                 contPID[i]->Sumw2();
289                 contPID[i]->Divide(contPID[i],rawspectramc[i],1,1,"B");
290         
291                 if(options&kuseprimaryPIDcont)
292                         confinal[i]=(TH1F*)contPIDpri[i]->Clone(Form(tmpname.Data(),"confinal"));
293                 else    
294                         confinal[i]=(TH1F*)contPID[i]->Clone(Form(tmpname.Data(),"confinal"));
295                 confinal[i]->SetTitle(Form(tmpname.Data(),"confinal"));
296
297                 contSecMC[i]=(TH1F*)contWD[i]->Clone(Form(tmpname.Data(),"conSecMC"));
298                 contSecMC[i]->Add(contMat[i]);
299                 contWD[i]->Divide(contWD[i],rawspectramc[i],1,1,"B");
300                 contMat[i]->Divide(contMat[i],rawspectramc[i],1,1,"B");
301                 contSecMC[i]->Divide(contSecMC[i],rawspectramc[i],1,1,"B");
302         
303         
304                 rawspectradata[i]->Scale(1./neventsdata,"width");
305                 rawspectramc[i]->Scale(1./neventsmc,"width");
306                 MCTruth[i]->Scale(1./neventsmcall,"width");
307                 spectraLeonardo[i]->Scale(1./neventsdata,"width");
308         
309         
310         
311                 lout->Add(rawspectradata[i]);
312                 lout->Add(rawspectramc[i]);
313                 lout->Add(MCTruth[i]);
314                 lout->Add(eff[i]);
315                 lout->Add(contallMC[i]);
316                 lout->Add(contPID[i]);
317                 lout->Add(contWD[i]);
318                 lout->Add(contMat[i]);
319                 lout->Add(contfit[i]);
320                 lout->Add(contWDfit[i]);
321                 lout->Add(contMatfit[i]);
322                 lout->Add(primaryfit[i]);       
323                 lout->Add(contfit[i+6]);
324                 lout->Add(contWDfit[i+6]);
325                 lout->Add(contMatfit[i+6]);
326                 lout->Add(primaryfit[i+6]);     
327                 lout->Add(spectra[i]);
328                 lout->Add(spectraLeonardo[i]);
329                 lout->Add(confinal[i]);
330                 lout->Add(contPIDpri[i]);
331                 lout->Add(contSecMC[i]);
332         }
333         outdate.ReplaceAll("/","_");
334         configfile.ReplaceAll(".","_");
335         TFile* fout=0x0;
336         if(configfile.Length()>0&&(options&kuserangeonfigfile))
337                 fout=new TFile(Form("./results/ResMY_%s_%s_%#X_%s.root",outnamemc.Data(),outdate.Data(),options,configfile.Data()),"RECREATE");
338         else
339                 fout=new TFile(Form("./results/ResMY_%s_%s_%#X.root",outnamemc.Data(),outdate.Data(),options),"RECREATE");
340         if (options&kdodca)
341                 DCACorrectionMarek(managerdata,managermc,dcacutxy,fout,contfit,contWDfit,contMatfit,primaryfit);
342         for (int i=0;i<6;i++)
343         {
344                         if(options&kdodca)
345                         {
346                                 if((options&kuseprimaryPIDcont)||(i!=1&&i!=4)) //if we do not use cont PId only for primary  and this is a kaon that do not use fit
347                                         confinal[i]->Add(contfit[i]);
348                                 GetCorrectedSpectra(spectra[i],rawspectradata[i],eff[i],confinal[i]);
349                         }
350                         else
351                         {
352                                 GetCorrectedSpectra(spectra[i],rawspectradata[i],eff[i],contallMC[i]);  
353                         }
354                         GetCorrectedSpectraLeonardo(spectraLeonardo[i],corrLeonardo[i],primaryfit[i],primaryfit[i+6]);
355                         if(options&kskipconcutonspectra)
356                                 continue;
357                         if(options&kuseprimaryPIDcont)
358                         {
359                                 CleanHisto(spectra[i],-1,100,contPIDpri[i]);
360                                 CleanHisto(spectraLeonardo[i],-1,100,contPIDpri[i]);            
361                         }
362                         else
363                         {
364                                 CleanHisto(spectra[i],-1,100,contPID[i]);
365                                 CleanHisto(spectraLeonardo[i],-1,100,contPID[i]);
366                         }                               
367         }
368         
369         GFCorrection(spectra,tcutsdata->GetPtTOFMatching(),options);
370         GFCorrection(spectraLeonardo,tcutsdata->GetPtTOFMatching(),options);
371                 TH1F* allch=GetSumAllCh(spectra,mass);
372         lout->Add(allch);       
373         if(options&kuseTOFmatchingcorrection)
374         {       
375                 MatchingTOFEff(spectra,lout);
376                 MatchingTOFEff(spectraLeonardo);
377                 TOFMatchingForNch(spectraall);
378         }
379 //      lout->ls();
380         fout->cd();     
381         TList* listqa=new TList();
382         TList* canvaslist=new TList();
383         Float_t vertexcorrection=1.0;
384         Float_t corrbadchunksvtx=1.0;
385         if ((options&knormalizationwithbin0integralsdata)||(options&knormalizationwithbin0integralsMC))
386                 corrbadchunksvtx=QAPlotsBoth(managerdata,managermc,ecutsdata,ecutsmc,tcutsdata,tcutsmc,listqa,canvaslist,0);
387         else
388                 corrbadchunksvtx=QAPlotsBoth(managerdata,managermc,ecutsdata,ecutsmc,tcutsdata,tcutsmc,listqa,canvaslist,1);
389         if (options&kveretxcorrectionandbadchunkscorr)
390                 vertexcorrection=corrbadchunksvtx;
391         cout<<" VTX corr="<<vertexcorrection<<endl;
392         Double_t ycut=tcutsdata->GetY();
393         if(TMath::Abs(ycut)>0.0)
394                 vertexcorrection=vertexcorrection/(2.0*ycut);
395         for (int i=0;i<6;i++)
396         {
397                 spectra[i]->Scale(vertexcorrection);
398                 spectraLeonardo[i]->Scale(vertexcorrection);
399                 if(TMath::Abs(ycut)>0.0)
400                 {
401                         rawspectradata[i]->Scale(1.0/(2.0*ycut));
402                         rawspectramc[i]->Scale(1.0/(2.0*ycut));
403                         MCTruth[i]->Scale(1.0/(2.0*ycut));
404                 }
405         }       
406         allch->Scale(vertexcorrection);
407         spectraall->Scale(vertexcorrection/1.6);
408
409         //spectraall->Scale(1.0/1.6);
410         lout->Write("output",TObject::kSingleKey);      
411         listqa->Write("outputQA",TObject::kSingleKey);
412         canvaslist->Write("outputcanvas",TObject::kSingleKey);
413
414         fout->Close();
415         //Normaliztionwithbin0integrals();
416
417 }
418
419 Bool_t   OpenFile(TString dirname,TString outputname, Bool_t mcflag, Bool_t mcasdata)
420 {
421         
422
423         TString nameFile = Form("./%s/AnalysisResults%s.root",dirname.Data(),(mcflag?"MC":"DATA"));
424         TFile *file = TFile::Open(nameFile.Data());
425         if(!file)
426         {
427                 cout<<"no file"<<endl;
428                 return false;
429         }       
430         TString sname=Form("OutputBothSpectraTask_%s_%s",(mcflag?"MC":"Data"),outputname.Data());
431         if(mcasdata)
432         {
433                 cout<<"using MC as data "<<endl;
434                 sname=Form("OutputBothSpectraTask_%s_%s","MC",outputname.Data());
435         }
436         file->ls();
437         TDirectoryFile *dir=(TDirectoryFile*)file->Get(sname.Data());
438         if(!dir)
439         {
440         //      cout<<"no dir "<<sname.Data()<<endl;    
441                 if(mcasdata)
442                 {
443                         cout<<"using MC as data "<<endl;
444                         sname=Form("OutputAODSpectraTask_%s_%s","MC",outputname.Data());
445                 }
446                 else    
447                         sname=Form("OutputAODSpectraTask_%s_%s",(mcflag?"MC":"Data"),outputname.Data());
448         //      cout<<"trying "<<sname.Data()<<endl;
449                 dir=(TDirectoryFile*)file->Get(sname.Data());
450                 if(!dir)
451                 {
452                         cout<<"no dir "<<sname.Data()<<endl;
453                         return false;
454                 }
455         }
456         cout << " -- Info about " <<(mcflag?"MC":"DATA") <<" -- "<< endl;
457         if(mcflag)
458         {
459                 managermc= (AliSpectraBothHistoManager*) dir->Get("SpectraHistos");     
460                 ecutsmc = (AliSpectraBothEventCuts*) dir->Get("Event Cuts");
461                 tcutsmc = (AliSpectraBothTrackCuts*) dir->Get("Track Cuts");
462                 ecutsmc->PrintCuts();
463                 tcutsmc->PrintCuts();
464                 if(!managermc||!ecutsmc||!tcutsmc)
465                         return false;
466                 if(managermc->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()!=ecutsmc->GetHistoCuts()->GetBinContent(3))
467                         cout<<"Please check MC file possible problem with merging"<<" "<<managermc->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()<<" "<<ecutsmc->GetHistoCuts()->GetBinContent(3)<<endl;
468         }
469         else
470         {
471                 managerdata= (AliSpectraBothHistoManager*) dir->Get("SpectraHistos");   
472                 ecutsdata = (AliSpectraBothEventCuts*) dir->Get("Event Cuts");
473                 tcutsdata = (AliSpectraBothTrackCuts*) dir->Get("Track Cuts");
474                 ecutsdata->PrintCuts();
475                 tcutsdata->PrintCuts();
476                 if(!managerdata||!ecutsdata||!tcutsdata)
477                         return false;
478                 if(managerdata->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()!=ecutsdata->GetHistoCuts()->GetBinContent(3))
479                         cout<<"Please check DATA file possible problem with merging"<<" "<<anagerdata->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()<<" "<<ecutsdata->GetHistoCuts()->GetBinContent(3)<<endl;
480
481         }
482         return true;
483 }
484
485  void GetMCTruth(TH1F** MCTruth)
486  {
487         for(Int_t icharge=0;icharge<2;icharge++)
488         {
489                 for(Int_t ipart=0;ipart<3;ipart++)
490                 {
491                         Int_t index=ipart+3*icharge;
492                         TString hname=Form("hHistPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
493                         MCTruth[index]=(TH1F*)((TH1F*)managermc->GetPtHistogram1D(hname.Data(),1,1))->Clone();
494                         MCTruth[index]->SetName(Form("MCTruth_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
495                         MCTruth[index]->SetTitle(Form("MCTruth_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
496                         MCTruth[index]->Sumw2(); 
497                 }
498         }
499 }
500
501 TH1F* GetOneHistFromPtDCAhisto(TString name,TString hnameout,AliSpectraBothHistoManager* hman,TFormula* dcacutxy)
502 {
503                         histo =(TH1F*)((TH1F*) hman->GetPtHistogram1D(name.Data(),-1,-1))->Clone();
504                         histo->SetName(hnameout.Data());
505                         histo->SetTitle(hnameout.Data());
506                   
507                         if(dcacutxy)
508                         {
509                                 for(int ibin=1;ibin<histo->GetNbinsX();ibin++)
510                                 {
511                                         Double_t lowedge=histo->GetBinLowEdge(ibin);
512                                         Float_t cut=dcacutxy->Eval(lowedge);
513                                         TH1F* dcahist=(TH1F*)hman->GetDCAHistogram1D(name.Data(),lowedge,lowedge);
514 //                                      Float_t inyield=dcahist->Integral(dcahist->GetXaxis()->FindBin(-1.0*cut),dcahist->GetXaxis()->FindBin(cut));
515                                         Float_t testyield=0.0;
516                                         Float_t testerror=0.0;  
517                                         for (int itest=dcahist->GetXaxis()->FindBin(-1.0*cut);itest<=dcahist->GetXaxis()->FindBin(cut);itest++)
518                                         {
519                                                 testyield+=dcahist->GetBinContent(itest);
520                                                 testerror+=dcahist->GetBinError(itest)*dcahist->GetBinError(itest);
521                                         }
522 //                                      cout<<"corr data "<<histo->GetBinContent(ibin)<<" "<<inyield<<" "<<dcahist->Integral()<<" "<<hnameout.Data()<<endl;
523 //                                      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;
524
525 //                                      cout<<testyield<<" "<<TMath::Sqrt(testerror)<<" error2 "<<inyield<<" "<<TMath::Sqrt(inyield)<<endl;
526                                                 
527                                         //histo->SetBinContent(ibin,inyield);
528                                         //histo->SetBinError(ibin,TMath::Sqrt(inyield));
529                                         histo->SetBinContent(ibin,testyield);
530                                         histo->SetBinError(ibin,TMath::Sqrt(testerror));
531
532                                 }
533                         }
534                         histo->Sumw2();
535                         return histo;
536 }
537
538
539
540         
541 void GetPtHistFromPtDCAhisto(TString hnamein, TString hnameout, AliSpectraBothHistoManager* hman,TH1F** histo,TFormula* dcacutxy)
542 {
543         //Float_t min[3]={0.3,0.3,0.45};
544         //Float_t max[3]={1.5,1.2,2.2};
545         for(Int_t icharge=0;icharge<2;icharge++)
546         {
547                 for(Int_t ipart=0;ipart<3;ipart++)
548                 {
549                         Int_t index=ipart+3*icharge;
550                         Printf("Getting %s",hnamein.Data());
551                         TString nameinfinal=Form("%s%s%s",hnamein.Data(),Particle[ipart].Data(),Sign[icharge].Data());
552                         TString nameoutfinal=Form("%s%s%s",hnameout.Data(),Particle[ipart].Data(),Sign[icharge].Data());
553                         
554                         
555                         histo[index]=GetOneHistFromPtDCAhisto(nameinfinal,nameoutfinal,hman,dcacutxy);
556                         CleanHisto(histo[index],minRanges[ipart],maxRanges[ipart]);
557                 }
558         } 
559 }
560 void CleanHisto(TH1F* h, Float_t minV, Float_t maxV,TH1* contpid=0x0)
561 {
562         for (int i=0;i<=h->GetNbinsX();i++)
563         {       
564                 if(h->GetXaxis()->GetBinCenter(i)<minV||h->GetXaxis()->GetBinCenter(i)>maxV)
565                 {
566                         h->SetBinContent(i,0);
567                         h->SetBinError(i,0);
568                 }       
569                 if(contpid)
570                 {
571                         if(contpid->GetBinContent(i)>fMaxContaminationPIDMC)
572                         {
573                                 h->SetBinContent(i,0);
574                                 h->SetBinError(i,0);
575                         }
576                 }
577         }
578 }
579
580
581 void DCACorrectionMarek(AliSpectraBothHistoManager* hman_data, AliSpectraBothHistoManager* hman_mc,TFormula* fun,TFile *fout,TH1F** hcon,TH1F** hconWD,TH1F** hconMat,TH1F** hprimary)
582 {
583   printf("\n\n-> DCA Correction");  
584   
585  
586   Printf("\DCACorr");
587   TString sample[2]={"data","mc"};
588   ofstream debug("debugDCA.txt");
589   TList* listofdcafits=new TList();
590   for(Int_t icharge=0;icharge<2;icharge++)
591   {
592                 for(Int_t ipart=0;ipart<3;ipart++)
593                 {
594                         Int_t index=ipart+3*icharge;
595                         for(Int_t isample=0;isample<2;isample++)
596                         {
597
598                                 for(Int_t ibin_data=7;ibin_data<40;ibin_data++)
599                                 {       
600                                                 
601                                         Double_t lowedge=hcon[index]->GetBinLowEdge(ibin_data);
602                                         Double_t binwidth=hcon[index]->GetBinWidth(ibin_data);
603                                         debug<<"NEW "<<Particle[ipart].Data()<<" "<<Sign[icharge].Data()<<" "<<lowedge<<endl;
604                                         if(fun)
605                                         {                                               
606                                                 CutRange[1]=fun->Eval(lowedge);
607                                                 CutRange[0]=-1.0*CutRange[1];
608                                         }       
609                                         debug<<"cut  "<<CutRange[1]<<" "<<CutRange[0]<<endl;            
610                                         Short_t fitsettings=DCAfitsettings(lowedge+0.5*binwidth,index);
611                                         debug<<"settings "<< fitsettings<<endl;
612                                         if(fitsettings==0)
613                                                 continue;       
614                                                 
615                                         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);
616                                         TLegend* Leg1 = new TLegend(0.6,0.6,0.85,0.85,"","NDC");
617                                         Leg1->SetFillStyle(kFALSE);
618                                         Leg1->SetLineColor(kWhite);
619                                         Leg1->SetBorderSize(0);
620
621                                         
622                                         if(isample==0)
623                                                 TH1F *hToFit =(TH1F*) ((TH1F*)hman_data->GetDCAHistogram1D(Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
624                                         if(isample==1)
625                                                 TH1F *hToFit =(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
626                                         debug<<Particle[ipart].Data()<<" "<<Sign[icharge].Data()<<" "<<lowedge<<endl;
627                                         TH1F *hmc1=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
628                                         TH1F *hmc2=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaSecondaryWeakDecay%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
629                                         TH1F *hmc3=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaSecondaryMaterial%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
630                                         Double_t minentries=1;
631                                         debug<<hToFit->GetEntries()<<" "<<hmc1->GetEntries()<<" "<<hmc2->GetEntries()<<" "<<hmc3->GetEntries()<<endl;
632                                         debug<<((fitsettings&0x1)&&hmc2->GetEntries()<=minentries)<<" "<<((fitsettings&0x2)&&hmc3->GetEntries()<=minentries)<<endl;
633                                         if(hToFit->GetEntries()<=minentries || hmc1->GetEntries()<=minentries || ((fitsettings&0x1)&&hmc2->GetEntries()<=minentries) || ((fitsettings&0x2)&&hmc3->GetEntries()<=minentries))
634                                                 continue;
635                                         hToFit->Sumw2();
636                                         hmc1->Sumw2();
637                                         hmc2->Sumw2();
638                                         hmc3->Sumw2();
639                                         
640                                         Float_t corrforrebinning[4]={1.0,1.0,1.0,1.0};  
641                                         Int_t binCutRange[]={hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1])};                    
642
643                                         
644                                         if(hmc3->GetNbinsX()>300)
645                                         {
646                                         
647                                                 corrforrebinning[0]=hToFit->Integral(binCutRange[0],binCutRange[1]);
648                                                 corrforrebinning[1]=hmc1->Integral(binCutRange[0],binCutRange[1]);
649                                                 corrforrebinning[2]=hmc2->Integral(binCutRange[0],binCutRange[1]);
650                                                 corrforrebinning[3]=hmc3->Integral(binCutRange[0],binCutRange[1]);
651
652                                                 hToFit->Rebin(30);
653                                                 hmc1->Rebin(30);
654                                                 hmc2->Rebin(30);
655                                                 hmc3->Rebin(30);
656
657                                                 binCutRange[0]=hmc1->GetXaxis()->FindBin(CutRange[0]);
658                                                 binCutRange[1]=hmc1->GetXaxis()->FindBin(CutRange[1]);
659
660
661                                                 //after rebbing we lose resolution of the dca this correction also us to do obtain inside used dca
662
663                                                 if(hToFit->Integral(binCutRange[0],binCutRange[1])>0.0)
664                                                         corrforrebinning[0]=corrforrebinning[0]/hToFit->Integral(binCutRange[0],binCutRange[1]);
665                                                 else    
666                                                         corrforrebinning[0]=1.0;
667                                                 if(hmc1->Integral(binCutRange[0],binCutRange[1])>0.0)
668                                                         corrforrebinning[1]=corrforrebinning[1]/hmc1->Integral(binCutRange[0],binCutRange[1]);
669                                                 else    
670                                                         corrforrebinning[1]=1.0;
671                                                 if(hmc2->Integral(binCutRange[0],binCutRange[1])>0.0)
672                                                         corrforrebinning[2]=corrforrebinning[2]/hmc2->Integral(binCutRange[0],binCutRange[1]);
673                                                 else    
674                                                         corrforrebinning[2]=1.0;
675                                                 if(hmc3->Integral(binCutRange[0],binCutRange[1])>0.0)
676                                                         corrforrebinning[3]=corrforrebinning[3]/hmc3->Integral(binCutRange[0],binCutRange[1]);
677                                                 else    
678                                                         corrforrebinning[3]=1.0;
679                                                         
680
681
682                                         }
683
684                                         cDCA->cd();
685                                         gPad->SetGridy();
686                                         gPad->SetGridx();
687                                         gPad->SetLogy();
688          
689                                         TObjArray *mc=0x0;
690                                         Int_t Npar=3;   
691                                         if(fitsettings==3)
692                                                 mc = new TObjArray(3);        // MC histograms are put in this array
693                                         else
694                                         {
695                                                 mc = new TObjArray(2);
696                                                 Npar=2;
697                                         }
698                                                 
699                                         mc->Add(hmc1);
700                                         if(fitsettings&0x1)
701                                                 mc->Add(hmc2);
702                                         if(fitsettings&0x2)
703                                                 mc->Add(hmc3);
704                                         TFractionFitter* fit = new TFractionFitter(hToFit,mc); // initialise
705                                         fit->Constrain(0,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
706                                         Double_t defaultStep = 0.01;
707                                         for (int iparm = 0; iparm < Npar; ++iparm) 
708                                         {
709                                                 TString name("frac"); name += iparm;
710                                                 if(iparm==0)
711                                                         fit->GetFitter()->SetParameter(iparm,name.Data(),0.85,0.01,0.0,1.0);
712                                                 else if (Npar==2)
713                                                         fit->GetFitter()->SetParameter(iparm,name.Data(),0.15,0.01,0.0,1.0);
714                                                 else
715                                                         fit->GetFitter()->SetParameter(iparm,name.Data(),0.075,0.01,0.0,1.0);
716         
717                                         }
718                                         if(fitsettings&0x1)
719                                                 fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
720                                         if(fitsettings&0x2)
721                                                 fit->Constrain(1+(fitsettings&0x1),0.0,1.0);               // constrain fraction 1 to be between 0 and 1
722                                         
723                                         Int_t binFitRange[]={hmc1->GetXaxis()->FindBin(FitRange[0]),hmc1->GetXaxis()->FindBin(FitRange[1])};                    
724                                         fit->SetRangeX(binFitRange[0],binFitRange[1]);
725                                         hToFit->GetXaxis()->SetRange(binFitRange[0],binFitRange[1]);
726                                         hToFit->SetTitle(Form("DCA distr - %s %s %s %lf",sample[isample].Data(),Particle[ipart].Data(),Sign[icharge].Data(),lowedge));
727                                         Int_t status = fit->Fit();               // perform the fit
728                                         cout << "fit status: " << status << endl;
729                                         debug<<"fit status: " << status << endl;
730                 
731                                         if (status == 0)
732                                         {       
733                                                 Double_t v1=0.0,v2=0.0,v3=0.0;
734                                                 Double_t ev1=0.0,ev2=0.0,ev3=0.0;
735                                                 Double_t cov=0.0;
736                                                                    
737                                                 // check on fit status
738                                                 TH1F* result = (TH1F*) fit->GetPlot();
739                                                 TH1F* PrimMCPred=(TH1F*)fit->GetMCPrediction(0);
740                                                 
741                                                 TH1F* secStMCPred=0X0;
742                                         
743                                                 TH1F* secMCPred=0x0;
744             
745                                                 //first method, use directly the fit result
746                                                 fit->GetResult(0,v1,ev1);
747
748                                                 if(fitsettings&0x1)
749                                                 {
750                                                         fit->GetResult(1,v2,ev2);
751                                                         secStMCPred=(TH1F*)fit->GetMCPrediction(1);
752
753                                                 }
754                                                 if(fitsettings&0x2)
755                                                 {
756                                                         fit->GetResult(1+(fitsettings&0x1),v3,ev3);
757                                                         secMCPred=(TH1F*)fit->GetMCPrediction(1+(fitsettings&0x1));
758                                                         if(fitsettings&0x1)
759                                                                 cov=fit->GetFitter()->GetCovarianceMatrixElement(1,2);
760
761
762                                                 }
763                                                 
764                                         
765                                                 debug<<v1<<" "<<ev1<<" "<<v2<<" "<<ev2<<" "<<v3<<" "<<ev3<<" "<<" "<<cov<<endl;
766                                         
767                                                 // becuase dca cut range is not a fit range the results from TFractionFitter should be rescale
768                                                 // This can be done in two ways or use input histograms or output histograms
769                                                 // The difference between those two methods should be on the level of statistical error         
770                                                 // I use output histograms 
771                                                 
772                                                 // Method 1 input histo         
773
774                                                 Float_t normalizationdata=hToFit->Integral(hToFit->GetXaxis()->FindBin(CutRange[0]),hToFit->GetXaxis()->FindBin(CutRange[1]))/hToFit->Integral(binFitRange[0],binFitRange[1]);
775                                                 normalizationdata*=corrforrebinning[0];
776
777                                                 Float_t normalizationmc1=(hmc1->Integral(binCutRange[0],binCutRange[1])/hmc1->Integral(binFitRange[0],binFitRange[1]))/normalizationdata;
778                                                 Float_t normalizationmc2=0.0;
779                                                 if(fitsettings&0x1)
780                                                         normalizationmc2=(hmc2->Integral(binCutRange[0],binCutRange[1])/hmc2->Integral(binFitRange[0],binFitRange[1]))/normalizationdata;
781                                                 Float_t normalizationmc3=0.0;
782                                                 if(fitsettings&0x2)
783                                                         normalizationmc3=(hmc3->Integral(binCutRange[0],binCutRange[1])/hmc3->Integral(binFitRange[0],binFitRange[1]))/normalizationdata;
784
785                                                 normalizationmc1*=corrforrebinning[1];
786                                                 normalizationmc2*=corrforrebinning[2];
787                                                 normalizationmc3*=corrforrebinning[3];
788
789                                                 debug<<"After Nor"<<endl;
790                                                 debug<<v1*normalizationmc1<<" "<<ev1*normalizationmc1<<" "<<v2*normalizationmc2<<" "<<ev2*normalizationmc2<<" "<<v3*normalizationmc3<<" "<<ev3*normalizationmc3<<" "<<endl;
791                                                 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;
792                                                 debug<<"addtional info"<<endl;
793
794
795                                                //Method 2 output histo          
796
797                                                 Float_t normalizationdata1=result->Integral(binCutRange[0],binCutRange[1])/result->Integral(binFitRange[0],binFitRange[1]);
798                                                 
799
800                                                 // if the cut range is bigger the fit range we should calculate the normalization factor for data using the data histogram 
801                                                 // because result histogram has entries only in fits range       
802                                                 if(FitRange[0]>CutRange[0]||FitRange[1]<CutRange[1])    
803                                                         normalizationdata1=normalizationdata;
804                                         
805                                                 normalizationdata1*=corrforrebinning[0];
806
807
808                                                 Float_t normalizationmc11=(PrimMCPred->Integral(binCutRange[0],binCutRange[1])/PrimMCPred->Integral(binFitRange[0],binFitRange[1]))/normalizationdata1;
809                                                 Float_t normalizationmc21=0.0;
810                                                 if(fitsettings&0x1)
811                                                         normalizationmc21=(secStMCPred->Integral(binCutRange[0],binCutRange[1])/secStMCPred->Integral(binFitRange[0],binFitRange[1]))/normalizationdata1;
812                                                 Float_t normalizationmc31=0.0;
813                                                 if(fitsettings&0x2)
814                                                         normalizationmc31=(secMCPred->Integral(binCutRange[0],binCutRange[1])/secMCPred->Integral(binFitRange[0],binFitRange[1]))/normalizationdata1;
815                                                 
816                                                 normalizationmc11*=corrforrebinning[1];
817                                                 normalizationmc21*=corrforrebinning[2];
818                                                 normalizationmc31*=corrforrebinning[3];
819
820                                                 debug<<"After Nor 2"<<endl;
821                                                 debug<<v1*normalizationmc11<<" "<<ev1*normalizationmc11<<" "<<v2*normalizationmc21<<" "<<ev2*normalizationmc21<<" "<<v3*normalizationmc31<<" "<<ev3*normalizationmc31<<endl;
822                                                 
823                                                 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;
824                                         
825                                                 
826                                                 hconWD[index+6*isample]->SetBinContent(ibin_data,v2*normalizationmc21);
827                                                 hconWD[index+6*isample]->SetBinError(ibin_data,ev2*normalizationmc21);
828                                                 hconMat[index+6*isample]->SetBinContent(ibin_data,v3*normalizationmc31);
829                                                 hconMat[index+6*isample]->SetBinError(ibin_data,ev3*normalizationmc31);
830                                                 hprimary[index+6*isample]->SetBinContent(ibin_data,v1*normalizationmc11);
831                                                 hprimary[index+6*isample]->SetBinError(ibin_data,ev1*normalizationmc11);
832                                                 hcon[index+6*isample]->SetBinContent(ibin_data,v2*normalizationmc21+v3*normalizationmc31);
833                                                 hcon[index+6*isample]->SetBinError(ibin_data,TMath::Sqrt(ev2*ev2*normalizationmc21*normalizationmc21+ev3*ev3*normalizationmc31*normalizationmc31+cov*normalizationmc31*normalizationmc21));
834                                                 
835                                                 
836                                                 
837                                                 //Drawing section
838                                                 result->Scale(1.0/result->Integral(result->GetXaxis()->FindBin(FitRange[0]),result->GetXaxis()->FindBin(FitRange[1])));
839                                                 hToFit->Scale(1.0/hToFit->Integral(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1])));
840                                                 PrimMCPred->Scale(v1/PrimMCPred->Integral(PrimMCPred->GetXaxis()->FindBin(FitRange[0]),PrimMCPred->GetXaxis()->FindBin(FitRange[1])));
841
842                                                 hToFit->SetMinimum(0.0001);
843                                                 hToFit->DrawClone("E1x0");
844                                                 result->SetTitle("Fit result");
845                                                 result->SetLineColor(kBlack);
846                                                 Leg1->AddEntry(result,"result","lp");
847                                                 result->DrawClone("lhistsame");
848                                         
849                                                 PrimMCPred->SetLineColor(kGreen+2);
850                                                 PrimMCPred->SetLineStyle(2);
851                                                  PrimMCPred->SetLineWidth(3.0);
852                                                 Leg1->AddEntry(PrimMCPred,"Prmi.","l");
853                                                 PrimMCPred->DrawClone("lhistsame");
854                                                 if(fitsettings&0x1)
855                                                 {
856
857                                                         secStMCPred->Scale(v2/secStMCPred->Integral(secStMCPred->GetXaxis()->FindBin(FitRange[0]),secStMCPred->GetXaxis()->FindBin(FitRange[1])));
858                                                         secStMCPred->SetLineColor(kRed);
859                                                         secStMCPred->SetLineWidth(3.0);
860
861                                                         secStMCPred->SetLineStyle(3);
862                                                         Leg1->AddEntry(secStMCPred,"Sec.WD","l");
863                                                         secStMCPred->DrawClone("lhistsame");
864
865                                                 }
866                                                 if(fitsettings&0x2)
867                                                 {
868                                                         
869                                                         secMCPred->Scale(v3/secMCPred->Integral(secMCPred->GetXaxis()->FindBin(FitRange[0]),secMCPred->GetXaxis()->FindBin(FitRange[1])));
870                                                         secMCPred->SetLineColor(kBlue);
871                                                         secMCPred->SetLineWidth(3.0);
872
873                                                         secMCPred->SetLineStyle(4);     
874                                                         Leg1->AddEntry(secMCPred,"Sec.Mat","l");
875                                                         secMCPred->DrawClone("lhistsame");
876             
877                                                 }   
878                                         }                               
879                                         else
880                                         {
881                                                 hconWD[index+6*isample]->SetBinContent(ibin_data,0.0);
882                                                 hconWD[index+6*isample]->SetBinError(ibin_data,0.0);
883                                                 hconMat[index+6*isample]->SetBinContent(ibin_data,0.0);
884                                                 hconMat[index+6*isample]->SetBinError(ibin_data,0.0);                                   
885                                                 hcon[index+6*isample]->SetBinContent(ibin_data,0.0);
886                                                 hcon[index+6*isample]->SetBinError(ibin_data,0.0);
887                                                 hprimary[index+6*isample]->SetBinContent(ibin_data,1.0);
888                                                 hprimary[index+6*isample]->SetBinError(ibin_data,0.0);
889                                         }
890                                         Leg1->Draw();
891                                         listofdcafits->Add(cDCA);
892                                         
893                                         //cDCA->Write();
894                                         delete hToFit;
895                                 }
896         
897
898                         }
899
900                 }
901         }
902         fout->cd();
903         listofdcafits->Write("DCAfits",TObject::kSingleKey);    
904 }
905
906 void RecomputeErrors(TH1* h)
907 {
908         for (int i=0; i<=h->GetXaxis()->GetNbins(); i++)
909         {
910                 cout<<h->GetBinContent(i)<<" "<<h->GetBinError(i)<<" error "<<TMath::Sqrt(h->GetBinContent(i))<<endl;
911                 h->SetBinError(i,TMath::Sqrt(h->GetBinContent(i)));
912         }
913         h->Sumw2();     
914 }
915 void SetBintoOne(TH1* h)
916 {
917         for (int i=0;i<=h->GetXaxis()->GetNbins(); i++) 
918         {
919                 h->SetBinContent(i,1);
920                 h->SetBinError(i,0);
921         }
922 }
923
924
925 void GetCorrectedSpectra(TH1F* corr,TH1F* raw,TH1F* eff, TH1F* con)
926 {
927         for (int i=0;i<=corr->GetXaxis()->GetNbins(); i++) 
928         {
929                 corr->SetBinContent(i,1);
930                 corr->SetBinError(i,0);
931         }
932         corr->Sumw2(); 
933         corr->Add(con,-1);
934         corr->Sumw2();  
935         corr->Divide(eff);
936         corr->Sumw2(); 
937         corr->Multiply(raw);
938         corr->Sumw2(); 
939 }
940 void GetCorrectedSpectraLeonardo(TH1F* spectra,TH1F* correction, TH1F* hprimaryData,TH1F* hprimaryMC)
941 {
942         spectra->Sumw2(); 
943         spectra->Multiply(correction);
944         spectra->Sumw2(); 
945         hprimaryData->Sumw2(); 
946         spectra->Multiply(hprimaryData);
947         hprimaryMC->Sumw2(); 
948         spectra->Divide(hprimaryMC);
949 }
950
951 void GFCorrection(TH1F **Spectra,Float_t tofpt,UInt_t options)
952 {
953         if (options&kgeantflukaKaon)
954         {               
955                  //Geant/Fluka Correction
956                 Printf("\nGF correction for Kaons");
957                 //Getting GF For Kaons in TPC
958                 TGraph *gGFCorrectionKaonPlus=new TGraph();
959                 gGFCorrectionKaonPlus->SetName("gGFCorrectionKaonPlus");
960                 gGFCorrectionKaonPlus->SetTitle("gGFCorrectionKaonPlus");
961                 TGraph *gGFCorrectionKaonMinus=new TGraph();
962                 gGFCorrectionKaonMinus->SetName("gGFCorrectionKaonMinus");
963                 gGFCorrectionKaonMinus->SetTitle("gGFCorrectionKaonMinus");
964                  TString fnameGeanFlukaK="GFCorrection/correctionForCrossSection.321.root";
965                 TFile *fGeanFlukaK= new TFile(fnameGeanFlukaK.Data());
966                 if (!fGeanFlukaK)
967                 {
968                         fnameGeanFlukaK="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/GFCorrection/correctionForCrossSection.321.root";
969                         fGeanFlukaK= new TFile(fnameGeanFlukaK.Data());
970                         if (!fGeanFlukaK)
971                                 return;
972                 }
973                 TH1F *hGeantFlukaKPos=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles");
974                 TH1F *hGeantFlukaKNeg=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles");
975                 //getting GF func for Kaons with TOF
976                 TF1 *fGFKPosTracking;
977                 fGFKPosTracking = TrackingEff_geantflukaCorrection(3,kPositive);
978                 TF1 *fGFKNegTracking;
979                  fGFKNegTracking = TrackingEff_geantflukaCorrection(3,kNegative);
980                  TF1 *fGFKPosMatching;
981                  fGFKPosMatching = TOFmatchMC_geantflukaCorrection(3,kPositive);
982                  TF1 *fGFKNegMatching;
983                  fGFKNegMatching = TOFmatchMC_geantflukaCorrection(3,kNegative);
984                  for(Int_t binK=0;binK<=Spectra[1]->GetNbinsX();binK++)
985                 {
986                         if(Spectra[1]->GetBinCenter(binK)<tofpt)
987                         {//use TPC GeantFlukaCorrection
988                                 Float_t FlukaCorrKPos=hGeantFlukaKPos->GetBinContent(hGeantFlukaKPos->FindBin(Spectra[1]->GetBinCenter(binK)));
989                                 Float_t FlukaCorrKNeg=hGeantFlukaKNeg->GetBinContent(hGeantFlukaKNeg->FindBin(Spectra[4]->GetBinCenter(binK)));
990                         //      Printf("TPC Geant/Fluka: pt:%f  Pos:%f  Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPos,FlukaCorrKNeg);
991                                 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPos);
992                                 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNeg);
993                                 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPos);
994                                 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNeg);
995                                 gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),FlukaCorrKPos);
996                                 gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),FlukaCorrKNeg);
997                         }
998                         else
999                         {
1000                                 gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),0);
1001                                 gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),0);
1002                                 Float_t FlukaCorrKPosTracking=fGFKPosTracking->Eval(Spectra[1]->GetBinCenter(binK));
1003                                 Float_t FlukaCorrKNegTracking=fGFKNegTracking->Eval(Spectra[1]->GetBinCenter(binK));
1004                         //      Printf("TPC/TOF Geant/Fluka Tracking: pt:%f  Pos:%f  Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosTracking,FlukaCorrKNegTracking);
1005                                 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosTracking);
1006                                 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegTracking);
1007                                 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosTracking);
1008                                 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegTracking);
1009                                 Float_t FlukaCorrKPosMatching=fGFKPosMatching->Eval(Spectra[1]->GetBinCenter(binK));
1010                                 Float_t FlukaCorrKNegMatching=fGFKNegMatching->Eval(Spectra[1]->GetBinCenter(binK));
1011                         //      Printf("TPC/TOF Geant/Fluka Matching: pt:%f  Pos:%f  Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosMatching,FlukaCorrKNegMatching);
1012                                 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosMatching);
1013                                 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegMatching);
1014                                 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosMatching);
1015                                 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegMatching);
1016                         }
1017                 }
1018           }
1019           if(!(options&kgeantflukaProton))      
1020                 return;
1021           //Geant Fluka for P in TPC
1022           Printf("\nGF correction for Protons");
1023           const Int_t kNCharge=2;
1024           Int_t kPos=0;
1025           Int_t kNeg=1;
1026           TString fnameGFProtons= "GFCorrection/correctionForCrossSection.root";
1027           TFile* fGFProtons = new TFile (fnameGFProtons.Data());
1028           if (!fGFProtons)
1029           { 
1030                 fnameGFProtons=="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/GFCorrection/correctionForCrossSection.root";
1031                 TFile* fGFProtons = new TFile (fnameGFProtons.Data());
1032                 if (!fGFProtons)
1033                         return;
1034           }
1035                 
1036
1037
1038           TH2D * hCorrFluka[kNCharge];
1039           TH2D * hCorrFluka[2];
1040           hCorrFluka[kPos] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionProtons");
1041           hCorrFluka[kNeg] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionAntiProtons");
1042           //getting GF func for Kaons with TPCTOF
1043           TF1 *fGFpPosTracking;
1044           fGFpPosTracking = TrackingEff_geantflukaCorrection(4,kPositive);
1045           TF1 *fGFpNegTracking;
1046           fGFpNegTracking = TrackingEff_geantflukaCorrection(4,kNegative);
1047           TF1 *fGFpPosMatching;
1048           fGFpPosMatching = TOFmatchMC_geantflukaCorrection(4,kPositive);
1049           TF1 *fGFpNegMatching;
1050           fGFpNegMatching = TOFmatchMC_geantflukaCorrection(4,kNegative);
1051           
1052          
1053                 Int_t nbins = Spectra[2]->GetNbinsX();
1054                 
1055                 for(Int_t ibin = 0; ibin < nbins; ibin++)
1056                 {
1057                         if(Spectra[2]->GetBinCenter(ibin)<tofpt)
1058                         {//use TPC GeantFlukaCorrection
1059                                 for(Int_t icharge = 0; icharge < kNCharge; icharge++)
1060                                 {
1061                                         Int_t nbinsy=hCorrFluka[icharge]->GetNbinsY();
1062                                         Float_t pt = Spectra[2]->GetBinCenter(ibin);
1063                                         Float_t minPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(1);
1064                                         Float_t maxPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1);
1065                                         if (pt < minPtCorrection) 
1066                                                 pt = minPtCorrection+0.0001;
1067                                         if (pt > maxPtCorrection) 
1068                                                 pt = maxPtCorrection;
1069                                         Float_t correction = hCorrFluka[icharge]->GetBinContent(1,hCorrFluka[icharge]->GetYaxis()->FindBin(pt));
1070                                 
1071                                         if (correction > 0.0) 
1072                                         {// If the bin is empty this is a  0
1073                         //                      cout<<icharge<<" "<<ibin<<" "<<correction<<endl;
1074                                                 Spectra[icharge*3+2]->SetBinContent(ibin,Spectra[icharge*3+2]->GetBinContent(ibin)*correction);
1075                                                 Spectra[icharge*3+2]->SetBinError(ibin,Spectra[icharge*3+2]->GetBinError(ibin)*correction);
1076                                         }
1077                                         else if (Spectra[icharge*3+2]->GetBinContent(ibin) > 0.0) 
1078                                         { 
1079                                                 // If we are skipping a non-empty bin, we notify the user
1080                                                 cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for " <<"protons Pt:"<< pt<< endl;
1081                                                 cout << " Bin content: " << Spectra[icharge*3+2]->GetBinContent(ibin)  << endl;
1082                                         }
1083                                 }
1084                         }
1085                         else
1086                         {
1087                                 Float_t FlukaCorrpPosTracking=fGFpPosTracking->Eval(Spectra[2]->GetBinCenter(ibin));
1088                                 Float_t FlukaCorrpNegTracking=fGFpNegTracking->Eval(Spectra[2]->GetBinCenter(ibin));
1089                         //      Printf("TPC/TOF Geant/Fluka Tracking: pt:%f  Pos:%f  Neg:%f",Spectra[2]->GetBinCenter(ibin),FlukaCorrpPosTracking,FlukaCorrpNegTracking);
1090                                 Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*FlukaCorrpPosTracking);
1091                                 Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*FlukaCorrpNegTracking);
1092                                 Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError(ibin)*FlukaCorrpPosTracking);
1093                                 Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError(ibin)*FlukaCorrpNegTracking);
1094                                 Float_t FlukaCorrpPosMatching=fGFpPosMatching->Eval(Spectra[2]->GetBinCenter(ibin));
1095                                 Float_t FlukaCorrpNegMatching=fGFpNegMatching->Eval(Spectra[2]->GetBinCenter(ibin));
1096                         //      Printf("TPC/TOF Geant/Fluka Matching: pt:%f  Pos:%f  Neg:%f",Spectra[2]->GetBinCenter(ibin),FlukaCorrpPosMatching,FlukaCorrpNegMatching);
1097                                 Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*FlukaCorrpPosMatching);
1098                                 Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*FlukaCorrpNegMatching);
1099                                 Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError(ibin)*FlukaCorrpPosMatching);
1100                                 Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError(ibin)*FlukaCorrpNegMatching);
1101                         }               
1102                 }
1103          fGeanFlukaK->Close();
1104          delete fGeanFlukaK;
1105 }
1106
1107
1108 ///////////
1109 TF1 *
1110 TrackingEff_geantflukaCorrection(Int_t ipart, Int_t icharge)
1111 {
1112
1113   if (ipart == 3 && icharge == kNegative) {
1114     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.);
1115     return f;
1116   }
1117   else if (ipart == 4 && icharge == kNegative) {
1118     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.);
1119   }
1120   else
1121     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionNull(x)", 0., 5.);
1122
1123   return f;
1124 }
1125
1126 Double_t
1127 TrackingPtGeantFlukaCorrectionNull(Double_t pTmc)
1128 {
1129   return 1.;
1130 }
1131
1132 Double_t
1133 TrackingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
1134 {
1135   return (1 - 0.129758 *TMath::Exp(-pTmc*0.679612));
1136 }
1137
1138 Double_t
1139 TrackingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
1140 {
1141   return TMath::Min((0.972865 + 0.0117093*pTmc), 1.);
1142 }
1143 ///////////////////////////////////////////
1144 TF1 *
1145 TOFmatchMC_geantflukaCorrection(Int_t ipart, Int_t icharge)
1146 {
1147
1148   if (ipart == 3 && icharge == kNegative) {
1149     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.);
1150     return f;
1151   }
1152   else if (ipart == 4 && icharge == kNegative) {
1153     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.);
1154   }
1155   else
1156     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionNull(x)", 0., 5.);
1157
1158   return f;
1159 }
1160
1161
1162 Double_t
1163 MatchingPtGeantFlukaCorrectionNull(Double_t pTmc)
1164 {
1165   return 1.;
1166 }
1167
1168 Double_t 
1169 MatchingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
1170 {
1171   Float_t ptTPCoutP =pTmc*(1-6.81059e-01*TMath::Exp(-pTmc*4.20094));
1172   return (TMath::Power(1 - 0.129758*TMath::Exp(-ptTPCoutP*0.679612),0.07162/0.03471));
1173 }
1174
1175 Double_t
1176 MatchingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
1177 {
1178   Float_t ptTPCoutK=pTmc*(1- 3.37297e-03/pTmc/pTmc - 3.26544e-03/pTmc);
1179   return TMath::Min((TMath::Power(0.972865 + 0.0117093*ptTPCoutK,0.07162/0.03471)), 1.);
1180 }
1181 void MatchingTOFEff(TH1F** Spectra, TList* list=0x0)
1182 {
1183           if(TOFMatchingScalling[0]<0.0&&TOFMatchingScalling[1]<0.0)
1184           {
1185                 TH1F *hMatcEffPos_data=(TH1F*)tcutsdata->GetHistoNMatchedPos();
1186                   hMatcEffPos_data->Sumw2();
1187                   //hMatcEffPos_data->Divide((TH1F*)tcutsdata->GetHistoNSelectedPos());
1188                   hMatcEffPos_data->Divide(hMatcEffPos_data,(TH1F*)tcutsdata->GetHistoNSelectedPos(),1,1,"B");
1189                   hMatcEffPos_data->SetTitle("Matching Eff Pos - data");
1190                   TH1F *hMatcEffNeg_data=(TH1F*)tcutsdata->GetHistoNMatchedNeg();
1191                   hMatcEffNeg_data->Sumw2();
1192                   //hMatcEffNeg_data->Divide((TH1F*)tcutsdata->GetHistoNSelectedNeg());
1193                   hMatcEffNeg_data->Divide(hMatcEffNeg_data,(TH1F*)tcutsdata->GetHistoNSelectedNeg(),1,1,"B");
1194                   hMatcEffNeg_data->SetTitle("Matching Eff Neg - data");
1195                   TH1F *hMatcEffPos_mc=(TH1F*)tcutsmc->GetHistoNMatchedPos();
1196                   hMatcEffPos_mc->Sumw2();
1197                   //hMatcEffPos_mc->Divide((TH1F*)tcutsmc->GetHistoNSelectedPos());
1198                   hMatcEffPos_mc->Divide(hMatcEffPos_mc,(TH1F*)tcutsmc->GetHistoNSelectedPos(),1,1,"B");
1199                   hMatcEffPos_mc->SetTitle("Matching Eff Pos - mc");
1200                   TH1F *hMatcEffNeg_mc=(TH1F*)tcutsmc->GetHistoNMatchedNeg();
1201                   hMatcEffNeg_mc->Sumw2();
1202                   //hMatcEffNeg_mc->Divide((TH1F*)tcutsmc->GetHistoNSelectedNeg());
1203                   hMatcEffNeg_mc->Divide(hMatcEffNeg_mc,(TH1F*)tcutsmc->GetHistoNSelectedNeg(),1,1,"B");
1204                   hMatcEffNeg_mc->SetTitle("Matching Eff Neg - mc");
1205
1206
1207                   hMatcEffPos_data->Divide(hMatcEffPos_mc);
1208                   hMatcEffNeg_data->Divide(hMatcEffNeg_mc);
1209                   hMatcEffPos_data->SetName("MatchingTOFPos");
1210                   hMatcEffNeg_data->SetName("MatchingTOFNeg");
1211                   
1212                   
1213                   TF1 *pol0MatchPos_data=new TF1("pol0MatchPos_data","pol0",.6,5);
1214                   hMatcEffPos_data->Fit("pol0MatchPos_data","MNR");
1215                   TF1 *pol0MatchNeg_data=new TF1("pol0MatchNeg_data","pol0",.6,5);
1216                   hMatcEffNeg_data->Fit("pol0MatchNeg_data","MNR");
1217                 
1218                 list->Add(hMatcEffPos_data);
1219                   list->Add(hMatcEffNeg_data);
1220                 
1221                         
1222                   TOFMatchingScalling[0]=pol0MatchPos_data->GetParameter(0);
1223                   TOFMatchingScalling[1]=pol0MatchNeg_data->GetParameter(0);
1224
1225           }
1226           //Correction spectra for matching efficiency
1227           //For the moment I'm using the inclusive correction
1228           for(Int_t ipart=0;ipart<3;ipart++)
1229           {
1230                   
1231                 for(Int_t ibin=1;ibin<Spectra[ipart]->GetNbinsX();ibin++)
1232                 {
1233                         Float_t ptspectra=Spectra[ipart]->GetBinCenter(ibin);
1234                         if(ptspectra<tcutsdata->GetPtTOFMatching())
1235                                 continue;
1236                   //Spectra[ipart]->SetBinContent(ibin,( Spectra[ipart]->GetBinContent(ibin)/hMatcEffPos_data->GetBinContent(hMatcEffPos_data->FindBin(ptspectra))));
1237                   //Spectra[ipart+3]->SetBinContent(ibin,( Spectra[ipart+3]->GetBinContent(ibin)/hMatcEffNeg_data->GetBinContent(hMatcEffNeg_data->FindBin(ptspectra))));
1238                         Spectra[ipart]->SetBinContent(ibin,( Spectra[ipart]->GetBinContent(ibin)/TOFMatchingScalling[0]));
1239                         Spectra[ipart+3]->SetBinContent(ibin,( Spectra[ipart+3]->GetBinContent(ibin)/TOFMatchingScalling[1]));
1240                 }
1241           }
1242 }
1243 Double_t eta2y(Double_t pt, Double_t mass, Double_t eta)
1244 {
1245   Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
1246   return TMath::ASinH(pt / mt * TMath::SinH(eta));
1247 }
1248
1249 TH1* GetSumAllCh(TH1F** spectra, Double_t* mass  )
1250 {
1251         TH1F* allch=(((TH1F*))spectra[0]->Clone("allCh"));
1252         allch->Reset();
1253         for (int i=0;i<6;i++)
1254         {
1255                 Double_t masstmp=mass[i%3];
1256                 for (int j=1;j<spectra[i]->GetXaxis()->GetNbins();j++)
1257                 {
1258                         Float_t value=spectra[i]->GetBinContent(j);
1259                         Float_t error=spectra[i]->GetBinError(j);
1260                         if(value>0.0)
1261                         {
1262                                 Float_t pt=spectra[i]->GetXaxis()->GetBinCenter(j);
1263                                 Float_t mt2=masstmp*masstmp+pt*pt;              
1264                                 Float_t mt=TMath::Sqrt(mt2);
1265                                 Float_t maxy=eta2y(pt,masstmp,0.8);
1266                                 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))));
1267                                 conver=conver/1.6;
1268                                 //cout<<maxy<<" "<<conver<<" "<<masstmp<<""<<spectra[i]->GetName()<<endl;
1269                                 Float_t bincontent=allch->GetBinContent(j);
1270                                 Float_t binerror=allch->GetBinError(j);
1271                                 bincontent+=conver*value;
1272                                 binerror=TMath::Sqrt(binerror*binerror+conver*conver*error*error);
1273                                 allch->SetBinContent(j,bincontent);
1274                                 allch->SetBinError(j,binerror);
1275                         }
1276                         
1277                 }
1278         }
1279         Divideby2pipt(allch);
1280         allch->SetTitle("N_{ch};p_{T} (GeV/c);1/(2 #pi p_{T})dN/p_{T} |#eta|<0.8");             
1281         return allch;
1282 }
1283
1284 void Divideby2pipt(TH1* hist)
1285 {
1286
1287         for (int i=1;i<hist->GetXaxis()->GetNbins();i++)
1288         {
1289                 Float_t value=hist->GetBinContent(i);
1290                 Float_t error=hist->GetBinError(i);
1291                 Float_t pt=hist->GetXaxis()->GetBinCenter(i);
1292                 hist->SetBinContent(i,value/(2.0*TMath::Pi()*pt));
1293                 hist->SetBinError(i,error/(2.0*TMath::Pi()*pt));
1294
1295         }
1296 }
1297
1298 Short_t DCAfitsettings (Float_t pt, Int_t type)
1299 {
1300         Short_t value=0x0;
1301         if (pt<maxptformaterial[type]&&pt>minptformaterial[type])
1302                 value=value+2;
1303         if (pt<maxptforWD[type]&&pt>minptforWD[type])
1304                 value=value+1;
1305         return value;   
1306
1307
1308
1309 Float_t Normaliztionwithbin0integrals(UInt_t options)
1310 {
1311         
1312         TH1F* bin0mcRec=(TH1F*)ecutsmc->GetHistoVtxGenerated()->Clone("Bin0_rec");
1313         TH1F* bin0mcMC=(TH1F*)ecutsmc->GetHistoVtxGenerated()->Clone("Bin0_MC");
1314
1315         TH1F* vertexmc=ecutsmc->GetHistoVtxAftSelwithoutZvertexCut(); 
1316         TH1F* vertexmcMCz=ecutsmc->GetHistoVtxAftSelwithoutZvertexCutusingMCz(); 
1317         TH1F* vertexdata=ecutsdata->GetHistoVtxAftSelwithoutZvertexCut();
1318
1319         TH1I* histodata=ecutsdata->GetHistoCuts();
1320         TH1I* histomc=ecutsmc->GetHistoCuts();
1321
1322         Float_t dataevents=(Float_t)histodata->GetBinContent(3);
1323         //cout<<histodata->GetBinContent(2)<<endl;
1324         Float_t databin0events=((Float_t)histodata->GetBinContent(2))-((Float_t)histodata->GetBinContent(4));   
1325
1326         bin0mcRec->Sumw2();
1327         bin0mcMC->Sumw2();
1328                 
1329         bin0mcRec->Add(vertexmc,-1);
1330         bin0mcMC->Add(vertexmcMCz,-1);
1331         
1332         bin0mcRec->Divide(vertexmc);
1333         bin0mcMC->Divide(vertexmcMCz);
1334         
1335         bin0mcRec->Multiply(vertexdata);
1336         bin0mcMC->Multiply(vertexdata);
1337         
1338         Float_t bin0mcRecN=0.0;
1339         Float_t bin0mcMCN=0.0;
1340
1341         for (int i=0;i<=bin0mcRec->GetXaxis()->GetNbins();i++)
1342         {
1343                 bin0mcRecN+=bin0mcRec->GetBinContent(i);
1344                 bin0mcMCN+=bin0mcMC->GetBinContent(i);
1345
1346         }
1347         bin0mcRec->Scale(databin0events/bin0mcRecN);
1348         bin0mcMC->Scale(databin0events/bin0mcMCN);              
1349         
1350         Int_t binmin=bin0mcRec->GetXaxis()->FindBin(-10);
1351         Int_t binmax=bin0mcRec->GetXaxis()->FindBin(10)-1;
1352         cout<<  bin0mcRec->GetXaxis()->GetBinLowEdge(binmin)<<" "<<bin0mcRec->GetXaxis()->GetBinUpEdge(binmax)<<endl;
1353         cout<<bin0mcRecN<<" "<<bin0mcMCN<<" "<<databin0events<<endl;    
1354         cout<<dataevents<<" normalization "<<dataevents+bin0mcRec->Integral(binmin,binmax)<<" "<<dataevents+bin0mcMC->Integral(binmin,binmax)<<endl;
1355         cout<<histodata->GetBinContent(2)<<" "<<histodata->GetBinContent(4)<<endl;
1356         if ((options&knormalizationwithbin0integralsdata)==knormalizationwithbin0integralsdata)
1357                 return  dataevents+bin0mcRec->Integral(binmin,binmax);
1358         else if ((options&kknormalizationwithbin0integralsMC)==knormalizationwithbin0integralsMC)   
1359                 return dataevents+bin0mcMC->Integral(binmin,binmax) ;
1360         else
1361                 return 1;               
1362 }
1363  
1364
1365 Bool_t ReadConfigFile(TString configfile)
1366 {
1367         ifstream infile(configfile.Data());
1368         if(infile.is_open()==false)
1369                 return false;
1370         TString namesofSetting[35]={"CutRangeMin","CutRangeMax","FitRangeMin","FitRangeMax","MinMatPionPlus","MaxMatPionPlus","MinMatKaonPlus","MaxMatKaonPlus","MinMatProtonPlus","MaxMatProtonPlus","MinMatPionMinus","MaxMatPionMinus","MinMatKaonMinus","MaxMatKaonMinus","MinMatProtonMinus","MaxMatProtonMinus","MinWDPionPlus","MaxWDPionPlus","MinWDKaonPlus","MaxWDKaonPlus","MinWDProtonPlus","MaxWDProtonPlus","MinWDPionMinus","MaxWDPionMinus","MinWDKaonMinus","MaxWDKaonMinus","MinWDProtonMinus","MaxWDProtonMinus","MaxContaminationPIDMC","MinPions","MaxPions","MinKaons","MaxKaons","MinProtons","MaxProtons"};     
1371
1372         char buffer[256];
1373         while (infile.eof()==false)
1374         {
1375                 buffer[0]='#'; 
1376                 while (buffer[0]=='#'&&infile.eof()==false)
1377                         infile.getline(buffer,256);
1378                 TString tmpstring(buffer);
1379                 cout<<buffer<<endl;
1380                 if(tmpstring.Contains(namesofSetting[0]))
1381                         CutRange[0]=(tmpstring.Remove(0,namesofSetting[0].Length()+1)).Atof();
1382                 else if (tmpstring.Contains(namesofSetting[1]))
1383                         CutRange[1]=(tmpstring.Remove(0,namesofSetting[1].Length()+1)).Atof();  
1384                 else if (tmpstring.Contains(namesofSetting[2]))
1385                         FitRange[0]=(tmpstring.Remove(0,namesofSetting[2].Length()+1)).Atof();
1386                 else if (tmpstring.Contains(namesofSetting[3]))
1387                         FitRange[1]=(tmpstring.Remove(0,namesofSetting[3].Length()+1)).Atof();
1388                 else if (tmpstring.Contains(namesofSetting[4]))
1389                         minptformaterial[0]=(tmpstring.Remove(0,namesofSetting[4].Length()+1)).Atof();          
1390                 else if (tmpstring.Contains(namesofSetting[5]))
1391                         maxptformaterial[0]=(tmpstring.Remove(0,namesofSetting[5].Length()+1)).Atof();  
1392                 else if (tmpstring.Contains(namesofSetting[6]))
1393                         minptformaterial[1]=(tmpstring.Remove(0,namesofSetting[6].Length()+1)).Atof();          
1394                 else if (tmpstring.Contains(namesofSetting[7]))
1395                         maxptformaterial[1]=(tmpstring.Remove(0,namesofSetting[7].Length()+1)).Atof();  
1396                 else if (tmpstring.Contains(namesofSetting[8]))
1397                         minptformaterial[2]=(tmpstring.Remove(0,namesofSetting[8].Length()+1)).Atof();          
1398                 else if (tmpstring.Contains(namesofSetting[9]))
1399                         maxptformaterial[2]=(tmpstring.Remove(0,namesofSetting[9].Length()+1)).Atof();
1400                 else if (tmpstring.Contains(namesofSetting[10]))
1401                         minptformaterial[3]=(tmpstring.Remove(0,namesofSetting[10].Length()+1)).Atof();         
1402                 else if (tmpstring.Contains(namesofSetting[11]))
1403                         maxptformaterial[3]=(tmpstring.Remove(0,namesofSetting[11].Length()+1)).Atof();
1404                 else if (tmpstring.Contains(namesofSetting[12]))
1405                         minptformaterial[4]=(tmpstring.Remove(0,namesofSetting[12].Length()+1)).Atof();         
1406                 else if (tmpstring.Contains(namesofSetting[13]))
1407                         maxptformaterial[4]=(tmpstring.Remove(0,namesofSetting[13].Length()+1)).Atof();
1408                 else if (tmpstring.Contains(namesofSetting[14]))
1409                         minptformaterial[5]=(tmpstring.Remove(0,namesofSetting[14].Length()+1)).Atof();         
1410                 else if (tmpstring.Contains(namesofSetting[15]))
1411                         maxptformaterial[5]=(tmpstring.Remove(0,namesofSetting[15].Length()+1)).Atof();
1412                 else if (tmpstring.Contains(namesofSetting[16]))
1413                         minptforWD[0]=(tmpstring.Remove(0,namesofSetting[16].Length()+1)).Atof();               
1414                 else if (tmpstring.Contains(namesofSetting[17]))
1415                         maxptforWD[0]=(tmpstring.Remove(0,namesofSetting[17].Length()+1)).Atof();       
1416                 else if (tmpstring.Contains(namesofSetting[18]))
1417                         minptforWD[1]=(tmpstring.Remove(0,namesofSetting[18].Length()+1)).Atof();               
1418                 else if (tmpstring.Contains(namesofSetting[19]))
1419                         maxptforWD[1]=(tmpstring.Remove(0,namesofSetting[19].Length()+1)).Atof();       
1420                 else if (tmpstring.Contains(namesofSetting[20]))
1421                         minptforWD[2]=(tmpstring.Remove(0,namesofSetting[20].Length()+1)).Atof();               
1422                 else if (tmpstring.Contains(namesofSetting[21]))
1423                         maxptforWD[2]=(tmpstring.Remove(0,namesofSetting[21].Length()+1)).Atof();
1424                 else if (tmpstring.Contains(namesofSetting[22]))
1425                         minptforWD[3]=(tmpstring.Remove(0,namesofSetting[22].Length()+1)).Atof();               
1426                 else if (tmpstring.Contains(namesofSetting[23]))
1427                         maxptforWD[3]=(tmpstring.Remove(0,namesofSetting[23].Length()+1)).Atof();
1428                 else if (tmpstring.Contains(namesofSetting[24]))
1429                         minptforWD[4]=(tmpstring.Remove(0,namesofSetting[24].Length()+1)).Atof();               
1430                 else if (tmpstring.Contains(namesofSetting[25]))
1431                         maxptforWD[4]=(tmpstring.Remove(0,namesofSetting[25].Length()+1)).Atof();
1432                 else if (tmpstring.Contains(namesofSetting[26]))
1433                         minptforWD[5]=(tmpstring.Remove(0,namesofSetting[26].Length()+1)).Atof();               
1434                 else if (tmpstring.Contains(namesofSetting[27]))
1435                         maxptforWD[5]=(tmpstring.Remove(0,namesofSetting[27].Length()+1)).Atof();                               
1436                 else if (tmpstring.Contains(namesofSetting[28]))
1437                         fMaxContaminationPIDMC=(tmpstring.Remove(0,namesofSetting[28].Length()+1)).Atof();              
1438                 else if (tmpstring.Contains(namesofSetting[29]))
1439                         minRanges[0]=(tmpstring.Remove(0,namesofSetting[29].Length()+1)).Atof();
1440                 else if (tmpstring.Contains(namesofSetting[30]))
1441                         maxRanges[0]=(tmpstring.Remove(0,namesofSetting[30].Length()+1)).Atof();                
1442                 else if (tmpstring.Contains(namesofSetting[31]))
1443                         minRanges[1]=(tmpstring.Remove(0,namesofSetting[31].Length()+1)).Atof();
1444                 else if (tmpstring.Contains(namesofSetting[32]))
1445                         maxRanges[1]=(tmpstring.Remove(0,namesofSetting[32].Length()+1)).Atof();                
1446                 else if (tmpstring.Contains(namesofSetting[33]))
1447                         minRanges[2]=(tmpstring.Remove(0,namesofSetting[33].Length()+1)).Atof();
1448                 else if (tmpstring.Contains(namesofSetting[34]))
1449                         maxRanges[2]=(tmpstring.Remove(0,namesofSetting[34].Length()+1)).Atof();                
1450                 else
1451                         continue;
1452
1453
1454 //      Double_t minRanges[3]={0.3,0.3,0.45};
1455 //Double_t maxRanges[3]={1.5,1.2,2.2};
1456 //Double_t fMaxContaminationPIDMC=0.2;
1457
1458
1459
1460         }
1461         for(int i=0;i<6;i++)
1462                 cout<<minptformaterial[i]<<" "<<maxptformaterial[i]<<" "<<minptforWD[i]<<" "<<maxptforWD[i]<<endl;
1463         cout<<FitRange[0]<<" "<<FitRange[1]<<" "<<CutRange[0]<<CutRange[1]<<endl;
1464         if(FitRange[0]>=FitRange[1])    
1465         {
1466                 cout<<"A"<<endl;                                
1467                 return false;
1468         }
1469         if(CutRange[0]>=CutRange[1])
1470         {       
1471                 cout<<"B"<<endl;                                
1472                 return false;
1473         }
1474         for(int i=0;i<6;i++)
1475         {
1476                 if((minptformaterial[i]>maxptformaterial[i]&&minptformaterial[i]>0.0)||minptformaterial[i]<0.0||maxptformaterial[i]<0.0)
1477                 {
1478                         cout<<"C"<<endl;
1479                         return false;
1480                 }
1481                 if((minptforWD[i]>maxptforWD[i]&&minptforWD[i]>0.0)||minptforWD[i]<0.0||maxptforWD[i]<0.0)
1482                 {
1483                         cout<<"D"<<endl;
1484                         return false;
1485                 }
1486         }
1487         for(int i=0;i<3;i++)
1488                 if(minRanges[i]>maxRanges[i])
1489                         return false;
1490                 
1491
1492         return true;
1493 }
1494
1495 void SubHistWithFullCorr(TH1F* h1, TH1F* h2, Float_t factor1=1.0, Float_t factor2=1.0)
1496 {
1497         if(h1->GetNbinsX()!=h2->GetNbinsX())
1498                 return;
1499         for (int i=0;i<=h1->GetNbinsX();i++)
1500         {
1501                 Float_t tmpvalue=factor1*h1->GetBinContent(i)-factor2*h2->GetBinContent(i);
1502                 Float_t tmperror=TMath::Abs(factor1*factor1*h1->GetBinError(i)*h1->GetBinError(i)-factor2*factor2*h2->GetBinError(i)*h2->GetBinError(i));
1503                 h1->SetBinContent(i,tmpvalue);
1504                 h1->SetBinError(i,TMath::Sqrt(tmperror));
1505         }               
1506         
1507 }
1508
1509 void TOFMatchingForNch(TH1* h)
1510 {
1511          if(TOFMatchingScalling[0]>0.0&&TOFMatchingScalling[1]>0.0)
1512          {
1513                 Float_t factor=0.5*TOFMatchingScalling[0]+0.5*TOFMatchingScalling[1];
1514                 for(Int_t ibin=1;ibin<h->GetNbinsX();ibin++)
1515                 {
1516                         Float_t ptspectra=h->GetBinCenter(ibin);
1517                         if(ptspectra<tcutsdata->GetPtTOFMatching())
1518                                 continue;
1519                         h->SetBinContent(ibin,(h->GetBinContent(ibin)/factor));
1520                 }
1521
1522          }
1523           else
1524                 return;                 
1525
1526
1527 }