]>
Commit | Line | Data |
---|---|---|
88d02ae4 | 1 | class AliSpectraBothHistoManager; |
2 | class AliSpectraBothEventCuts; | |
3 | class AliSpectraBothTrackCuts; | |
4 | TString Charge[]={"Pos","Neg"}; | |
5 | TString Sign[]={"Plus","Minus"}; | |
6 | TString Particle[]={"Pion","Kaon","Proton"}; | |
7 | AliSpectraBothHistoManager* managerdata=0x0; | |
8 | AliSpectraBothEventCuts* ecutsdata=0x0; | |
9 | AliSpectraBothTrackCuts* tcutsdata=0x0; | |
10 | ||
11 | AliSpectraBothHistoManager* managermc=0x0; | |
12 | AliSpectraBothEventCuts* ecutsmc=0x0; | |
13 | AliSpectraBothTrackCuts* tcutsmc=0x0; | |
14 | ||
15 | Float_t TOFMatchingScalling[2]={-1,-1}; | |
16 | Int_t Color[3]={1,2,4}; | |
17 | Int_t Marker[6]={20,21,22,24,25,26}; | |
18 | Double_t Range[3]={0.3,0.3,0.5}; // LowPt range for pi k p | |
19 | ||
20 | enum ECharge_t { | |
21 | kPositive, | |
22 | kNegative, | |
23 | kNCharges | |
24 | }; | |
b3ea73e1 | 25 | enum { |
26 | kdodca=0x1, | |
27 | kgeantflukaKaon=0x2, | |
28 | kgeantflukaProton=0x4, | |
29 | knormalizationtoeventspassingPhySel=0x8, | |
0c1898ba | 30 | kveretxcorrectionandbadchunkscorr=0x10, |
8cec8e63 | 31 | kmcisusedasdata=0x20, |
32 | kdonotusedcacuts=0x40 | |
b3ea73e1 | 33 | }; |
88d02ae4 | 34 | |
0c1898ba | 35 | Bool_t OpenFile(TString dirname, TString outputname, Bool_t mcflag,Bool_t mcasdata=false); |
b3ea73e1 | 36 | void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TString outnamemc="" ) |
88d02ae4 | 37 | { |
38 | TH1::AddDirectory(kFALSE); | |
39 | gSystem->Load("libCore.so"); | |
40 | gSystem->Load("libPhysics.so"); | |
41 | gSystem->Load("libTree"); | |
42 | gSystem->Load("libMatrix"); | |
43 | gSystem->Load("libSTEERBase"); | |
44 | gSystem->Load("libESD"); | |
45 | gSystem->Load("libAOD"); | |
46 | gSystem->Load("libANALYSIS"); | |
47 | gSystem->Load("libOADB"); | |
48 | gSystem->Load("libANALYSISalice"); | |
49 | gSystem->Load("libTENDER"); | |
50 | gSystem->Load("libCORRFW"); | |
51 | gSystem->Load("libPWGTools"); | |
52 | gSystem->Load("libPWGLFspectra"); | |
53 | ||
54 | gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/QAPlotsBoth.C"); | |
55 | Double_t mass[3]; | |
56 | mass[0] = TDatabasePDG::Instance()->GetParticle("pi+")->Mass(); | |
57 | mass[1] = TDatabasePDG::Instance()->GetParticle("K+")->Mass(); | |
58 | mass[2] = TDatabasePDG::Instance()->GetParticle("proton")->Mass(); | |
59 | ||
8cec8e63 | 60 | TFormula* dcacutxy=0x0; |
61 | if(!(options&kdonotusedcacuts)) | |
62 | { | |
63 | ||
64 | AliESDtrackCuts* esdtrackcuts= AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE); | |
65 | TString formulastring(esdtrackcuts->GetMaxDCAToVertexXYPtDep()); | |
66 | formulastring.ReplaceAll("pt","x"); | |
67 | dcacutxy=new TFormula("dcacutxy",formulastring.Data()); | |
68 | } | |
88d02ae4 | 69 | TList* lout=new TList(); |
70 | ||
71 | ||
72 | TString indirname=Form("/output/train%s",outdate.Data()); | |
73 | //TString indirname("/output/train24072012"); | |
74 | if(outnamemc.Length()==0) | |
75 | outnamemc=outnamedata; | |
76 | cout<<indirname.Data()<<" "<<outnamemc.Data()<<endl; | |
77 | // Do the job | |
78 | ||
79 | ||
80 | OpenFile(indirname,outnamemc,true); | |
0c1898ba | 81 | OpenFile(indirname,outnamedata,false,((Bool_t)(options&kmcisusedasdata))); |
88d02ae4 | 82 | if(!managermc||!managerdata) |
83 | { | |
84 | cout<<managermc<<" "<<managerdata<<endl; | |
85 | return; | |
86 | } | |
87 | TH1F* rawspectradata[6]; | |
88 | TH1F* rawspectramc[6]; | |
89 | TH1F* MCTruth[6]; | |
90 | TH1F* eff[6]; | |
91 | TH1F* contallMC[6]; | |
92 | TH1F* contPID[6]; | |
93 | TH1F* contWD[6]; | |
94 | TH1F* contMat[6]; | |
95 | TH1F* confinal[6]; | |
96 | ||
97 | TH1F* contfit[12]; | |
98 | TH1F* contWDfit[12]; | |
99 | TH1F* contMatfit[12]; | |
100 | TH1F* primaryfit[12]; | |
101 | ||
102 | ||
103 | ||
104 | TH1F* spectra[6]; | |
105 | TH1F* spectraLeonardo[6]; | |
106 | ||
107 | TH1F* corrLeonardo[6]; | |
108 | //GetSpectra(managerdata,rawspectradata,true); | |
109 | //GetSpectra(managermc,rawspectramc,true,true); | |
110 | ||
111 | GetPtHistFromPtDCAhisto("hHistPtRecSigma","SpectraMC",managermc,rawspectramc,dcacutxy); | |
112 | GetPtHistFromPtDCAhisto("hHistPtRecSigma","SpectraDATA",managerdata,rawspectradata,dcacutxy); | |
113 | GetPtHistFromPtDCAhisto("hHistPtRecTruePrimary","eff",managermc,eff,dcacutxy); | |
114 | GetPtHistFromPtDCAhisto("hHistPtRecTrue","conPID",managermc,contPID,dcacutxy); | |
115 | GetPtHistFromPtDCAhisto("hHistPtRecSigmaSecondaryWeakDecay","conWD",managermc,contWD,dcacutxy); | |
116 | GetPtHistFromPtDCAhisto("hHistPtRecSigmaSecondaryMaterial","conMat",managermc,contMat,dcacutxy); | |
117 | ||
118 | ||
119 | // Double_t neventsdata = ecutsdata->NumberOfEvents(); | |
b3ea73e1 | 120 | Double_t neventsmcall = 1 ; //if loop over MC is done after or befor events cuts this will be changed |
121 | Double_t neventsdata = 1; | |
122 | Double_t neventsmc = 1; | |
123 | ||
124 | if(options&knormalizationtoeventspassingPhySel) | |
125 | { | |
126 | neventsmcall= ecutsmc->NumberOfProcessedEvents(); | |
127 | neventsdata=ecutsdata->NumberOfPhysSelEvents(); | |
128 | neventsmc=ecutsmc->NumberOfPhysSelEvents(); | |
129 | } | |
130 | else | |
131 | { | |
132 | neventsdata=ecutsdata->NumberOfEvents(); | |
133 | neventsmc=ecutsmc->NumberOfEvents(); | |
134 | neventsmcall= ecutsmc->NumberOfEvents(); | |
135 | ||
136 | ||
137 | } | |
88d02ae4 | 138 | GetMCTruth(MCTruth); |
139 | ||
140 | ||
141 | ||
142 | TH1F* allgen=((TH1F*)managermc->GetPtHistogram1D("hHistPtGen",1,1))->Clone(); | |
143 | allgen->SetName("AllGen"); | |
144 | TH1F* allrecMC=GetOneHistFromPtDCAhisto("hHistPtRec","rawallMC",managermc,dcacutxy); | |
145 | TH1F* alleff=GetOneHistFromPtDCAhisto("hHistPtRecPrimary","effall",managermc,dcacutxy); | |
146 | TH1F* allrecdata=GetOneHistFromPtDCAhisto("hHistPtRec","rawalldata",managerdata,dcacutxy); | |
147 | ||
148 | ||
149 | ||
150 | ||
151 | TH1F* spectraall=(TH1F*)allrecdata->Clone("recNch"); | |
152 | TH1F* contall=(TH1F*)allrecMC->Clone("contall"); | |
153 | contall->Add(alleff,-1); | |
154 | alleff->Divide(alleff,allgen,1,1,"B"); | |
b3ea73e1 | 155 | contall->Divide(contall,allrecMC,1,1,"B"); |
88d02ae4 | 156 | |
157 | GetCorrectedSpectra(spectraall,allrecdata,alleff,contall); | |
158 | Divideby2pipt(spectraall); | |
159 | ||
160 | allrecdata->Scale(1./neventsdata,"width"); | |
161 | allgen->Scale(1./neventsmcall,"width"); | |
162 | allrecMC->Scale(1./neventsmc,"width"); | |
163 | spectraall->Scale(1./neventsdata,"width"); | |
164 | ||
165 | ||
166 | lout->Add(allgen); | |
167 | lout->Add(allrecMC); | |
168 | lout->Add(alleff); | |
169 | lout->Add(allrecdata); | |
170 | lout->Add(spectraall); | |
1e6bff53 | 171 | lout->Add(contall); |
172 | ||
88d02ae4 | 173 | for (int i=0;i<6;i++) |
174 | { | |
175 | ||
176 | ||
177 | TString tmpname(rawspectramc[i]->GetTitle()); | |
178 | tmpname.ReplaceAll("SpectraMC","%s"); | |
179 | contallMC[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contallMC")); | |
180 | contfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contfit")); | |
181 | contWDfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contWDfit")); | |
182 | contMatfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contMatfit")); | |
183 | primaryfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"primaryfit")); | |
184 | ||
185 | contfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contfitonMC")); | |
186 | contWDfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contWDfitonMC")); | |
187 | contMatfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contMatfitonMC")); | |
188 | primaryfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"primaryfitMC")); | |
189 | ||
190 | contfit[i]->Reset(); | |
191 | contWDfit[i]->Reset(); | |
192 | contMatfit[i]->Reset(); | |
193 | primaryfit[i]->Reset(); | |
194 | ||
195 | ||
196 | contfit[i+6]->Reset(); | |
197 | contWDfit[i+6]->Reset(); | |
198 | contMatfit[i+6]->Reset(); | |
199 | primaryfit[i+6]->Reset(); | |
200 | ||
201 | SetBintoOne(primaryfit[i]); | |
202 | SetBintoOne(primaryfit[i+6]); | |
203 | spectra[i]=(TH1F*)rawspectradata[i]->Clone(Form(tmpname.Data(),"SpectraFinal")); | |
204 | ||
205 | spectraLeonardo[i]=(TH1F*)rawspectradata[i]->Clone(Form(tmpname.Data(),"SpectraFinalLeonardo")); | |
206 | corrLeonardo[i]=(TH1F*)MCTruth[i]->Clone(Form(tmpname.Data(),"CorrFactLeonardo")); | |
207 | ||
208 | corrLeonardo[i]->Divide(corrLeonardo[i],rawspectramc[i],1,1,"B"); | |
209 | ||
210 | ||
211 | ||
212 | contallMC[i]->Add(eff[i],-1.0); | |
213 | RecomputeErrors(contallMC[i]); | |
214 | contallMC[i]->Sumw2(); | |
215 | contallMC[i]->Divide(contallMC[i],rawspectramc[i],1,1,"B"); | |
216 | ||
217 | eff[i]->Divide(eff[i],MCTruth[i],1,1,"B"); | |
218 | ||
219 | ||
220 | contPID[i]->Sumw2(); | |
221 | rawspectramc[i]->Sumw2(); | |
222 | contPID[i]->Add(contPID[i],rawspectramc[i],-1,1); | |
223 | RecomputeErrors(contPID[i]); | |
224 | contPID[i]->ResetStats(); | |
225 | contPID[i]->Sumw2(); | |
226 | contPID[i]->Divide(contPID[i],rawspectramc[i],1,1,"B"); | |
227 | ||
228 | confinal[i]=(TH1F*)contPID[i]->Clone(Form(tmpname.Data(),"confinal")); | |
229 | ||
230 | ||
231 | contWD[i]->Divide(contWD[i],rawspectramc[i],1,1,"B"); | |
232 | contMat[i]->Divide(contMat[i],rawspectramc[i],1,1,"B"); | |
233 | ||
234 | ||
235 | ||
236 | rawspectradata[i]->Scale(1./neventsdata,"width"); | |
237 | rawspectramc[i]->Scale(1./neventsmc,"width"); | |
238 | MCTruth[i]->Scale(1./neventsmcall,"width"); | |
239 | spectraLeonardo[i]->Scale(1./neventsdata,"width"); | |
240 | ||
241 | ||
242 | ||
243 | lout->Add(rawspectradata[i]); | |
244 | lout->Add(rawspectramc[i]); | |
245 | lout->Add(MCTruth[i]); | |
246 | lout->Add(eff[i]); | |
247 | lout->Add(contallMC[i]); | |
248 | lout->Add(contPID[i]); | |
249 | lout->Add(contWD[i]); | |
250 | lout->Add(contMat[i]); | |
251 | lout->Add(contfit[i]); | |
252 | lout->Add(contWDfit[i]); | |
253 | lout->Add(contMatfit[i]); | |
1e6bff53 | 254 | lout->Add(primaryfit[i]); |
88d02ae4 | 255 | lout->Add(contfit[i+6]); |
256 | lout->Add(contWDfit[i+6]); | |
257 | lout->Add(contMatfit[i+6]); | |
1e6bff53 | 258 | lout->Add(primaryfit[i+6]); |
88d02ae4 | 259 | lout->Add(spectra[i]); |
260 | lout->Add(spectraLeonardo[i]); | |
261 | lout->Add(confinal[i]); | |
262 | } | |
263 | TFile* fout=new TFile(Form("./results/ResMY_%s_%s.root",outnamemc.Data(),outdate.Data()),"RECREATE"); | |
b3ea73e1 | 264 | if (options&kdodca) |
88d02ae4 | 265 | DCACorrectionMarek(managerdata,managermc,dcacutxy,fout,contfit,contWDfit,contMatfit,primaryfit); |
266 | for (int i=0;i<6;i++) | |
267 | { | |
b3ea73e1 | 268 | if(options&kdodca) |
88d02ae4 | 269 | { |
270 | confinal[i]->Add(contfit[i]); | |
271 | GetCorrectedSpectra(spectra[i],rawspectradata[i],eff[i],confinal[i]); | |
272 | } | |
273 | else | |
274 | { | |
275 | GetCorrectedSpectra(spectra[i],rawspectradata[i],eff[i],contallMC[i]); | |
276 | } | |
277 | GetCorrectedSpectraLeonardo(spectraLeonardo[i],corrLeonardo[i],primaryfit[i],primaryfit[i+6]); | |
278 | CleanHisto(spectra[i],-1,100,contPID[i]); | |
279 | CleanHisto(spectraLeonardo[i],-1,100,contPID[i]); | |
280 | } | |
281 | ||
b3ea73e1 | 282 | GFCorrection(spectra,tcutsdata->GetPtTOFMatching(),options); |
283 | GFCorrection(spectraLeonardo,tcutsdata->GetPtTOFMatching(),options); | |
88d02ae4 | 284 | MatchingTOFEff(spectra,lout); |
285 | MatchingTOFEff(spectraLeonardo); | |
286 | TH1F* allch=GetSumAllCh(spectra,mass); | |
287 | lout->Add(allch); | |
288 | ||
289 | // lout->ls(); | |
290 | fout->cd(); | |
291 | TList* listqa=new TList(); | |
33d374d5 | 292 | TList* canvaslist=new TList(); |
b3ea73e1 | 293 | Float_t vertexcorrection=1.0; |
294 | Float_t corrbadchunksvtx=QAPlotsBoth(managerdata,managermc,ecutsdata,ecutsmc,tcutsdata,tcutsmc,listqa,canvaslist); | |
295 | if (options&kveretxcorrectionandbadchunkscorr) | |
296 | vertexcorrection=corrbadchunksvtx; | |
88d02ae4 | 297 | cout<<" VTX corr="<<vertexcorrection<<endl; |
b3ea73e1 | 298 | Double_t ycut=tcutsdata->GetY(); |
299 | if(TMath::Abs(ycut)>0.0) | |
300 | vertexcorrection=vertexcorrection/(2.0*ycut); | |
88d02ae4 | 301 | for (int i=0;i<6;i++) |
302 | { | |
303 | spectra[i]->Scale(vertexcorrection); | |
304 | spectraLeonardo[i]->Scale(vertexcorrection); | |
890fd855 | 305 | if(TMath::Abs(ycut)>0.0) |
306 | { | |
307 | rawspectradata[i]->Scale(1.0/(2.0*ycut)); | |
308 | rawspectramc[i]->Scale(1.0/(2.0*ycut)); | |
309 | MCTruth[i]->Scale(1.0/(2.0*ycut)); | |
310 | } | |
88d02ae4 | 311 | } |
312 | allch->Scale(vertexcorrection); | |
313 | spectraall->Scale(vertexcorrection/1.6); | |
314 | ||
315 | //spectraall->Scale(1.0/1.6); | |
316 | lout->Write("output",TObject::kSingleKey); | |
317 | listqa->Write("outputQA",TObject::kSingleKey); | |
33d374d5 | 318 | canvaslist->Write("outputcanvas",TObject::kSingleKey); |
319 | ||
88d02ae4 | 320 | fout->Close(); |
321 | ||
322 | } | |
323 | ||
0c1898ba | 324 | Bool_t OpenFile(TString dirname,TString outputname, Bool_t mcflag, Bool_t mcasdata) |
88d02ae4 | 325 | { |
0c1898ba | 326 | |
88d02ae4 | 327 | |
328 | TString nameFile = Form("./%s/AnalysisResults%s.root",dirname.Data(),(mcflag?"MC":"DATA")); | |
329 | TFile *file = TFile::Open(nameFile.Data()); | |
330 | if(!file) | |
331 | { | |
332 | cout<<"no file"<<endl; | |
333 | return false; | |
334 | } | |
335 | TString sname=Form("OutputBothSpectraTask_%s_%s",(mcflag?"MC":"Data"),outputname.Data()); | |
0c1898ba | 336 | if(mcasdata) |
337 | { | |
338 | cout<<"using MC as data "<<endl; | |
339 | sname=Form("OutputBothSpectraTask_%s_%s","MC",outputname.Data()); | |
340 | } | |
88d02ae4 | 341 | file->ls(); |
342 | TDirectoryFile *dir=(TDirectoryFile*)file->Get(sname.Data()); | |
343 | if(!dir) | |
344 | { | |
ab13cf18 | 345 | // cout<<"no dir "<<sname.Data()<<endl; |
346 | if(mcasdata) | |
0c1898ba | 347 | { |
348 | cout<<"using MC as data "<<endl; | |
349 | sname=Form("OutputAODSpectraTask_%s_%s","MC",outputname.Data()); | |
350 | } | |
351 | else | |
352 | sname=Form("OutputAODSpectraTask_%s_%s",(mcflag?"MC":"Data"),outputname.Data()); | |
88d02ae4 | 353 | // cout<<"trying "<<sname.Data()<<endl; |
354 | dir=(TDirectoryFile*)file->Get(sname.Data()); | |
355 | if(!dir) | |
356 | { | |
357 | cout<<"no dir "<<sname.Data()<<endl; | |
358 | return false; | |
359 | } | |
360 | } | |
361 | cout << " -- Info about " <<(mcflag?"MC":"DATA") <<" -- "<< endl; | |
362 | if(mcflag) | |
363 | { | |
364 | managermc= (AliSpectraBothHistoManager*) dir->Get("SpectraHistos"); | |
365 | ecutsmc = (AliSpectraBothEventCuts*) dir->Get("Event Cuts"); | |
366 | tcutsmc = (AliSpectraBothTrackCuts*) dir->Get("Track Cuts"); | |
367 | ecutsmc->PrintCuts(); | |
368 | tcutsmc->PrintCuts(); | |
369 | if(!managermc||!ecutsmc||!tcutsmc) | |
370 | return false; | |
371 | } | |
372 | else | |
373 | { | |
374 | managerdata= (AliSpectraBothHistoManager*) dir->Get("SpectraHistos"); | |
375 | ecutsdata = (AliSpectraBothEventCuts*) dir->Get("Event Cuts"); | |
376 | tcutsdata = (AliSpectraBothTrackCuts*) dir->Get("Track Cuts"); | |
377 | ecutsdata->PrintCuts(); | |
378 | tcutsdata->PrintCuts(); | |
379 | if(!managerdata||!ecutsdata||!tcutsdata) | |
380 | return false; | |
381 | } | |
382 | return true; | |
383 | } | |
384 | ||
385 | void GetMCTruth(TH1F** MCTruth) | |
386 | { | |
387 | for(Int_t icharge=0;icharge<2;icharge++) | |
388 | { | |
389 | for(Int_t ipart=0;ipart<3;ipart++) | |
390 | { | |
391 | Int_t index=ipart+3*icharge; | |
392 | TString hname=Form("hHistPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()); | |
393 | MCTruth[index]=(TH1F*)((TH1F*)managermc->GetPtHistogram1D(hname.Data(),1,1))->Clone(); | |
394 | MCTruth[index]->SetName(Form("MCTruth_%s%s",Particle[ipart].Data(),Sign[icharge].Data())); | |
395 | MCTruth[index]->SetTitle(Form("MCTruth_%s%s",Particle[ipart].Data(),Sign[icharge].Data())); | |
396 | MCTruth[index]->Sumw2(); | |
397 | } | |
398 | } | |
399 | } | |
400 | ||
401 | TH1F* GetOneHistFromPtDCAhisto(TString name,TString hnameout,AliSpectraBothHistoManager* hman,TFormula* dcacutxy) | |
402 | { | |
403 | histo =(TH1F*)((TH1F*) hman->GetPtHistogram1D(name.Data(),-1,-1))->Clone(); | |
404 | histo->SetName(hnameout.Data()); | |
405 | histo->SetTitle(hnameout.Data()); | |
406 | ||
407 | if(dcacutxy) | |
408 | { | |
409 | for(int ibin=1;ibin<histo->GetNbinsX();ibin++) | |
410 | { | |
411 | Double_t lowedge=histo->GetBinLowEdge(ibin); | |
412 | Float_t cut=dcacutxy->Eval(lowedge); | |
413 | TH1F* dcahist=(TH1F*)hman->GetDCAHistogram1D(name.Data(),lowedge,lowedge)); | |
414 | Float_t inyield=dcahist->Integral(dcahist->GetXaxis()->FindBin(-1.0*cut),dcahist->GetXaxis()->FindBin(cut)); | |
415 | cout<<"corr data "<<histo->GetBinContent(ibin)<<" "<<inyield<<" "<<dcahist->Integral()<<" "<<hnameout.Data()<<endl; | |
416 | 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; | |
417 | histo->SetBinContent(ibin,inyield); | |
418 | histo->SetBinError(ibin,TMath::Sqrt(inyield)); | |
419 | } | |
420 | } | |
421 | histo->Sumw2(); | |
422 | return histo; | |
423 | } | |
424 | ||
425 | ||
426 | ||
427 | ||
428 | void GetPtHistFromPtDCAhisto(TString hnamein, TString hnameout, AliSpectraBothHistoManager* hman,TH1F** histo,TFormula* dcacutxy) | |
429 | { | |
430 | Float_t min[3]={0.3,0.3,0.4}; | |
431 | Float_t max[3]={1.5,1.2,2.2}; | |
432 | for(Int_t icharge=0;icharge<2;icharge++) | |
433 | { | |
434 | for(Int_t ipart=0;ipart<3;ipart++) | |
435 | { | |
436 | Int_t index=ipart+3*icharge; | |
437 | //TString hname=Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()); | |
438 | printf("Getting %s",hnamein.Data()); | |
439 | TString nameinfinal=Form("%s%s%s",hnamein.Data(),Particle[ipart].Data(),Sign[icharge].Data()); | |
440 | TString nameoutfinal=Form("%s%s%s",hnameout.Data(),Particle[ipart].Data(),Sign[icharge].Data()); | |
441 | ||
442 | ||
443 | histo[index]=GetOneHistFromPtDCAhisto(nameinfinal,nameoutfinal,hman,dcacutxy); | |
444 | /*histo[index] =(TH1F*)((TH1F*) hman->GetPtHistogram1D(Form("%s%s%s",hnamein.Data(),Particle[ipart].Data(),Sign[icharge].Data()),-1,-1))->Clone(); | |
445 | histo[index]->SetName(Form("%s%s%s",hnameout.Data(),Particle[ipart].Data(),Sign[icharge].Data())); | |
446 | histo[index]->SetTitle(Form("%s%s%s",hnameout.Data(),Particle[ipart].Data(),Sign[icharge].Data())); | |
447 | ||
448 | if(dcacutxy) | |
449 | { | |
450 | for(int ibin=1;ibin<histo[index]->GetNbinsX();ibin++) | |
451 | { | |
452 | Double_t lowedge=histo[index]->GetBinLowEdge(ibin); | |
453 | Float_t cut=dcacutxy->Eval(lowedge); | |
454 | TH1F* dcahist=(TH1F*)hman->GetDCAHistogram1D(Form("%s%s%s",hnamein.Data(),Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge)); | |
455 | Float_t inyield=dcahist->Integral(dcahist->GetXaxis()->FindBin(-1.0*cut),dcahist->GetXaxis()->FindBin(cut)); | |
456 | cout<<"corr data "<<histo[index]->GetBinContent(ibin)<<" "<<inyield<<" "<<dcahist->Integral()<<endl; | |
457 | histo[index]->SetBinContent(ibin,inyield); | |
458 | histo[index]->SetBinError(ibin,TMath::Sqrt(inyield)); | |
459 | } | |
460 | }*/ | |
461 | CleanHisto(histo[index],min[ipart],max[ipart]); | |
462 | // histo[index]->Sumw2(); | |
463 | } | |
464 | } | |
465 | } | |
466 | void CleanHisto(TH1F* h, Float_t minV, Float_t maxV,TH1* contpid=0x0) | |
467 | { | |
468 | for (int i=0;i<=h->GetNbinsX();i++) | |
469 | { | |
470 | if(h->GetXaxis()->GetBinCenter(i)<minV||h->GetXaxis()->GetBinCenter(i)>maxV) | |
471 | { | |
472 | h->SetBinContent(i,0); | |
473 | h->SetBinError(i,0); | |
474 | } | |
475 | if(contpid) | |
476 | { | |
477 | if(contpid->GetBinContent(i)>0.2) | |
478 | { | |
479 | h->SetBinContent(i,0); | |
480 | h->SetBinError(i,0); | |
481 | } | |
482 | } | |
483 | } | |
484 | } | |
485 | ||
486 | ||
487 | void DCACorrectionMarek(AliSpectraBothHistoManager* hman_data, AliSpectraBothHistoManager* hman_mc,TFormula* fun,TFile *fout,TH1F** hcon,TH1F** hconWD,TH1F** hconMat,TH1F** hprimary) | |
488 | { | |
489 | printf("\n\n-> DCA Correction"); | |
b3ea73e1 | 490 | Double_t FitRange[2]={-2.0,2.0}; |
8cec8e63 | 491 | Double_t CutRange[2]={-3.0,3.0}; |
88d02ae4 | 492 | Double_t minptformaterial[6]={0.0,100.0,0.0,0.0,100.0,0.0}; |
493 | Double_t maxptformaterial[6]={0.0,-100.0,1.3,0.0,-100.0,0.0}; | |
494 | Bool_t usefit[3]={true,false,true}; | |
495 | Double_t maxbinforfit[6]={1.5,0,2.0,1.5,0,2.0}; | |
496 | Printf("\DCACorr"); | |
497 | // TH1F *hcorrection[2]; | |
498 | /* TCanvas *ccorrection=new TCanvas("DCAcorrection","DCAcorrection",700,500); | |
499 | TCanvas *cRatiocorrection=new TCanvas("DCARatiocorrection","DCARatiocorrection",700,500); | |
500 | cRatiocorrection->Divide(2,1); | |
501 | ccorrection->Divide(2,1);*/ | |
502 | TString sample[2]={"data","mc"}; | |
503 | ofstream debug("debugDCA.txt"); | |
504 | TList* listofdcafits=new TList(); | |
505 | for(Int_t icharge=0;icharge<2;icharge++) | |
506 | { | |
507 | for(Int_t ipart=0;ipart<3;ipart++) | |
508 | { | |
509 | Int_t index=ipart+3*icharge; | |
510 | for(Int_t isample=0;isample<2;isample++) | |
511 | { | |
512 | ||
513 | /*hcorrection[isample]=(TH1F*)Spectra[index]->Clone(); | |
514 | hcorrection[isample]->Reset("all");*/ | |
515 | for(Int_t ibin_data=7;ibin_data<40;ibin_data++) | |
516 | { | |
517 | ||
518 | Double_t lowedge=hcon[index]->GetBinLowEdge(ibin_data); | |
519 | Double_t binwidth=hcon[index]->GetBinWidth(ibin_data); | |
8cec8e63 | 520 | debug<<"NEW "<<Particle[ipart].Data()<<" "<<Sign[icharge].Data()<<" "<<lowedge<<endl; |
521 | if(fun) | |
522 | { | |
523 | CutRange[1]=fun->Eval(lowedge); | |
524 | CutRange[0]=-1.0*CutRange[1]; | |
525 | } | |
88d02ae4 | 526 | debug<<"cut "<<CutRange[1]<<" "<<CutRange[0]<<endl; |
527 | Bool_t useMaterial=kFALSE; | |
528 | cout<<"try to fit "<<lowedge<<" "<<maxbinforfit[index]<<endl; | |
529 | if(lowedge>maxbinforfit[index]) | |
530 | continue; | |
531 | if(lowedge>minptformaterial[index]&&lowedge<maxptformaterial[index]) | |
532 | useMaterial=kTRUE; | |
533 | ||
534 | 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); | |
535 | if(isample==0) | |
536 | TH1F *hToFit =(TH1F*) ((TH1F*)hman_data->GetDCAHistogram1D(Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone(); | |
537 | if(isample==1) | |
538 | TH1F *hToFit =(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone(); | |
539 | debug<<Particle[ipart].Data()<<" "<<Sign[icharge].Data()<<" "<<lowedge<<endl; | |
540 | TH1F *hmc1=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone(); | |
541 | TH1F *hmc2=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaSecondaryWeakDecay%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone(); | |
542 | TH1F *hmc3=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaSecondaryMaterial%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone(); | |
543 | Double_t minentries=1; | |
544 | debug<<" Entries "<<isample<<" "<<hToFit->GetEntries()<<" "<<hmc1->GetEntries()<<" "<<hmc2->GetEntries()<<" "<<hmc3->GetEntries()<<endl; | |
545 | debug<<"2 Entries "<<isample<<" "<<hToFit->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))<<" "<<hmc1->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))<<" "<<hmc2->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))<<" "<<hmc3->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))<<endl; | |
546 | debug<< CutRange[0]<<" "<<CutRange[1]<<" "<<lowedge<<endl; | |
547 | if(hToFit->GetEntries()<=minentries || hmc1->GetEntries()<=minentries || hmc2->GetEntries()<=minentries || hmc3->GetEntries()<=minentries) | |
548 | continue; | |
549 | //Data and MC can have different stat | |
550 | hToFit->Sumw2(); | |
551 | hmc1->Sumw2(); | |
552 | hmc2->Sumw2(); | |
553 | hmc3->Sumw2(); | |
554 | ||
555 | Float_t corrforrebinning[4]={1.0,1.0,1.0,1.0}; | |
556 | ||
557 | ||
558 | if(hmc3->GetNbinsX()>300) | |
559 | { | |
560 | ||
561 | corrforrebinning[0]=hToFit->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1])); | |
562 | corrforrebinning[1]=hmc1->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1])); | |
563 | corrforrebinning[2]=hmc2->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1])); | |
564 | corrforrebinning[3]=hmc3->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1])); | |
565 | ||
566 | hToFit->Rebin(30); | |
567 | hmc1->Rebin(30); | |
568 | hmc2->Rebin(30); | |
569 | hmc3->Rebin(30); | |
570 | if(hToFit->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))>0.0) | |
571 | corrforrebinning[0]=corrforrebinning[0]/hToFit->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1])); | |
572 | else | |
573 | corrforrebinning[0]=1.0; | |
574 | if(hmc1->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))>0.0) | |
575 | corrforrebinning[1]=corrforrebinning[1]/hmc1->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1])); | |
576 | else | |
577 | corrforrebinning[1]=1.0; | |
578 | if(hmc2->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))>0.0) | |
579 | corrforrebinning[2]=corrforrebinning[2]/hmc2->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1])); | |
580 | else | |
581 | corrforrebinning[2]=1.0; | |
582 | if(hmc3->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))>0.0) | |
583 | corrforrebinning[3]=corrforrebinning[3]/hmc3->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1])); | |
584 | else | |
585 | corrforrebinning[3]=1.0; | |
586 | ||
587 | debug<<" cor bin "<<corrforrebinning[0]<<" "<<corrforrebinning[1]<<" "<<corrforrebinning[2]<<" "<<corrforrebinning[3]<<endl; | |
588 | ||
589 | ||
590 | } | |
591 | ||
592 | cDCA->cd(); | |
593 | gPad->SetGridy(); | |
594 | gPad->SetGridx(); | |
595 | gPad->SetLogy(); | |
596 | ||
597 | TObjArray *mc=0x0; | |
598 | if(useMaterial) | |
599 | mc = new TObjArray(3); // MC histograms are put in this array | |
600 | else | |
601 | mc = new TObjArray(2); | |
602 | mc->Add(hmc1); | |
603 | mc->Add(hmc2); | |
604 | if(useMaterial) | |
605 | mc->Add(hmc3); | |
606 | TFractionFitter* fit = new TFractionFitter(hToFit,mc); // initialise | |
607 | fit->Constrain(0,0.0,1.0); // constrain fraction 1 to be between 0 and 1 | |
608 | fit->Constrain(1,0.0,1.0); // constrain fraction 1 to be between 0 and 1 | |
609 | if(useMaterial) | |
610 | fit->Constrain(2,0.0,1.0); // constrain fraction 1 to be between 0 and 1 | |
611 | fit->SetRangeX(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1])); | |
612 | hToFit->GetXaxis()->SetRange(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1])); | |
613 | hToFit->SetTitle(Form("DCA distr - %s %s %s %lf",sample[isample].Data(),Particle[ipart].Data(),Sign[icharge].Data(),lowedge)); | |
614 | Int_t status = fit->Fit(); // perform the fit | |
615 | cout << "fit status: " << status << endl; | |
616 | debug<<"fit status: " << status << endl; | |
617 | ||
618 | if (status == 0 && usefit[ipart]) | |
619 | { // check on fit status | |
620 | TH1F* result = (TH1F*) fit->GetPlot(); | |
621 | TH1F* PrimMCPred=(TH1F*)fit->GetMCPrediction(0); | |
622 | TH1F* secStMCPred=(TH1F*)fit->GetMCPrediction(1); | |
623 | ||
624 | TH1F* secMCPred=0x0; | |
625 | if(useMaterial) | |
626 | secMCPred=(TH1F*)fit->GetMCPrediction(2); | |
627 | ||
628 | Double_t v1=0,v2=0,v3=0; | |
629 | Double_t ev1=0,ev2=0,ev3=0; | |
630 | Double_t cov=0.0; | |
631 | //first method, use directly the fit result | |
632 | fit->GetResult(0,v1,ev1); | |
633 | fit->GetResult(1,v2,ev2); | |
634 | if(useMaterial) | |
635 | { | |
636 | fit->GetResult(2,v3,ev3); | |
637 | fit->GetFitter()->GetCovarianceMatrixElement(1,2); | |
638 | } | |
639 | debug<<v1<<" "<<ev1<<" "<<v2<<" "<<ev2<<" "<<v3<<" "<<ev3<<" "<<endl; | |
640 | ||
641 | ||
642 | ||
643 | Float_t normalizationdata=hToFit->Integral(hToFit->GetXaxis()->FindBin(CutRange[0]),hToFit->GetXaxis()->FindBin(CutRange[1]))/hToFit->Integral(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1])); | |
644 | ||
645 | Float_t normalizationmc1=hmc1->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))/hmc1->Integral(hmc1->GetXaxis()->FindBin(FitRange[0]),hmc1->GetXaxis()->FindBin(FitRange[1]))/normalizationdata; | |
646 | Float_t normalizationmc2=hmc2->Integral(hmc2->GetXaxis()->FindBin(CutRange[0]),hmc2->GetXaxis()->FindBin(CutRange[1]))/hmc2->Integral(hmc1->GetXaxis()->FindBin(FitRange[0]),hmc2->GetXaxis()->FindBin(FitRange[1]))/normalizationdata; | |
647 | Float_t normalizationmc3=hmc3->Integral(hmc3->GetXaxis()->FindBin(CutRange[0]),hmc3->GetXaxis()->FindBin(CutRange[1]))/hmc3->Integral(hmc1->GetXaxis()->FindBin(FitRange[0]),hmc3->GetXaxis()->FindBin(FitRange[1]))/normalizationdata; | |
648 | debug<<"After Nor"<<endl; | |
649 | debug<<v1*normalizationmc1<<" "<<ev1*normalizationmc1<<" "<<v2*normalizationmc2<<" "<<ev2*normalizationmc2<<" "<<v3*normalizationmc3<<" "<<ev3*normalizationmc3<<" "<<endl; | |
650 | 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 | 651 | debug<<"addtional info"<<endl; |
652 | Float_t normalizationmc1b=(hmc1->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))/hmc1->Integral(hmc1->GetXaxis()->FindBin(FitRange[0]),hmc1->GetXaxis()->FindBin(FitRange[1])))/normalizationdata; | |
653 | ||
654 | debug<<normalizationmc1<<" "<<normalizationmc1b<<endl; | |
655 | debug<<hmc1->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))/hmc1->Integral(hmc1->GetXaxis()->FindBin(FitRange[0]),hmc1->GetXaxis()->FindBin(FitRange[1]))<<" "<<secStMCPred->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))/secStMCPred->Integral(hmc1->GetXaxis()->FindBin(FitRange[0]),hmc1->GetXaxis()->FindBin(FitRange[1]))<<endl; | |
656 | //debug<<hmc1->Integral(hmc1->GetXaxis()->FindBin(FitRange[0]),hmc1->GetXaxis()->FindBin(FitRange[1]))<<" "<<secStMCPred->Integral(hmc1->GetXaxis()->FindBin(FitRange[0]),hmc1->GetXaxis()->FindBin(FitRange[1]))<<endl; | |
657 | debug<<hmc2->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))/hmc2->Integral(hmc1->GetXaxis()->FindBin(FitRange[0]),hmc1->GetXaxis()->FindBin(FitRange[1]))<<" "<<PrimMCPred->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))/PrimMCPred->Integral(hmc1->GetXaxis()->FindBin(FitRange[0]),hmc1->GetXaxis()->FindBin(FitRange[1]))<<endl; | |
658 | //debug<<hmc2->Integral(hmc1->GetXaxis()->FindBin(FitRange[0]),hmc1->GetXaxis()->FindBin(FitRange[1]))<<" "<<PrimMCPred->Integral(hmc1->GetXaxis()->FindBin(FitRange[0]),hmc1->GetXaxis()->FindBin(FitRange[1]))<<endl; | |
659 | if(PrimMCPred->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))>0.0) | |
660 | { | |
661 | if(useMaterial) | |
662 | { | |
663 | debug<<" ITSsa "<<endl; | |
664 | debug<<PrimMCPred->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))/(PrimMCPred->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))+secStMCPred->Integral(hmc1->GetXaxis()->FindBin(FitRange[0]),hmc1->GetXaxis()->FindBin(FitRange[1]))+secMCPred->Integral(hmc1->GetXaxis()->FindBin(FitRange[0]),hmc1->GetXaxis()->FindBin(FitRange[1])))<<endl; | |
665 | } | |
666 | else | |
667 | { | |
668 | debug<<" ITSsa "<<endl; | |
669 | debug<<PrimMCPred->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))/(PrimMCPred->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))+secStMCPred->Integral(hmc1->GetXaxis()->FindBin(FitRange[0]),hmc1->GetXaxis()->FindBin(FitRange[1])))<<endl; | |
670 | ||
671 | } | |
672 | } | |
673 | else | |
674 | { | |
675 | debug<<" ITSsa "<<endl; | |
676 | debug<<" NO "<<endl; | |
677 | } | |
678 | ||
679 | ||
680 | ||
88d02ae4 | 681 | Float_t normalizationdata1=result->Integral(result->GetXaxis()->FindBin(CutRange[0]),result->GetXaxis()->FindBin(CutRange[1]))/result->Integral(result->GetXaxis()->FindBin(FitRange[0]),result->GetXaxis()->FindBin(FitRange[1])); |
682 | ||
683 | ||
684 | normalizationdata1*=corrforrebinning[0]; | |
685 | ||
686 | ||
687 | Float_t normalizationmc11=PrimMCPred->Integral(PrimMCPred->GetXaxis()->FindBin(CutRange[0]),PrimMCPred->GetXaxis()->FindBin(CutRange[1]))/PrimMCPred->Integral(PrimMCPred->GetXaxis()->FindBin(FitRange[0]),PrimMCPred->GetXaxis()->FindBin(FitRange[1]))/normalizationdata1; | |
688 | Float_t normalizationmc21=secStMCPred->Integral(secStMCPred->GetXaxis()->FindBin(CutRange[0]),secStMCPred->GetXaxis()->FindBin(CutRange[1]))/secStMCPred->Integral(secStMCPred->GetXaxis()->FindBin(FitRange[0]),secStMCPred->GetXaxis()->FindBin(FitRange[1]))/normalizationdata1; | |
689 | Float_t normalizationmc31=0; | |
690 | if(useMaterial) | |
691 | normalizationmc31=secMCPred->Integral(secMCPred->GetXaxis()->FindBin(CutRange[0]),secMCPred->GetXaxis()->FindBin(CutRange[1]))/secMCPred->Integral(secMCPred->GetXaxis()->FindBin(FitRange[0]),secMCPred->GetXaxis()->FindBin(FitRange[1]))/normalizationdata1; | |
692 | ||
693 | normalizationmc11*=corrforrebinning[1]; | |
694 | normalizationmc21*=corrforrebinning[2]; | |
695 | normalizationmc31*=corrforrebinning[3]; | |
696 | ||
697 | debug<<"After Nor 2"<<endl; | |
698 | debug<<v1*normalizationmc11<<" "<<ev1*normalizationmc11<<" "<<v2*normalizationmc21<<" "<<ev2*normalizationmc21<<" "<<v3*normalizationmc31<<" "<<ev3*normalizationmc31<<endl; | |
699 | ||
700 | 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; | |
701 | ||
702 | debug<<CutRange[0]<<" "<<CutRange[1]<<endl; | |
703 | debug<<" Entries "<<isample<<" "<<hToFit->GetEntries()<<" "<<hmc1->GetEntries()<<" "<<hmc2->GetEntries()<<" "<<hmc3->GetEntries()<<endl; | |
704 | debug<<"2 Entries "<<isample<<" "<<hToFit->Integral(result->GetXaxis()->FindBin(CutRange[0]),result->GetXaxis()->FindBin(CutRange[1]))<<" "<<hmc1->Integral(result->GetXaxis()->FindBin(CutRange[0]),result->GetXaxis()->FindBin(CutRange[1]))<<" "<<hmc2->Integral(result->GetXaxis()->FindBin(CutRange[0]),result->GetXaxis()->FindBin(CutRange[1]))<<" "<<hmc3->Integral(result->GetXaxis()->FindBin(CutRange[0]),result->GetXaxis()->FindBin(CutRange[1]))<<endl; | |
705 | ||
706 | ||
707 | hconWD[index+6*isample]->SetBinContent(ibin_data,v2*normalizationmc21); | |
708 | hconWD[index+6*isample]->SetBinError(ibin_data,ev2*normalizationmc21); | |
709 | hconMat[index+6*isample]->SetBinContent(ibin_data,v3*normalizationmc31); | |
710 | hconMat[index+6*isample]->SetBinError(ibin_data,ev3*normalizationmc31); | |
711 | hprimary[index+6*isample]->SetBinContent(ibin_data,v1*normalizationmc11); | |
1e6bff53 | 712 | hprimary[index+6*isample]->SetBinError(ibin_data,ev1*normalizationmc11); |
88d02ae4 | 713 | if(useMaterial) |
714 | { | |
715 | hcon[index+6*isample]->SetBinContent(ibin_data,v2*normalizationmc21+v3*normalizationmc31); | |
716 | hcon[index+6*isample]->SetBinError(ibin_data,TMath::Sqrt(ev2*ev2*normalizationmc21*normalizationmc21+ev3*ev3*normalizationmc31*normalizationmc31+cov*normalizationmc31*normalizationmc21)); | |
717 | } | |
718 | else | |
719 | { | |
720 | hcon[index+6*isample]->SetBinContent(ibin_data,v2*normalizationmc21); | |
721 | hcon[index+6*isample]->SetBinError(ibin_data,ev2*normalizationmc21); | |
722 | } | |
723 | ||
724 | ||
725 | ||
726 | //Drawing section | |
727 | result->Scale(1.0/result->Integral(result->GetXaxis()->FindBin(FitRange[0]),result->GetXaxis()->FindBin(FitRange[1]))); | |
728 | hToFit->Scale(1.0/hToFit->Integral(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1]))); | |
729 | PrimMCPred->Scale(v1/PrimMCPred->Integral(PrimMCPred->GetXaxis()->FindBin(FitRange[0]),PrimMCPred->GetXaxis()->FindBin(FitRange[1]))); | |
730 | secStMCPred->Scale(v2/secStMCPred->Integral(secStMCPred->GetXaxis()->FindBin(FitRange[0]),secStMCPred->GetXaxis()->FindBin(FitRange[1]))); | |
731 | if(useMaterial) | |
732 | secMCPred->Scale(v3/secMCPred->Integral(secMCPred->GetXaxis()->FindBin(FitRange[0]),secMCPred->GetXaxis()->FindBin(FitRange[1]))); | |
733 | ||
734 | result->SetLineColor(kBlack); | |
735 | PrimMCPred->SetLineColor(kGreen+2); | |
736 | secStMCPred->SetLineColor(kRed); | |
737 | ||
738 | hToFit->SetMinimum(0.0001); | |
739 | hToFit->DrawClone("E1x0"); | |
740 | result->SetTitle("Fit result"); | |
741 | result->DrawClone("lhistsame"); | |
742 | PrimMCPred->DrawClone("lhistsame"); | |
743 | secStMCPred->DrawClone("lhistsame"); | |
744 | if(useMaterial) | |
745 | { | |
746 | secMCPred->SetLineColor(kBlue); | |
747 | secMCPred->DrawClone("lhistsame"); | |
748 | } | |
749 | } | |
750 | else | |
751 | { | |
752 | hconWD[index+6*isample]->SetBinContent(ibin_data,0.0); | |
753 | hconWD[index+6*isample]->SetBinError(ibin_data,0.0); | |
754 | hconMat[index+6*isample]->SetBinContent(ibin_data,0.0); | |
755 | hconMat[index+6*isample]->SetBinError(ibin_data,0.0); | |
756 | hcon[index+6*isample]->SetBinContent(ibin_data,0.0); | |
757 | hcon[index+6*isample]->SetBinError(ibin_data,0.0); | |
758 | hprimary[index+6*isample]->SetBinContent(ibin_data,1.0); | |
759 | hprimary[index+6*isample]->SetBinError(ibin_data,0.0); | |
760 | } | |
761 | listofdcafits->Add(cDCA); | |
762 | ||
763 | //cDCA->Write(); | |
764 | delete hToFit; | |
765 | } | |
766 | ||
767 | ||
768 | } | |
769 | ||
770 | } | |
771 | } | |
772 | fout->cd(); | |
773 | listofdcafits->Write("DCAfits",TObject::kSingleKey); | |
774 | } | |
775 | ||
776 | void RecomputeErrors(TH1* h) | |
777 | { | |
778 | for (int i=0; i<=h->GetXaxis()->GetNbins(); i++) | |
779 | h->SetBinError(i,TMath::Sqrt(h->GetBinContent(i))); | |
780 | h->Sumw2(); | |
781 | } | |
782 | void SetBintoOne(TH1* h) | |
783 | { | |
784 | for (int i=0;i<=h->GetXaxis()->GetNbins(); i++) | |
785 | { | |
786 | h->SetBinContent(i,1); | |
787 | h->SetBinError(i,0); | |
788 | } | |
789 | } | |
790 | ||
791 | ||
792 | void GetCorrectedSpectra(TH1F* corr,TH1F* raw,TH1F* eff, TH1F* con) | |
793 | { | |
794 | for (int i=0;i<=corr->GetXaxis()->GetNbins(); i++) | |
795 | { | |
796 | corr->SetBinContent(i,1); | |
797 | corr->SetBinError(i,0); | |
798 | } | |
799 | corr->Sumw2(); | |
800 | corr->Add(con,-1); | |
801 | corr->Sumw2(); | |
802 | corr->Divide(eff); | |
803 | corr->Sumw2(); | |
804 | corr->Multiply(raw); | |
805 | corr->Sumw2(); | |
806 | } | |
807 | void GetCorrectedSpectraLeonardo(TH1F* spectra,TH1F* correction, TH1F* hprimaryData,TH1F* hprimaryMC) | |
808 | { | |
809 | spectra->Sumw2(); | |
810 | spectra->Multiply(correction); | |
811 | spectra->Sumw2(); | |
812 | hprimaryData->Sumw2(); | |
813 | spectra->Multiply(hprimaryData); | |
814 | hprimaryMC->Sumw2(); | |
815 | spectra->Divide(hprimaryMC); | |
816 | } | |
817 | ||
b3ea73e1 | 818 | void GFCorrection(TH1F **Spectra,Float_t tofpt,UInt_t options) |
88d02ae4 | 819 | { |
b3ea73e1 | 820 | if (options&kgeantflukaKaon) |
821 | { | |
822 | //Geant/Fluka Correction | |
823 | Printf("\nGF correction for Kaons"); | |
824 | //Getting GF For Kaons in TPC | |
825 | TGraph *gGFCorrectionKaonPlus=new TGraph(); | |
826 | gGFCorrectionKaonPlus->SetName("gGFCorrectionKaonPlus"); | |
827 | gGFCorrectionKaonPlus->SetTitle("gGFCorrectionKaonPlus"); | |
828 | TGraph *gGFCorrectionKaonMinus=new TGraph(); | |
829 | gGFCorrectionKaonMinus->SetName("gGFCorrectionKaonMinus"); | |
830 | gGFCorrectionKaonMinus->SetTitle("gGFCorrectionKaonMinus"); | |
831 | TString fnameGeanFlukaK="GFCorrection/correctionForCrossSection.321.root"; | |
832 | TFile *fGeanFlukaK= new TFile(fnameGeanFlukaK.Data()); | |
833 | if (!fGeanFlukaK) | |
834 | return; | |
835 | TH1F *hGeantFlukaKPos=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles"); | |
836 | TH1F *hGeantFlukaKNeg=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles"); | |
837 | //getting GF func for Kaons with TOF | |
838 | TF1 *fGFKPosTracking; | |
839 | fGFKPosTracking = TrackingEff_geantflukaCorrection(3,kPositive); | |
840 | TF1 *fGFKNegTracking; | |
841 | fGFKNegTracking = TrackingEff_geantflukaCorrection(3,kNegative); | |
842 | TF1 *fGFKPosMatching; | |
843 | fGFKPosMatching = TOFmatchMC_geantflukaCorrection(3,kPositive); | |
844 | TF1 *fGFKNegMatching; | |
845 | fGFKNegMatching = TOFmatchMC_geantflukaCorrection(3,kNegative); | |
846 | for(Int_t binK=0;binK<=Spectra[1]->GetNbinsX();binK++) | |
847 | { | |
88d02ae4 | 848 | if(Spectra[1]->GetBinCenter(binK)<tofpt) |
849 | {//use TPC GeantFlukaCorrection | |
850 | Float_t FlukaCorrKPos=hGeantFlukaKPos->GetBinContent(hGeantFlukaKPos->FindBin(Spectra[1]->GetBinCenter(binK))); | |
851 | Float_t FlukaCorrKNeg=hGeantFlukaKNeg->GetBinContent(hGeantFlukaKNeg->FindBin(Spectra[4]->GetBinCenter(binK))); | |
852 | // Printf("TPC Geant/Fluka: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPos,FlukaCorrKNeg); | |
853 | Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPos); | |
854 | Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNeg); | |
855 | Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPos); | |
856 | Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNeg); | |
857 | gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),FlukaCorrKPos); | |
858 | gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),FlukaCorrKNeg); | |
859 | } | |
860 | else | |
861 | { | |
862 | gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),0); | |
863 | gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),0); | |
864 | Float_t FlukaCorrKPosTracking=fGFKPosTracking->Eval(Spectra[1]->GetBinCenter(binK)); | |
865 | Float_t FlukaCorrKNegTracking=fGFKNegTracking->Eval(Spectra[1]->GetBinCenter(binK)); | |
866 | // Printf("TPC/TOF Geant/Fluka Tracking: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosTracking,FlukaCorrKNegTracking); | |
867 | Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosTracking); | |
868 | Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegTracking); | |
869 | Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosTracking); | |
870 | Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegTracking); | |
871 | Float_t FlukaCorrKPosMatching=fGFKPosMatching->Eval(Spectra[1]->GetBinCenter(binK)); | |
872 | Float_t FlukaCorrKNegMatching=fGFKNegMatching->Eval(Spectra[1]->GetBinCenter(binK)); | |
873 | // Printf("TPC/TOF Geant/Fluka Matching: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosMatching,FlukaCorrKNegMatching); | |
874 | Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosMatching); | |
875 | Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegMatching); | |
876 | Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosMatching); | |
877 | Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegMatching); | |
878 | } | |
b3ea73e1 | 879 | } |
88d02ae4 | 880 | } |
b3ea73e1 | 881 | if(!(options&kgeantflukaProton)) |
882 | return; | |
88d02ae4 | 883 | //Geant Fluka for P in TPC |
884 | Printf("\nGF correction for Protons"); | |
885 | const Int_t kNCharge=2; | |
886 | Int_t kPos=0; | |
887 | Int_t kNeg=1; | |
888 | TFile* fGFProtons = new TFile ("GFCorrection/correctionForCrossSection.root"); | |
889 | TH2D * hCorrFluka[kNCharge]; | |
890 | TH2D * hCorrFluka[2]; | |
891 | hCorrFluka[kPos] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionProtons"); | |
892 | hCorrFluka[kNeg] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionAntiProtons"); | |
893 | //getting GF func for Kaons with TPCTOF | |
894 | TF1 *fGFpPosTracking; | |
895 | fGFpPosTracking = TrackingEff_geantflukaCorrection(4,kPositive); | |
896 | TF1 *fGFpNegTracking; | |
897 | fGFpNegTracking = TrackingEff_geantflukaCorrection(4,kNegative); | |
898 | TF1 *fGFpPosMatching; | |
899 | fGFpPosMatching = TOFmatchMC_geantflukaCorrection(4,kPositive); | |
900 | TF1 *fGFpNegMatching; | |
901 | fGFpNegMatching = TOFmatchMC_geantflukaCorrection(4,kNegative); | |
902 | ||
903 | ||
904 | Int_t nbins = Spectra[2]->GetNbinsX(); | |
905 | ||
906 | for(Int_t ibin = 0; ibin < nbins; ibin++) | |
907 | { | |
908 | if(Spectra[2]->GetBinCenter(ibin)<tofpt) | |
909 | {//use TPC GeantFlukaCorrection | |
910 | for(Int_t icharge = 0; icharge < kNCharge; icharge++) | |
911 | { | |
912 | Int_t nbinsy=hCorrFluka[icharge]->GetNbinsY(); | |
913 | Float_t pt = Spectra[2]->GetBinCenter(ibin); | |
914 | Float_t minPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(1); | |
915 | Float_t maxPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1); | |
916 | if (pt < minPtCorrection) | |
917 | pt = minPtCorrection+0.0001; | |
918 | if (pt > maxPtCorrection) | |
919 | pt = maxPtCorrection; | |
920 | Float_t correction = hCorrFluka[icharge]->GetBinContent(1,hCorrFluka[icharge]->GetYaxis()->FindBin(pt)); | |
921 | ||
922 | if (correction > 0.0) | |
923 | {// If the bin is empty this is a 0 | |
924 | // cout<<icharge<<" "<<ibin<<" "<<correction<<endl; | |
925 | Spectra[icharge*3+2]->SetBinContent(ibin,Spectra[icharge*3+2]->GetBinContent(ibin)*correction); | |
926 | Spectra[icharge*3+2]->SetBinError(ibin,Spectra[icharge*3+2]->GetBinError(ibin)*correction); | |
927 | } | |
928 | else if (Spectra[icharge*3+2]->GetBinContent(ibin) > 0.0) | |
929 | { | |
930 | // If we are skipping a non-empty bin, we notify the user | |
931 | cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for " <<"protons Pt:"<< pt<< endl; | |
932 | cout << " Bin content: " << Spectra[icharge*3+2]->GetBinContent(ibin) << endl; | |
933 | } | |
934 | } | |
935 | } | |
936 | else | |
937 | { | |
938 | Float_t FlukaCorrpPosTracking=fGFpPosTracking->Eval(Spectra[2]->GetBinCenter(ibin)); | |
939 | Float_t FlukaCorrpNegTracking=fGFpNegTracking->Eval(Spectra[2]->GetBinCenter(ibin)); | |
940 | // Printf("TPC/TOF Geant/Fluka Tracking: pt:%f Pos:%f Neg:%f",Spectra[2]->GetBinCenter(ibin),FlukaCorrpPosTracking,FlukaCorrpNegTracking); | |
941 | Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*FlukaCorrpPosTracking); | |
942 | Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*FlukaCorrpNegTracking); | |
943 | Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError(ibin)*FlukaCorrpPosTracking); | |
944 | Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError(ibin)*FlukaCorrpNegTracking); | |
945 | Float_t FlukaCorrpPosMatching=fGFpPosMatching->Eval(Spectra[2]->GetBinCenter(ibin)); | |
946 | Float_t FlukaCorrpNegMatching=fGFpNegMatching->Eval(Spectra[2]->GetBinCenter(ibin)); | |
947 | // Printf("TPC/TOF Geant/Fluka Matching: pt:%f Pos:%f Neg:%f",Spectra[2]->GetBinCenter(ibin),FlukaCorrpPosMatching,FlukaCorrpNegMatching); | |
948 | Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*FlukaCorrpPosMatching); | |
949 | Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*FlukaCorrpNegMatching); | |
950 | Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError(ibin)*FlukaCorrpPosMatching); | |
951 | Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError(ibin)*FlukaCorrpNegMatching); | |
952 | } | |
953 | } | |
954 | fGeanFlukaK->Close(); | |
955 | delete fGeanFlukaK; | |
956 | } | |
957 | ||
958 | ||
959 | /////////// | |
960 | TF1 * | |
961 | TrackingEff_geantflukaCorrection(Int_t ipart, Int_t icharge) | |
962 | { | |
963 | ||
964 | if (ipart == 3 && icharge == kNegative) { | |
965 | TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.); | |
966 | return f; | |
967 | } | |
968 | else if (ipart == 4 && icharge == kNegative) { | |
969 | TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.); | |
970 | } | |
971 | else | |
972 | TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionNull(x)", 0., 5.); | |
973 | ||
974 | return f; | |
975 | } | |
976 | ||
977 | Double_t | |
978 | TrackingPtGeantFlukaCorrectionNull(Double_t pTmc) | |
979 | { | |
980 | return 1.; | |
981 | } | |
982 | ||
983 | Double_t | |
984 | TrackingPtGeantFlukaCorrectionPrMinus(Double_t pTmc) | |
985 | { | |
986 | return (1 - 0.129758 *TMath::Exp(-pTmc*0.679612)); | |
987 | } | |
988 | ||
989 | Double_t | |
990 | TrackingPtGeantFlukaCorrectionKaMinus(Double_t pTmc) | |
991 | { | |
992 | return TMath::Min((0.972865 + 0.0117093*pTmc), 1.); | |
993 | } | |
994 | /////////////////////////////////////////// | |
995 | TF1 * | |
996 | TOFmatchMC_geantflukaCorrection(Int_t ipart, Int_t icharge) | |
997 | { | |
998 | ||
999 | if (ipart == 3 && icharge == kNegative) { | |
1000 | TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.); | |
1001 | return f; | |
1002 | } | |
1003 | else if (ipart == 4 && icharge == kNegative) { | |
1004 | TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.); | |
1005 | } | |
1006 | else | |
1007 | TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionNull(x)", 0., 5.); | |
1008 | ||
1009 | return f; | |
1010 | } | |
1011 | ||
1012 | ||
1013 | Double_t | |
1014 | MatchingPtGeantFlukaCorrectionNull(Double_t pTmc) | |
1015 | { | |
1016 | return 1.; | |
1017 | } | |
1018 | ||
1019 | Double_t | |
1020 | MatchingPtGeantFlukaCorrectionPrMinus(Double_t pTmc) | |
1021 | { | |
1022 | Float_t ptTPCoutP =pTmc*(1-6.81059e-01*TMath::Exp(-pTmc*4.20094)); | |
1023 | return (TMath::Power(1 - 0.129758*TMath::Exp(-ptTPCoutP*0.679612),0.07162/0.03471)); | |
1024 | } | |
1025 | ||
1026 | Double_t | |
1027 | MatchingPtGeantFlukaCorrectionKaMinus(Double_t pTmc) | |
1028 | { | |
1029 | Float_t ptTPCoutK=pTmc*(1- 3.37297e-03/pTmc/pTmc - 3.26544e-03/pTmc); | |
1030 | return TMath::Min((TMath::Power(0.972865 + 0.0117093*ptTPCoutK,0.07162/0.03471)), 1.); | |
1031 | } | |
1032 | void MatchingTOFEff(TH1F** Spectra, TList* list=0x0) | |
1033 | { | |
1034 | if(TOFMatchingScalling[0]<0.0&&TOFMatchingScalling[1]<0.0) | |
1035 | { | |
1036 | TH1F *hMatcEffPos_data=(TH1F*)tcutsdata->GetHistoNMatchedPos(); | |
1037 | hMatcEffPos_data->Divide((TH1F*)tcutsdata->GetHistoNSelectedPos()); | |
1038 | hMatcEffPos_data->SetTitle("Matching Eff Pos - data"); | |
1039 | TH1F *hMatcEffNeg_data=(TH1F*)tcutsdata->GetHistoNMatchedNeg(); | |
1040 | hMatcEffNeg_data->Divide((TH1F*)tcutsdata->GetHistoNSelectedNeg()); | |
1041 | hMatcEffNeg_data->SetTitle("Matching Eff Neg - data"); | |
1042 | TH1F *hMatcEffPos_mc=(TH1F*)tcutsmc->GetHistoNMatchedPos(); | |
1043 | hMatcEffPos_mc->Divide((TH1F*)tcutsmc->GetHistoNSelectedPos()); | |
1044 | hMatcEffPos_mc->SetTitle("Matching Eff Pos - mc"); | |
1045 | TH1F *hMatcEffNeg_mc=(TH1F*)tcutsmc->GetHistoNMatchedNeg(); | |
1046 | hMatcEffNeg_mc->Divide((TH1F*)tcutsmc->GetHistoNSelectedNeg()); | |
1047 | hMatcEffNeg_mc->SetTitle("Matching Eff Neg - mc"); | |
1048 | ||
1049 | ||
1050 | hMatcEffPos_data->Divide(hMatcEffPos_mc); | |
1051 | hMatcEffNeg_data->Divide(hMatcEffNeg_mc); | |
1052 | hMatcEffPos_data->SetName("MatchingTOFPos"); | |
1053 | hMatcEffNeg_data->SetName("MatchingTOFNeg"); | |
1054 | ||
1055 | ||
1056 | TF1 *pol0MatchPos_data=new TF1("pol0MatchPos_data","pol0",.6,5); | |
1057 | hMatcEffPos_data->Fit("pol0MatchPos_data","MNR"); | |
1058 | TF1 *pol0MatchNeg_data=new TF1("pol0MatchNeg_data","pol0",.6,5); | |
1059 | hMatcEffNeg_data->Fit("pol0MatchNeg_data","MNR"); | |
1060 | ||
1061 | list->Add(hMatcEffPos_data); | |
1062 | list->Add(hMatcEffNeg_data); | |
1063 | ||
1064 | ||
1065 | TOFMatchingScalling[0]=pol0MatchPos_data->GetParameter(0); | |
1066 | TOFMatchingScalling[1]=pol0MatchNeg_data->GetParameter(0); | |
1067 | } | |
1068 | //Correction spectra for matching efficiency | |
1069 | //For the moment I'm using the inclusive correction | |
1070 | for(Int_t ipart=0;ipart<3;ipart++) | |
1071 | { | |
1072 | ||
1073 | for(Int_t ibin=1;ibin<Spectra[ipart]->GetNbinsX();ibin++) | |
1074 | { | |
1075 | Float_t ptspectra=Spectra[ipart]->GetBinCenter(ibin); | |
1076 | if(ptspectra<tcutsdata->GetPtTOFMatching()) | |
1077 | continue; | |
1078 | //Spectra[ipart]->SetBinContent(ibin,( Spectra[ipart]->GetBinContent(ibin)/hMatcEffPos_data->GetBinContent(hMatcEffPos_data->FindBin(ptspectra)))); | |
1079 | //Spectra[ipart+3]->SetBinContent(ibin,( Spectra[ipart+3]->GetBinContent(ibin)/hMatcEffNeg_data->GetBinContent(hMatcEffNeg_data->FindBin(ptspectra)))); | |
1080 | Spectra[ipart]->SetBinContent(ibin,( Spectra[ipart]->GetBinContent(ibin)/TOFMatchingScalling[0])); | |
1081 | Spectra[ipart+3]->SetBinContent(ibin,( Spectra[ipart+3]->GetBinContent(ibin)/TOFMatchingScalling[1])); | |
1082 | } | |
1083 | } | |
1084 | } | |
1085 | Double_t eta2y(Double_t pt, Double_t mass, Double_t eta) | |
1086 | { | |
1087 | Double_t mt = TMath::Sqrt(mass * mass + pt * pt); | |
1088 | return TMath::ASinH(pt / mt * TMath::SinH(eta)); | |
1089 | } | |
1090 | ||
1091 | TH1* GetSumAllCh(TH1F** spectra, Double_t* mass ) | |
1092 | { | |
1093 | TH1F* allch=(((TH1F*))spectra[0]->Clone("allCh")); | |
1094 | allch->Reset(); | |
1095 | for (int i=0;i<6;i++) | |
1096 | { | |
1097 | Double_t masstmp=mass[i%3]; | |
1098 | for (int j=1;j<spectra[i]->GetXaxis()->GetNbins();j++) | |
1099 | { | |
1100 | Float_t value=spectra[i]->GetBinContent(j); | |
1101 | Float_t error=spectra[i]->GetBinError(j); | |
1102 | if(value>0.0) | |
1103 | { | |
1104 | Float_t pt=spectra[i]->GetXaxis()->GetBinCenter(j); | |
1105 | Float_t mt2=masstmp*masstmp+pt*pt; | |
1106 | Float_t mt=TMath::Sqrt(mt2); | |
1107 | Float_t maxy=eta2y(pt,masstmp,0.8); | |
1108 | 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)))); | |
1109 | conver=conver/1.6; | |
1110 | cout<<maxy<<" "<<conver<<" "<<masstmp<<""<<spectra[i]->GetName()<<endl; | |
1111 | Float_t bincontent=allch->GetBinContent(j); | |
1112 | Float_t binerror=allch->GetBinError(j); | |
1113 | bincontent+=conver*value; | |
1114 | binerror=TMath::Sqrt(binerror*binerror+conver*conver*error*error); | |
1115 | allch->SetBinContent(j,bincontent); | |
1116 | allch->SetBinError(j,binerror); | |
1117 | } | |
1118 | ||
1119 | } | |
1120 | } | |
1121 | Divideby2pipt(allch); | |
1122 | allch->SetTitle("N_{ch};p_{T} (GeV/c);1/(2 #pi p_{T})dN/p_{T} |#eta|<0.8"); | |
1123 | return allch; | |
1124 | } | |
1125 | ||
1126 | void Divideby2pipt(TH1* hist) | |
1127 | { | |
1128 | ||
1129 | for (int i=1;i<hist->GetXaxis()->GetNbins();i++) | |
1130 | { | |
1131 | Float_t value=hist->GetBinContent(i); | |
1132 | Float_t error=hist->GetBinError(i); | |
1133 | Float_t pt=hist->GetXaxis()->GetBinCenter(i); | |
1134 | hist->SetBinContent(i,value/(2.0*TMath::Pi()*pt)); | |
1135 | hist->SetBinError(i,error/(2.0*TMath::Pi()*pt)); | |
1136 | ||
1137 | } | |
1138 | } |