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