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