]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/SPECTRA/PiKaPr/TestAOD/AnalysisBoth.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
TENDER becomes Tender
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / AnalysisBoth.C
CommitLineData
1af4558e 1
88d02ae4 2class AliSpectraBothHistoManager;
3class AliSpectraBothEventCuts;
4class AliSpectraBothTrackCuts;
5TString Charge[]={"Pos","Neg"};
6TString Sign[]={"Plus","Minus"};
7TString Particle[]={"Pion","Kaon","Proton"};
c0a5ac55 8TString symboles[]={"#pi^{+}","K^{+}","p","pi^{-}","K^{-}","#bar{p}"};
88d02ae4 9AliSpectraBothHistoManager* managerdata=0x0;
10AliSpectraBothEventCuts* ecutsdata=0x0;
11AliSpectraBothTrackCuts* tcutsdata=0x0;
12
13AliSpectraBothHistoManager* managermc=0x0;
14AliSpectraBothEventCuts* ecutsmc=0x0;
15AliSpectraBothTrackCuts* tcutsmc=0x0;
16
17Float_t TOFMatchingScalling[2]={-1,-1};
18Int_t Color[3]={1,2,4};
19Int_t Marker[6]={20,21,22,24,25,26};
20Double_t Range[3]={0.3,0.3,0.5}; // LowPt range for pi k p
21
a4ca64e3 22
23Double_t FitRange[2]={-2.0,2.0};
24Double_t CutRange[2]={-3.0,3.0};
25Double_t minptformaterial[6]={0.0,0.2,0.0,0.0,0.2,0.0};
26Double_t maxptformaterial[6]={0.0,0.6,1.3,0.0,0.6,0.0};
27Double_t minptforWD[6]={0.2,100.0,0.3,0.2,100.0,0.3};
28Double_t maxptforWD[6]={1.5,-100.0,2.0,1.5,-100.0,2.0};
1af4558e 29Double_t minRanges[3]={0.3,0.3,0.45};
30Double_t maxRanges[3]={1.5,1.2,2.2};
cd28560b 31Double_t TOFPIDsignalmatching[]={-1.0,-1.0,-1.0};
1af4558e 32Double_t fMaxContaminationPIDMC=0.2;
93531558 33TString filenames[]={"eff.root","pid.root","sec.root"};
a4ca64e3 34
88d02ae4 35enum ECharge_t {
36 kPositive,
37 kNegative,
38 kNCharges
39};
b3ea73e1 40enum {
a4ca64e3 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
b6a74f8a 51 kuserangeonfigfile=0x400, // use of config file for dca fit settings
1af4558e 52 kskipconcutonspectra=0x800, //do not use conPID<02 cut useful for syst. studies
cd28560b 53 kuseTOFmatchingcorrection=0x1000, // if set tof matching correction is applied.
93531558 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
b3ea73e1 59};
88d02ae4 60
0c1898ba 61Bool_t OpenFile(TString dirname, TString outputname, Bool_t mcflag,Bool_t mcasdata=false);
a4ca64e3 62void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TString outnamemc="",TString configfile="" )
88d02ae4 63{
57a75b31 64gStyle->SetOptStat(0);
88d02ae4 65 TH1::AddDirectory(kFALSE);
66 gSystem->Load("libCore.so");
67 gSystem->Load("libPhysics.so");
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");
af472fff 76 gSystem->Load("libTender");
88d02ae4 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
8cec8e63 87 TFormula* dcacutxy=0x0;
88 if(!(options&kdonotusedcacuts))
89 {
90
91 AliESDtrackCuts* esdtrackcuts= AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE);
a4ca64e3 92 TString formulastring(esdtrackcuts->GetMaxDCAToVertexXYPtDep());
8cec8e63 93 formulastring.ReplaceAll("pt","x");
94 dcacutxy=new TFormula("dcacutxy",formulastring.Data());
95 }
a4ca64e3 96 if(options&kuserangeonfigfile)
97 if(!ReadConfigFile(configfile))
98 return;
88d02ae4 99 TList* lout=new TList();
a4ca64e3 100
88d02ae4 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);
0c1898ba 111 OpenFile(indirname,outnamedata,false,((Bool_t)(options&kmcisusedasdata)));
88d02ae4 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];
2245e1bd 126 TH1F* contPIDpri[6];
b6a74f8a 127 TH1F* contSecMC[6];
88d02ae4 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
c21424b2 139 TH1F* corrLeonardo[6];
140 TH1F* doubleconuntsdata[3];
141 TH1F* doubleconuntsMC[3];
88d02ae4 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);
2245e1bd 151 GetPtHistFromPtDCAhisto("hHistPtRecSigmaPrimary","conPIDprimary",managermc,contPIDpri,dcacutxy);
152
88d02ae4 153
b3ea73e1 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
3a71f081 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();
b3ea73e1 165 if(options&knormalizationtoeventspassingPhySel)
166 {
3a71f081 167 //neventsmcall= ecutsmc->NumberOfProcessedEvents();
b3ea73e1 168 neventsdata=ecutsdata->NumberOfPhysSelEvents();
169 neventsmc=ecutsmc->NumberOfPhysSelEvents();
170 }
3a71f081 171 else if ((options&knormalizationwithbin0integralsdata)||(options&knormalizationwithbin0integralsMC))
172 {
173 neventsdata=Normaliztionwithbin0integrals(options);
174 neventsmc=ecutsmc->NumberOfPhysSelEvents();
175 }
b3ea73e1 176 else
177 {
50c532c4 178 neventsdata=ecutsdata->NumberOfEvents(); //number of accepted events
b3ea73e1 179 neventsmc=ecutsmc->NumberOfEvents();
180 neventsmcall= ecutsmc->NumberOfEvents();
181
b3ea73e1 182 }
88d02ae4 183 GetMCTruth(MCTruth);
50c532c4 184 cout<<neventsdata<<" Events"<<endl;
57a75b31 185 cout<< neventsmc<<" Events "<<endl;
93531558 186 if(neventsdata<1)
57a75b31 187 {
93531558 188 cout<<"No events DATA"<<endl;
57a75b31 189 return;
190 }
93531558 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;
88d02ae4 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");
a4ca64e3 214 spectraall->SetTitle("recNch");
88d02ae4 215 TH1F* contall=(TH1F*)allrecMC->Clone("contall");
a4ca64e3 216 contall->SetTitle("contall");
50c532c4 217 //contall->Add(alleff,-1);
218 SubHistWithFullCorr(contall,alleff);
88d02ae4 219 alleff->Divide(alleff,allgen,1,1,"B");
b3ea73e1 220 contall->Divide(contall,allrecMC,1,1,"B");
88d02ae4 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);
1e6bff53 236 lout->Add(contall);
237
88d02ae4 238 for (int i=0;i<6;i++)
239 {
240
241
242 TString tmpname(rawspectramc[i]->GetTitle());
243 tmpname.ReplaceAll("SpectraMC","%s");
a4ca64e3 244 contallMC[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contallMC"));
245 contallMC[i]->SetTitle(Form(tmpname.Data(),"contallMC"));
88d02ae4 246 contfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contfit"));
a4ca64e3 247 contfit[i]->SetTitle(Form(tmpname.Data(),"contfit"));
88d02ae4 248 contWDfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contWDfit"));
a4ca64e3 249 contWDfit[i]->SetTitle(Form(tmpname.Data(),"contWDfit"));
88d02ae4 250 contMatfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contMatfit"));
a4ca64e3 251 contMatfit[i]->SetTitle(Form(tmpname.Data(),"contMatfit"));
88d02ae4 252 primaryfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"primaryfit"));
a4ca64e3 253 primaryfit[i]->SetTitle(Form(tmpname.Data(),"primaryfit"));
88d02ae4 254 contfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contfitonMC"));
a4ca64e3 255 contfit[i+6]->SetTitle(Form(tmpname.Data(),"contfitonMC"));
88d02ae4 256 contWDfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contWDfitonMC"));
a4ca64e3 257 contWDfit[i+6]->SetTitle(Form(tmpname.Data(),"contWDfitonMC"));
88d02ae4 258 contMatfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contMatfitonMC"));
a4ca64e3 259 contMatfit[i+6]->SetTitle(Form(tmpname.Data(),"contMatfitonMC"));
88d02ae4 260 primaryfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"primaryfitMC"));
a4ca64e3 261 primaryfit[i+6]->SetTitle(Form(tmpname.Data(),"primaryfitMC"));
88d02ae4 262
c21424b2 263
264
265
88d02ae4 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"));
a4ca64e3 280 spectra[i]->SetTitle(Form(tmpname.Data(),"SpectraFinal"));
88d02ae4 281 spectraLeonardo[i]=(TH1F*)rawspectradata[i]->Clone(Form(tmpname.Data(),"SpectraFinalLeonardo"));
a4ca64e3 282 spectraLeonardo[i]->SetTitle(Form(tmpname.Data(),"SpectraFinalLeonardo"));
283
88d02ae4 284 corrLeonardo[i]=(TH1F*)MCTruth[i]->Clone(Form(tmpname.Data(),"CorrFactLeonardo"));
a4ca64e3 285 corrLeonardo[i]->SetTitle(Form(tmpname.Data(),"CorrFactLeonardo"));
88d02ae4 286 corrLeonardo[i]->Divide(corrLeonardo[i],rawspectramc[i],1,1,"B");
287
288
93531558 289 if((options&kusePIDcontaminatiofromfile)==kusePIDcontaminatiofromfile)
290 {
291 CopyCorrectionFromFile(filenames[1],"conPID",contPID);
292 CopyCorrectionFromFile(filenames[1],"conPIDprimary",contPIDpri);
50c532c4 293
93531558 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");
88d02ae4 313
314
315 contPID[i]->Sumw2();
316 rawspectramc[i]->Sumw2();
50c532c4 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]);
88d02ae4 322 contPID[i]->ResetStats();
323 contPID[i]->Sumw2();
324 contPID[i]->Divide(contPID[i],rawspectramc[i],1,1,"B");
2245e1bd 325
93531558 326 if((options&kuseprimaryPIDcont)==kuseprimaryPIDcont)
2245e1bd 327 confinal[i]=(TH1F*)contPIDpri[i]->Clone(Form(tmpname.Data(),"confinal"));
328 else
329 confinal[i]=(TH1F*)contPID[i]->Clone(Form(tmpname.Data(),"confinal"));
a4ca64e3 330 confinal[i]->SetTitle(Form(tmpname.Data(),"confinal"));
88d02ae4 331
b6a74f8a 332 contSecMC[i]=(TH1F*)contWD[i]->Clone(Form(tmpname.Data(),"conSecMC"));
333 contSecMC[i]->Add(contMat[i]);
88d02ae4 334 contWD[i]->Divide(contWD[i],rawspectramc[i],1,1,"B");
335 contMat[i]->Divide(contMat[i],rawspectramc[i],1,1,"B");
b6a74f8a 336 contSecMC[i]->Divide(contSecMC[i],rawspectramc[i],1,1,"B");
c21424b2 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
88d02ae4 348
c21424b2 349 }
350
351
88d02ae4 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]);
1e6bff53 371 lout->Add(primaryfit[i]);
88d02ae4 372 lout->Add(contfit[i+6]);
373 lout->Add(contWDfit[i+6]);
374 lout->Add(contMatfit[i+6]);
1e6bff53 375 lout->Add(primaryfit[i+6]);
88d02ae4 376 lout->Add(spectra[i]);
377 lout->Add(spectraLeonardo[i]);
378 lout->Add(confinal[i]);
2245e1bd 379 lout->Add(contPIDpri[i]);
b6a74f8a 380 lout->Add(contSecMC[i]);
c21424b2 381 if(i>2)
382 {
383 lout->Add(doubleconuntsdata[i-3]);
384 lout->Add(doubleconuntsMC[i-3]);
385 }
88d02ae4 386 }
2245e1bd 387 outdate.ReplaceAll("/","_");
a4ca64e3 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");
93531558 394 if (((options&kdodca)==kdodca)&&((options&kuseseccontaminatiofromfile)!=kuseseccontaminatiofromfile))
88d02ae4 395 DCACorrectionMarek(managerdata,managermc,dcacutxy,fout,contfit,contWDfit,contMatfit,primaryfit);
93531558 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
88d02ae4 408 for (int i=0;i<6;i++)
409 {
93531558 410 if(((options&kdodca)==kdodca)||((options&kuseseccontaminatiofromfile)==kuseseccontaminatiofromfile))
88d02ae4 411 {
2245e1bd 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]);
88d02ae4 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]);
b6a74f8a 421 if(options&kskipconcutonspectra)
422 continue;
2245e1bd 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]);
cd28560b 432 }
433 // Apply correction for wrongly simulated TOF signal in MC
434 if(options&kuseTOFcorrforPIDsignalmatching)
435 TOFPIDsignalmatchingApply(spectra[i],TOFPIDsignalmatching[i%3]);
436
88d02ae4 437 }
438
b3ea73e1 439 GFCorrection(spectra,tcutsdata->GetPtTOFMatching(),options);
440 GFCorrection(spectraLeonardo,tcutsdata->GetPtTOFMatching(),options);
ab0298c1 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 }
cd28560b 450
ab0298c1 451 TH1F* allch=GetSumAllCh(spectra,mass,etacut);
88d02ae4 452 lout->Add(allch);
1af4558e 453 if(options&kuseTOFmatchingcorrection)
454 {
455 MatchingTOFEff(spectra,lout);
456 MatchingTOFEff(spectraLeonardo);
457 TOFMatchingForNch(spectraall);
458 }
88d02ae4 459// lout->ls();
460 fout->cd();
461 TList* listqa=new TList();
33d374d5 462 TList* canvaslist=new TList();
b3ea73e1 463 Float_t vertexcorrection=1.0;
3a71f081 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);
b3ea73e1 469 if (options&kveretxcorrectionandbadchunkscorr)
470 vertexcorrection=corrbadchunksvtx;
88d02ae4 471 cout<<" VTX corr="<<vertexcorrection<<endl;
b3ea73e1 472 if(TMath::Abs(ycut)>0.0)
473 vertexcorrection=vertexcorrection/(2.0*ycut);
88d02ae4 474 for (int i=0;i<6;i++)
475 {
476 spectra[i]->Scale(vertexcorrection);
477 spectraLeonardo[i]->Scale(vertexcorrection);
890fd855 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 }
88d02ae4 484 }
485 allch->Scale(vertexcorrection);
ab0298c1 486 spectraall->Scale(vertexcorrection/etacut);
88d02ae4 487
488 //spectraall->Scale(1.0/1.6);
489 lout->Write("output",TObject::kSingleKey);
490 listqa->Write("outputQA",TObject::kSingleKey);
33d374d5 491 canvaslist->Write("outputcanvas",TObject::kSingleKey);
492
88d02ae4 493 fout->Close();
3a71f081 494 //Normaliztionwithbin0integrals();
88d02ae4 495
496}
497
0c1898ba 498Bool_t OpenFile(TString dirname,TString outputname, Bool_t mcflag, Bool_t mcasdata)
88d02ae4 499{
0c1898ba 500
88d02ae4 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());
0c1898ba 510 if(mcasdata)
511 {
512 cout<<"using MC as data "<<endl;
513 sname=Form("OutputBothSpectraTask_%s_%s","MC",outputname.Data());
514 }
88d02ae4 515 file->ls();
516 TDirectoryFile *dir=(TDirectoryFile*)file->Get(sname.Data());
517 if(!dir)
518 {
ab13cf18 519 // cout<<"no dir "<<sname.Data()<<endl;
520 if(mcasdata)
0c1898ba 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());
88d02ae4 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;
b6a74f8a 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;
88d02ae4 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;
b6a74f8a 557 if(managerdata->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()!=ecutsdata->GetHistoCuts()->GetBinContent(3))
8af77493 558 cout<<"Please check DATA file possible problem with merging"<<" "<<managerdata->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()<<" "<<ecutsdata->GetHistoCuts()->GetBinContent(3)<<endl;
a4ca64e3 559
88d02ae4 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
580TH1F* 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);
50c532c4 592 TH1F* dcahist=(TH1F*)hman->GetDCAHistogram1D(name.Data(),lowedge,lowedge);
1af4558e 593// Float_t inyield=dcahist->Integral(dcahist->GetXaxis()->FindBin(-1.0*cut),dcahist->GetXaxis()->FindBin(cut));
50c532c4 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 }
1af4558e 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;
50c532c4 603
1af4558e 604// cout<<testyield<<" "<<TMath::Sqrt(testerror)<<" error2 "<<inyield<<" "<<TMath::Sqrt(inyield)<<endl;
50c532c4 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
88d02ae4 611 }
612 }
613 histo->Sumw2();
614 return histo;
615}
616
617
618
619
620void GetPtHistFromPtDCAhisto(TString hnamein, TString hnameout, AliSpectraBothHistoManager* hman,TH1F** histo,TFormula* dcacutxy)
621{
1af4558e 622 //Float_t min[3]={0.3,0.3,0.45};
623 //Float_t max[3]={1.5,1.2,2.2};
88d02ae4 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;
a4ca64e3 629 Printf("Getting %s",hnamein.Data());
88d02ae4 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);
1af4558e 635 CleanHisto(histo[index],minRanges[ipart],maxRanges[ipart]);
88d02ae4 636 }
637 }
638}
639void 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 {
1af4558e 650 if(contpid->GetBinContent(i)>fMaxContaminationPIDMC)
88d02ae4 651 {
652 h->SetBinContent(i,0);
653 h->SetBinError(i,0);
654 }
655 }
656 }
657}
658
659
660void 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");
a4ca64e3 663
2245e1bd 664
88d02ae4 665 Printf("\DCACorr");
88d02ae4 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
88d02ae4 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);
8cec8e63 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 }
88d02ae4 688 debug<<"cut "<<CutRange[1]<<" "<<CutRange[0]<<endl;
2245e1bd 689 Short_t fitsettings=DCAfitsettings(lowedge+0.5*binwidth,index);
690 debug<<"settings "<< fitsettings<<endl;
691 if(fitsettings==0)
692 continue;
693
88d02ae4 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);
c0a5ac55 695 cDCA->SetMargin(0.1,0.02,0.1,0.02);
b6a74f8a 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
88d02ae4 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();
57a75b31 710 Double_t minentries=100;
2245e1bd 711 debug<<hToFit->GetEntries()<<" "<<hmc1->GetEntries()<<" "<<hmc2->GetEntries()<<" "<<hmc3->GetEntries()<<endl;
712 debug<<((fitsettings&0x1)&&hmc2->GetEntries()<=minentries)<<" "<<((fitsettings&0x2)&&hmc3->GetEntries()<=minentries)<<endl;
57a75b31 713 cout<<hToFit->GetEntries()<<" "<<hmc1->GetEntries()<<" "<<hmc2->GetEntries()<<" "<<hmc3->GetEntries()<<endl;
714 cout<<((fitsettings&0x1)&&hmc2->GetEntries()<=minentries)<<" "<<((fitsettings&0x2)&&hmc3->GetEntries()<=minentries)<<endl;
2245e1bd 715 if(hToFit->GetEntries()<=minentries || hmc1->GetEntries()<=minentries || ((fitsettings&0x1)&&hmc2->GetEntries()<=minentries) || ((fitsettings&0x2)&&hmc3->GetEntries()<=minentries))
88d02ae4 716 continue;
88d02ae4 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};
2245e1bd 723 Int_t binCutRange[]={hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1])};
724
88d02ae4 725
88d02ae4 726 if(hmc3->GetNbinsX()>300)
727 {
728
2245e1bd 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]);
88d02ae4 733
734 hToFit->Rebin(30);
735 hmc1->Rebin(30);
736 hmc2->Rebin(30);
737 hmc3->Rebin(30);
2245e1bd 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]);
88d02ae4 747 else
748 corrforrebinning[0]=1.0;
2245e1bd 749 if(hmc1->Integral(binCutRange[0],binCutRange[1])>0.0)
750 corrforrebinning[1]=corrforrebinning[1]/hmc1->Integral(binCutRange[0],binCutRange[1]);
88d02ae4 751 else
752 corrforrebinning[1]=1.0;
2245e1bd 753 if(hmc2->Integral(binCutRange[0],binCutRange[1])>0.0)
754 corrforrebinning[2]=corrforrebinning[2]/hmc2->Integral(binCutRange[0],binCutRange[1]);
88d02ae4 755 else
756 corrforrebinning[2]=1.0;
2245e1bd 757 if(hmc3->Integral(binCutRange[0],binCutRange[1])>0.0)
758 corrforrebinning[3]=corrforrebinning[3]/hmc3->Integral(binCutRange[0],binCutRange[1]);
88d02ae4 759 else
760 corrforrebinning[3]=1.0;
761
88d02ae4 762
763
764 }
765
766 cDCA->cd();
767 gPad->SetGridy();
768 gPad->SetGridx();
769 gPad->SetLogy();
770
771 TObjArray *mc=0x0;
1af4558e 772 Int_t Npar=3;
2245e1bd 773 if(fitsettings==3)
88d02ae4 774 mc = new TObjArray(3); // MC histograms are put in this array
775 else
1af4558e 776 {
88d02ae4 777 mc = new TObjArray(2);
1af4558e 778 Npar=2;
779 }
2245e1bd 780
88d02ae4 781 mc->Add(hmc1);
2245e1bd 782 if(fitsettings&0x1)
783 mc->Add(hmc2);
784 if(fitsettings&0x2)
88d02ae4 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
1af4558e 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 }
2245e1bd 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]);
c0a5ac55 808 // hToFit->SetTitle(Form("DCA distr - %s %s %s %lf",sample[isample].Data(),Particle[ipart].Data(),Sign[icharge].Data(),lowedge));
809 hToFit->SetTitle("");
88d02ae4 810 Int_t status = fit->Fit(); // perform the fit
811 cout << "fit status: " << status << endl;
812 debug<<"fit status: " << status << endl;
813
2245e1bd 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
88d02ae4 821 TH1F* result = (TH1F*) fit->GetPlot();
822 TH1F* PrimMCPred=(TH1F*)fit->GetMCPrediction(0);
2245e1bd 823
824 TH1F* secStMCPred=0X0;
88d02ae4 825
826 TH1F* secMCPred=0x0;
88d02ae4 827
88d02ae4 828 //first method, use directly the fit result
829 fit->GetResult(0,v1,ev1);
2245e1bd 830
831 if(fitsettings&0x1)
832 {
833 fit->GetResult(1,v2,ev2);
834 secStMCPred=(TH1F*)fit->GetMCPrediction(1);
835
836 }
837 if(fitsettings&0x2)
88d02ae4 838 {
2245e1bd 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
88d02ae4 845 }
2245e1bd 846
88d02ae4 847
2245e1bd 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
88d02ae4 854
2245e1bd 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
88d02ae4 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;
8cec8e63 875 debug<<"addtional info"<<endl;
8cec8e63 876
877
2245e1bd 878 //Method 2 output histo
8cec8e63 879
2245e1bd 880 Float_t normalizationdata1=result->Integral(binCutRange[0],binCutRange[1])/result->Integral(binFitRange[0],binFitRange[1]);
88d02ae4 881
882
a4ca64e3 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
88d02ae4 888 normalizationdata1*=corrforrebinning[0];
889
890
2245e1bd 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;
88d02ae4 898
899 normalizationmc11*=corrforrebinning[1];
900 normalizationmc21*=corrforrebinning[2];
901 normalizationmc31*=corrforrebinning[3];
902
2245e1bd 903 debug<<"After Nor 2"<<endl;
88d02ae4 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
88d02ae4 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);
1e6bff53 914 hprimary[index+6*isample]->SetBinError(ibin_data,ev1*normalizationmc11);
2245e1bd 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));
88d02ae4 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])));
2245e1bd 924
88d02ae4 925 hToFit->SetMinimum(0.0001);
c0a5ac55 926 hToFit->GetXaxis()->SetTitle("DCA_{xy} (cm)");
927 hToFit->GetYaxis()->SetTitle("N_{counts}/N_{counts}(-3cm;3cm)");
928 hToFit->GetYaxis()->SetTitleOffset(1.3);
88d02ae4 929 hToFit->DrawClone("E1x0");
c0a5ac55 930 Leg1->AddEntry(hToFit,"data","p");
88d02ae4 931 result->SetTitle("Fit result");
2245e1bd 932 result->SetLineColor(kBlack);
c0a5ac55 933 Leg1->AddEntry(result,"fit result","l");
934 result->DrawClone("histsame");
2245e1bd 935
936 PrimMCPred->SetLineColor(kGreen+2);
b6a74f8a 937 PrimMCPred->SetLineStyle(2);
938 PrimMCPred->SetLineWidth(3.0);
c0a5ac55 939 Leg1->AddEntry(PrimMCPred,"primaries","l");
940 PrimMCPred->DrawClone("histsame");
2245e1bd 941 if(fitsettings&0x1)
88d02ae4 942 {
2245e1bd 943
944 secStMCPred->Scale(v2/secStMCPred->Integral(secStMCPred->GetXaxis()->FindBin(FitRange[0]),secStMCPred->GetXaxis()->FindBin(FitRange[1])));
945 secStMCPred->SetLineColor(kRed);
b6a74f8a 946 secStMCPred->SetLineWidth(3.0);
947
948 secStMCPred->SetLineStyle(3);
c0a5ac55 949 Leg1->AddEntry(secStMCPred,"weak decays","l");
950 secStMCPred->DrawClone("histsame");
2245e1bd 951
952 }
953 if(fitsettings&0x2)
954 {
955
956 secMCPred->Scale(v3/secMCPred->Integral(secMCPred->GetXaxis()->FindBin(FitRange[0]),secMCPred->GetXaxis()->FindBin(FitRange[1])));
b6a74f8a 957 secMCPred->SetLineColor(kBlue);
958 secMCPred->SetLineWidth(3.0);
959
960 secMCPred->SetLineStyle(4);
c0a5ac55 961 Leg1->AddEntry(secMCPred,"material","l");
962 secMCPred->DrawClone("histsame");
2245e1bd 963
964 }
965 }
88d02ae4 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 }
c0a5ac55 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));
88d02ae4 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
997void RecomputeErrors(TH1* h)
998{
999 for (int i=0; i<=h->GetXaxis()->GetNbins(); i++)
50c532c4 1000 {
1001 cout<<h->GetBinContent(i)<<" "<<h->GetBinError(i)<<" error "<<TMath::Sqrt(h->GetBinContent(i))<<endl;
88d02ae4 1002 h->SetBinError(i,TMath::Sqrt(h->GetBinContent(i)));
50c532c4 1003 }
88d02ae4 1004 h->Sumw2();
1005}
1006void 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
1016void 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}
1031void 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
b3ea73e1 1042void GFCorrection(TH1F **Spectra,Float_t tofpt,UInt_t options)
88d02ae4 1043{
b3ea73e1 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";
57a75b31 1056 TFile *fGeanFlukaK= TFile::Open(fnameGeanFlukaK.Data());
2245e1bd 1057 if (!fGeanFlukaK)
1058 {
ab0298c1 1059 fnameGeanFlukaK="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/correctionForCrossSection.321.root";
57a75b31 1060 fGeanFlukaK= TFile::Open(fnameGeanFlukaK.Data());
2245e1bd 1061 if (!fGeanFlukaK)
1062 return;
1063 }
b3ea73e1 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);
bdbba618 1075 if(tcutsdata->GetUseTypeDependedTOFCut())
1076 tofpt=tcutsdata->GetPtTOFMatchingKaon();
1077
b3ea73e1 1078 for(Int_t binK=0;binK<=Spectra[1]->GetNbinsX();binK++)
1079 {
88d02ae4 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 }
b3ea73e1 1111 }
88d02ae4 1112 }
b3ea73e1 1113 if(!(options&kgeantflukaProton))
1114 return;
88d02ae4 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;
50c532c4 1120 TString fnameGFProtons= "GFCorrection/correctionForCrossSection.root";
57a75b31 1121 TFile* fGFProtons = TFile::Open(fnameGFProtons.Data());
2245e1bd 1122 if (!fGFProtons)
1123 {
57a75b31 1124 fnameGFProtons="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/correctionForCrossSection.root";
1125 fGFProtons = TFile::Open(fnameGFProtons.Data());
2245e1bd 1126 if (!fGFProtons)
1127 return;
1128 }
1129
1130
1131
88d02ae4 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
bdbba618 1147 Int_t nbins = Spectra[2]->GetNbinsX();
1148 if(tcutsdata->GetUseTypeDependedTOFCut())
1149 tofpt=tcutsdata->GetPtTOFMatchingProton();
1150
88d02ae4 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///////////
1205TF1 *
1206TrackingEff_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
1222Double_t
1223TrackingPtGeantFlukaCorrectionNull(Double_t pTmc)
1224{
1225 return 1.;
1226}
1227
1228Double_t
1229TrackingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
1230{
1231 return (1 - 0.129758 *TMath::Exp(-pTmc*0.679612));
1232}
1233
1234Double_t
1235TrackingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
1236{
1237 return TMath::Min((0.972865 + 0.0117093*pTmc), 1.);
1238}
1239///////////////////////////////////////////
1240TF1 *
1241TOFmatchMC_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
1258Double_t
1259MatchingPtGeantFlukaCorrectionNull(Double_t pTmc)
1260{
1261 return 1.;
1262}
1263
1264Double_t
1265MatchingPtGeantFlukaCorrectionPrMinus(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
1271Double_t
1272MatchingPtGeantFlukaCorrectionKaMinus(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}
1277void MatchingTOFEff(TH1F** Spectra, TList* list=0x0)
1278{
1279 if(TOFMatchingScalling[0]<0.0&&TOFMatchingScalling[1]<0.0)
1280 {
50c532c4 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");
88d02ae4 1285 hMatcEffPos_data->SetTitle("Matching Eff Pos - data");
1286 TH1F *hMatcEffNeg_data=(TH1F*)tcutsdata->GetHistoNMatchedNeg();
50c532c4 1287 hMatcEffNeg_data->Sumw2();
1288 //hMatcEffNeg_data->Divide((TH1F*)tcutsdata->GetHistoNSelectedNeg());
1289 hMatcEffNeg_data->Divide(hMatcEffNeg_data,(TH1F*)tcutsdata->GetHistoNSelectedNeg(),1,1,"B");
88d02ae4 1290 hMatcEffNeg_data->SetTitle("Matching Eff Neg - data");
1291 TH1F *hMatcEffPos_mc=(TH1F*)tcutsmc->GetHistoNMatchedPos();
50c532c4 1292 hMatcEffPos_mc->Sumw2();
1293 //hMatcEffPos_mc->Divide((TH1F*)tcutsmc->GetHistoNSelectedPos());
1294 hMatcEffPos_mc->Divide(hMatcEffPos_mc,(TH1F*)tcutsmc->GetHistoNSelectedPos(),1,1,"B");
88d02ae4 1295 hMatcEffPos_mc->SetTitle("Matching Eff Pos - mc");
1296 TH1F *hMatcEffNeg_mc=(TH1F*)tcutsmc->GetHistoNMatchedNeg();
50c532c4 1297 hMatcEffNeg_mc->Sumw2();
1298 //hMatcEffNeg_mc->Divide((TH1F*)tcutsmc->GetHistoNSelectedNeg());
1299 hMatcEffNeg_mc->Divide(hMatcEffNeg_mc,(TH1F*)tcutsmc->GetHistoNSelectedNeg(),1,1,"B");
88d02ae4 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);
50c532c4 1320
88d02ae4 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}
1339Double_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
ab0298c1 1345TH1* GetSumAllCh(TH1F** spectra, Double_t* mass,Double_t etacut)
88d02ae4 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);
ab0298c1 1361 Float_t maxy=eta2y(pt,masstmp,etacut/2.0);
88d02ae4 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))));
ab0298c1 1363 conver=conver/etacut;
a4ca64e3 1364 //cout<<maxy<<" "<<conver<<" "<<masstmp<<""<<spectra[i]->GetName()<<endl;
88d02ae4 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
1380void 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}
2245e1bd 1393
1394Short_t DCAfitsettings (Float_t pt, Int_t type)
1395{
2245e1bd 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
3a71f081 1403}
1404
1405Float_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);
a4ca64e3 1419 //cout<<histodata->GetBinContent(2)<<endl;
3a71f081 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);
1af4558e 1454 else if ((options&kknormalizationwithbin0integralsMC)==knormalizationwithbin0integralsMC)
3a71f081 1455 return dataevents+bin0mcMC->Integral(binmin,binmax) ;
1456 else
1457 return 1;
1458}
1459
a4ca64e3 1460
1461Bool_t ReadConfigFile(TString configfile)
1462{
1463 ifstream infile(configfile.Data());
1464 if(infile.is_open()==false)
1465 return false;
93531558 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"};
a4ca64e3 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();
1af4558e 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();
cd28560b 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();
93531558 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
a4ca64e3 1559 continue;
1560
1af4558e 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
a4ca64e3 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 }
1af4558e 1595 for(int i=0;i<3;i++)
1596 if(minRanges[i]>maxRanges[i])
1597 return false;
a4ca64e3 1598
1599
1600 return true;
1601}
50c532c4 1602
1603void 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
1af4558e 1615}
1616
1617void 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
cd28560b 1635}
1636void 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
c21424b2 1649}
1650void 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
93531558 1674}
1675
1676void 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();
cd28560b 1704}