0d300b4f |
1 | //------------------------------------------------------------------------- |
2 | // Implementation of Class AliCollisionNormalization |
3 | // |
4 | // This class is used to store the vertex ditributions in the data |
5 | // and in Monte Carlo, needed to compute the real number of |
6 | // collisions a given sample is corresponding to. |
7 | // The strategy matches what described in CERN-THESIS-2009-033 p 119. |
8 | // |
9 | // Author: Michele Floris, CERN |
10 | //------------------------------------------------------------------------- |
11 | |
9b0cb3c3 |
12 | #include "AliCollisionNormalization.h" |
13 | #include "AliPhysicsSelection.h" |
14 | #include "AliLog.h" |
15 | #include "TFile.h" |
16 | #include "TCanvas.h" |
17 | #include "AliGenPythiaEventHeader.h" |
18 | #include "AliGenDPMjetEventHeader.h" |
19 | #include "AliGenEventHeader.h" |
20 | #include "AliMCEvent.h" |
21 | |
22 | ClassImp(AliCollisionNormalization) |
23 | |
24 | const char * AliCollisionNormalization::fProcLabel[] = {"SD","DD","ND", "Unknown"}; |
25 | |
26 | AliCollisionNormalization::AliCollisionNormalization() : |
27 | TObject(), |
28 | fNbinsVz(0), |
29 | fMinVz(0), |
30 | fMaxVz(0), |
31 | fZRange(9.99), |
32 | fIsMC(0), |
33 | fReferenceXS(0), |
34 | fVerbose(0), |
0d300b4f |
35 | fEnergy(900), |
9b0cb3c3 |
36 | fHistVzData (0), |
37 | fHistProcTypes (0), |
0d300b4f |
38 | fHistStatBin0 (0), |
7a0032a9 |
39 | fHistStat (0), |
40 | fInputEvents(0), |
41 | fPhysSelEvents(0), |
42 | fBgEvents(0), |
43 | fAllEvents(0), |
44 | fAllEventsZRange(0), |
45 | fAllEventsZRangeMult1(0), |
46 | fAllEventsInBin0ZRange(0), |
47 | fTrigEffBin0(0) |
9b0cb3c3 |
48 | { |
49 | |
50 | // ctor |
51 | |
52 | fNbinsVz = 30; |
53 | fMinVz = -15; |
54 | fMaxVz = +15; |
55 | |
56 | for(Int_t iproc = 0; iproc < kNProcs; iproc++){ |
57 | fHistVzMCGen[iproc] = 0; |
58 | fHistVzMCRec[iproc] = 0; |
59 | fHistVzMCTrg[iproc] = 0; |
60 | } |
61 | |
62 | |
63 | BookAllHistos(); |
64 | } |
65 | |
66 | AliCollisionNormalization::AliCollisionNormalization(Int_t nbinz, Float_t minz, Float_t maxz): |
67 | TObject(), |
68 | fNbinsVz(0), |
69 | fMinVz(0), |
70 | fMaxVz(0), |
71 | fZRange(9.99), |
72 | fIsMC(0), |
73 | fReferenceXS(0), |
74 | fVerbose(0), |
0d300b4f |
75 | fEnergy(900), |
9b0cb3c3 |
76 | fHistVzData (0), |
77 | fHistProcTypes (0), |
0d300b4f |
78 | fHistStatBin0 (0), |
7a0032a9 |
79 | fHistStat (0), |
80 | fInputEvents(0), |
81 | fPhysSelEvents(0), |
82 | fBgEvents(0), |
83 | fAllEvents(0), |
84 | fAllEventsZRange(0), |
85 | fAllEventsZRangeMult1(0), |
86 | fAllEventsInBin0ZRange(0), |
87 | fTrigEffBin0(0) |
9b0cb3c3 |
88 | { |
89 | |
90 | // ctor, allows setting binning |
91 | fNbinsVz = nbinz; |
92 | fMinVz = minz; |
93 | fMaxVz = maxz; |
94 | |
95 | for(Int_t iproc = 0; iproc < kNProcs; iproc++){ |
96 | fHistVzMCGen[iproc] = 0; |
97 | fHistVzMCRec[iproc] = 0; |
98 | fHistVzMCTrg[iproc] = 0; |
99 | } |
100 | |
101 | BookAllHistos(); |
102 | } |
103 | |
104 | AliCollisionNormalization::AliCollisionNormalization(const char * dataFile, const char * dataListName, |
105 | const char * mcFile, const char * mcListName, |
106 | const char * eventStatFile) : |
107 | TObject(), |
108 | fNbinsVz(0), |
109 | fMinVz(0), |
110 | fMaxVz(0), |
111 | fZRange(9.99), |
112 | fIsMC(0), |
113 | fReferenceXS(0), |
114 | fVerbose(0), |
0d300b4f |
115 | fEnergy(900), |
9b0cb3c3 |
116 | fHistVzData (0), |
117 | fHistProcTypes (0), |
0d300b4f |
118 | fHistStatBin0 (0), |
7a0032a9 |
119 | fHistStat (0), |
120 | fInputEvents(0), |
121 | fPhysSelEvents(0), |
122 | fBgEvents(0), |
123 | fAllEvents(0), |
124 | fAllEventsZRange(0), |
125 | fAllEventsZRangeMult1(0), |
126 | fAllEventsInBin0ZRange(0), |
127 | fTrigEffBin0(0) |
0d300b4f |
128 | |
9b0cb3c3 |
129 | { |
130 | |
131 | // ctor, loads histograms from file |
132 | for(Int_t iproc = 0; iproc < kNProcs; iproc++){ |
133 | fHistVzMCGen[iproc] = 0; |
134 | fHistVzMCRec[iproc] = 0; |
135 | fHistVzMCTrg[iproc] = 0; |
136 | } |
137 | |
138 | |
139 | TFile * fdata = new TFile (dataFile); |
140 | TFile * fmc = new TFile (mcFile ); |
0d300b4f |
141 | TFile * fstat = new TFile(eventStatFile); |
142 | |
143 | if (!fdata->IsOpen() || !fmc->IsOpen() || !fstat->IsOpen()) { |
144 | AliFatal("Cannot open input file(s)"); |
145 | } |
9b0cb3c3 |
146 | |
147 | TList * ldata = (TList*) fdata->Get(dataListName); |
148 | TList * lmc = (TList*) fmc ->Get(mcListName ); |
149 | |
150 | AliCollisionNormalization * cndata = (AliCollisionNormalization*) ldata->FindObject("AliCollisionNormalization"); |
151 | AliCollisionNormalization * cnmc = (AliCollisionNormalization*) lmc ->FindObject("AliCollisionNormalization"); |
152 | |
153 | |
154 | // Assign or book all histos |
155 | for(Int_t iproc = 0; iproc < kNProcs; iproc++){ |
156 | fHistVzMCGen[iproc]= cnmc->GetVzMCGen(iproc); |
157 | fHistVzMCRec[iproc]= cnmc->GetVzMCRec(iproc); |
158 | fHistVzMCTrg[iproc]= cnmc->GetVzMCTrg(iproc); |
159 | } |
160 | fHistVzData = cndata->GetVzData(); |
161 | fHistProcTypes = cnmc->GetHistProcTypes(); |
162 | |
9b0cb3c3 |
163 | fHistStatBin0 = (TH1F*) fstat->Get("fHistStatistics_Bin0"); |
0d300b4f |
164 | fHistStat = (TH1F*) fstat->Get("fHistStatistics"); |
9b0cb3c3 |
165 | |
166 | } |
167 | |
168 | |
169 | AliCollisionNormalization::~AliCollisionNormalization(){ |
170 | |
171 | // dtor |
172 | for(Int_t iproc = 0; iproc < kNProcs; iproc++){ |
173 | if(fHistVzMCGen[iproc]) { delete fHistVzMCGen[iproc] ; fHistVzMCGen[iproc] =0;} |
174 | if(fHistVzMCRec[iproc]) { delete fHistVzMCRec[iproc] ; fHistVzMCRec[iproc] =0;} |
175 | if(fHistVzMCTrg[iproc]) { delete fHistVzMCTrg[iproc] ; fHistVzMCTrg[iproc] =0;} |
176 | } |
177 | |
178 | if(fHistVzData ) { delete fHistVzData ; fHistVzData =0;} |
179 | if(fHistStatBin0 ) { delete fHistStatBin0 ; fHistStatBin0 =0;} |
0d300b4f |
180 | if(fHistStat ) { delete fHistStat ; fHistStat =0;} |
9b0cb3c3 |
181 | if(fHistProcTypes ) { delete fHistProcTypes ; fHistProcTypes =0;} |
182 | |
183 | } |
184 | |
185 | void AliCollisionNormalization::BookAllHistos(){ |
186 | // books all histos |
187 | // Book histos of vz distributions vs multiplicity |
188 | // if vzOnly == kTRUE, it returns a 1d histo with vz dist only |
189 | |
190 | // Do not attach histos to the current directory |
191 | Bool_t oldStatus = TH1::AddDirectoryStatus(); |
192 | TH1::AddDirectory(kFALSE); |
193 | |
194 | for(Int_t iproc = 0; iproc < kNProcs; iproc++){ |
195 | fHistVzMCGen [iproc] = (TH2F*) BookVzHisto(TString("fHistVzMCGen")+ fProcLabel[iproc] ,"Vz distribution of generated events vs rec multiplicity "); |
196 | fHistVzMCRec [iproc] = (TH2F*) BookVzHisto(TString("fHistVzMCRec")+ fProcLabel[iproc] ,"Vz distribution of reconstructed events vs rec multiplicity"); |
197 | fHistVzMCTrg [iproc] = (TH2F*) BookVzHisto(TString("fHistVzMCTrg")+ fProcLabel[iproc] ,"Vz distribution of triggered events vs rec multiplicity "); |
198 | } |
199 | fHistVzData = (TH2F*) BookVzHisto("fHistVzData" ,"Vz distribution of triggered events vs rec multiplicity "); |
200 | fHistProcTypes = new TH1F ("fHistProcTypes", "Number of events in the different process classes", kNProcs, -0.5 , kNProcs-0.5); |
201 | |
202 | fHistProcTypes->GetXaxis()->SetBinLabel(kProcSD+1,"SD"); |
203 | fHistProcTypes->GetXaxis()->SetBinLabel(kProcND+1,"ND"); |
204 | fHistProcTypes->GetXaxis()->SetBinLabel(kProcDD+1,"DD"); |
205 | fHistProcTypes->GetXaxis()->SetBinLabel(kProcUnknown+1,"Unknown"); |
206 | |
207 | TH1::AddDirectory(oldStatus); |
208 | |
209 | } |
210 | |
211 | TH1 * AliCollisionNormalization::BookVzHisto(const char * name , const char * title, Bool_t vzOnly) { |
212 | |
213 | // Book histos of vz distributions vs multiplicity |
214 | // if vzOnly == kTRUE, it returns a 1d histo with vz dist only |
215 | |
216 | |
217 | // Do not attach histos to the current directory |
218 | Bool_t oldStatus = TH1::AddDirectoryStatus(); |
219 | TH1::AddDirectory(kFALSE); |
220 | |
221 | TH1 * h; |
222 | Double_t binLimitsVtx[] = {-30,-25,-20,-15,-10,-7,-5.5,-4,-3,-2,-1,0,1,2,3,4,5.5,7,10,15,20,25,30}; |
223 | if (vzOnly) { |
224 | h = new TH1F(name,title,22,binLimitsVtx); |
225 | } else { |
226 | h = new TH2F(name,title,22,binLimitsVtx,100,Double_t(-0.5),Double_t(99.5)); |
227 | } |
228 | |
229 | h->SetXTitle("V_{z} (cm)"); |
230 | h->SetYTitle("n_{trk}"); |
231 | h->Sumw2(); |
232 | |
233 | TH1::AddDirectory(oldStatus); |
234 | |
235 | return h; |
236 | |
237 | } |
238 | |
239 | TH2F * AliCollisionNormalization::GetVzMCGen (Int_t procType) { |
240 | |
241 | // returns MC gen histo. If proc type < 0 sums over all proc types, reweighting the XS |
242 | |
243 | if(procType>=0) return fHistVzMCGen[procType] ; |
244 | |
245 | TH2F * sum = (TH2F*) fHistVzMCGen[kProcSD]->Clone(); |
246 | sum->Reset(); |
247 | |
248 | for(Int_t iproc = 0; iproc < kProcUnknown; iproc++){ |
249 | sum->Add(fHistVzMCGen[iproc],GetProcessWeight(iproc)); |
250 | } |
251 | |
252 | return sum; |
253 | } |
254 | TH2F * AliCollisionNormalization::GetVzMCRec (Int_t procType) { |
255 | |
256 | // returns MC rec histo. If proc type < 0 sums over all proc types, reweighting the XS |
257 | |
258 | if(procType>=0) return fHistVzMCRec[procType] ; |
259 | |
260 | TH2F * sum = (TH2F*) fHistVzMCRec[kProcSD]->Clone(); |
261 | sum->Reset(); |
262 | |
263 | for(Int_t iproc = 0; iproc < kProcUnknown; iproc++){ |
264 | sum->Add(fHistVzMCRec[iproc],GetProcessWeight(iproc)); |
265 | } |
266 | |
267 | return sum; |
268 | |
269 | } |
270 | |
271 | |
272 | TH2F * AliCollisionNormalization::GetVzMCTrg (Int_t procType) { |
273 | |
274 | // returns MC trg histo. If proc type < 0 sums over all proc types, reweighting the XS |
275 | |
276 | if(procType>=0) return fHistVzMCTrg[procType] ; |
277 | |
278 | TH2F * sum = (TH2F*) fHistVzMCTrg[kProcSD]->Clone(); |
279 | sum->Reset(); |
280 | |
281 | for(Int_t iproc = 0; iproc < kProcUnknown; iproc++){ |
282 | sum->Add(fHistVzMCTrg[iproc],GetProcessWeight(iproc)); |
283 | } |
284 | |
285 | return sum; |
286 | |
287 | } |
288 | |
289 | Double_t AliCollisionNormalization::ComputeNint() { |
290 | |
291 | // Compute the number of collisions applying all corrections |
0d300b4f |
292 | // TODO: check error propagation |
9b0cb3c3 |
293 | |
294 | TH1 * hVzData = fHistVzData->ProjectionX("_px",2,-1,"E"); // Skip zero bin |
295 | Int_t allEventsWithVertex = (Int_t) fHistVzData->Integral(0, fHistVzData->GetNbinsX()+1, |
296 | 2, fHistVzData->GetNbinsY()+1); // include under/overflow!, skip 0 bin |
297 | |
298 | // Assign histos reweighted with measured XS |
299 | TH2F * histVzMCGen = GetVzMCGen(-1); |
300 | TH2F * histVzMCTrg = GetVzMCTrg(-1); |
301 | TH2F * histVzMCRec = GetVzMCRec(-1); |
302 | |
0d300b4f |
303 | // Before we start: print number of input events to the physics |
304 | // selection: this allows to crosscheck that all runs were |
305 | // successfully processed (useful if running on grid: you may have a |
306 | // crash without noticing it). |
7a0032a9 |
307 | |
308 | fInputEvents = fHistStat->GetBinContent(1,1); |
309 | fPhysSelEvents = fHistStat->GetBinContent( fHistStat->GetNbinsX(),1); |
310 | |
0d300b4f |
311 | AliInfo(Form("Input Events (No cuts: %d, After Phys. Sel.:%d)", |
7a0032a9 |
312 | Int_t(fInputEvents), |
313 | Int_t(fPhysSelEvents))); // Fixme: will this change with a different trigger scheme? |
0d300b4f |
314 | |
9b0cb3c3 |
315 | // Get or compute BG. This assumes the CINT1B suite |
316 | Double_t triggeredEventsWith0MultWithBG = fHistVzData->Integral(0, fHistVzData->GetNbinsX()+1,1, 1); |
7a0032a9 |
317 | // Double_t bg = 0; // This will include beam gas + accidentals |
9b0cb3c3 |
318 | if (fHistStatBin0->GetNbinsY() > 4) { // FIXME: we need a better criterion to decide... |
319 | AliInfo("Using BG computed by Physics Selection"); |
46918b3b |
320 | // WARNING |
321 | // CHECK IF THE COMPUTATION OF BG OFFSET IS STILL OK, IN CASE OF CHANGES TO THE PHYSICS SELECTION |
322 | Int_t bgOffset = 0; |
323 | Int_t nbiny = fHistStatBin0->GetNbinsY(); |
324 | for(Int_t ibiny =1; ibiny <= nbiny; ibiny++){ |
325 | bgOffset++; |
326 | printf("%d, %s\n", bgOffset, fHistStatBin0->GetYaxis()->GetBinLabel(ibiny) ); |
327 | |
328 | if(!strncmp("All B", fHistStatBin0->GetYaxis()->GetBinLabel(ibiny),5)) break; |
329 | if((ibiny+1)==nbiny) AliFatal("Cannot compute bg offset"); |
330 | } |
331 | |
332 | if(fVerbose > 2) AliInfo(Form("BG Offset: %d",bgOffset)); |
333 | |
334 | |
7a0032a9 |
335 | fBgEvents = fHistStatBin0->GetBinContent(fHistStatBin0->GetNbinsX(), bgOffset+AliPhysicsSelection::kStatRowBG); |
336 | fBgEvents += fHistStatBin0->GetBinContent(fHistStatBin0->GetNbinsX(),bgOffset+AliPhysicsSelection::kStatRowAcc); |
0d300b4f |
337 | Int_t cint1B = (Int_t) fHistStatBin0->GetBinContent(fHistStatBin0->GetNbinsX(),1); |
338 | if (cint1B != Int_t(triggeredEventsWith0MultWithBG)) { |
339 | AliWarning(Form("Events in bin0 from physics selection and local counter not consistent: %d - %d", cint1B, Int_t(triggeredEventsWith0MultWithBG))); |
340 | } |
9b0cb3c3 |
341 | } else { |
342 | AliInfo("Computing BG using CINT1A/B/C/E, ignoring intensities"); |
46918b3b |
343 | AliError("This will only work for early runs!!! USE BG FROM PHYSICS SELECTION INSTEAD"); |
9b0cb3c3 |
344 | Int_t icol = fHistStatBin0->GetNbinsX(); |
345 | Int_t cint1B = (Int_t) fHistStatBin0->GetBinContent(icol,1); |
346 | Int_t cint1A = (Int_t) fHistStatBin0->GetBinContent(icol,2); |
347 | Int_t cint1C = (Int_t) fHistStatBin0->GetBinContent(icol,3); |
348 | Int_t cint1E = (Int_t) fHistStatBin0->GetBinContent(icol,4); |
7a0032a9 |
349 | fBgEvents = cint1A + cint1C-2*cint1E ; // FIXME: to be changed to take into account ratios of events |
9b0cb3c3 |
350 | if (cint1B != triggeredEventsWith0MultWithBG) { |
e8b839ab |
351 | AliWarning(Form("Events in bin0 from physics selection and local counter not consistent: %d - %d", cint1B, (Int_t)triggeredEventsWith0MultWithBG)); |
9b0cb3c3 |
352 | } |
353 | } |
354 | |
7a0032a9 |
355 | Double_t triggeredEventsWith0Mult = triggeredEventsWith0MultWithBG - fBgEvents; |
9b0cb3c3 |
356 | if(fVerbose > 0) AliInfo(Form("Measured events with vertex: %d",allEventsWithVertex)); |
357 | Double_t bin0 = fHistVzData->Integral(0, fHistVzData->GetNbinsX()+1, |
358 | 1, 1); |
7a0032a9 |
359 | if(fVerbose > 0) AliInfo(Form("Zero Bin, Meas: %2.2f, BG: %2.2f, Meas - BG: %2.2f", bin0, fBgEvents, triggeredEventsWith0Mult)); |
9b0cb3c3 |
360 | |
361 | // This pointers are here so that I use the same names used in Jan Fiete's code (allows for faster/easier comparison) |
362 | |
363 | TH2* eTrig = histVzMCTrg; |
364 | TH2* eTrigVtx = histVzMCRec; |
365 | TH1* eTrigVtx_projx = eTrigVtx->ProjectionX("eTrigVtx_projx", 2, eTrigVtx->GetNbinsX()+1); |
366 | |
367 | // compute trigger and vertex efficiency |
368 | TH2 * effVtxTrig = (TH2*) histVzMCRec->Clone("effVtxTrig"); |
369 | effVtxTrig->Reset(); |
370 | effVtxTrig->Divide(histVzMCRec,histVzMCGen,1,1,"B"); |
371 | // Apply correction to data to get corrected events |
372 | TH2 * correctedEvents = (TH2*) fHistVzData->Clone("correctedEvents"); |
373 | correctedEvents->Divide(effVtxTrig); |
374 | |
375 | // TH1 * correctedEvents = fHistVzData->ProjectionX("eTrigVtx_projx", 2, eTrigVtx->GetNbinsX()+1); |
376 | |
377 | // loop over vertex bins |
378 | for (Int_t i = 1; i <= fHistVzData->GetNbinsX(); i++) |
379 | { |
380 | Double_t alpha = (Double_t) hVzData->GetBinContent(i) / allEventsWithVertex; |
381 | Double_t events = alpha * triggeredEventsWith0Mult; |
382 | |
b2f8342e |
383 | // if (histVzMCRec->GetBinContent(i,1) == 0) |
384 | // continue; |
385 | if (!eTrigVtx_projx->GetBinContent(i) || !eTrig->Integral(0, eTrig->GetNbinsX()+1, 1, 1)) |
9b0cb3c3 |
386 | continue; |
387 | |
388 | Double_t fZ = eTrigVtx_projx->Integral(0, eTrigVtx_projx->GetNbinsX()+1) / eTrigVtx_projx->GetBinContent(i) * |
389 | eTrig->GetBinContent(i, 1) / eTrig->Integral(0, eTrig->GetNbinsX()+1, 1, 1); |
390 | |
391 | events *= fZ; |
392 | |
393 | // multiply with trigger correction |
394 | Double_t correction = histVzMCGen->GetBinContent(i,1)/histVzMCTrg->GetBinContent(i,1); |
395 | events *= correction; |
396 | |
397 | if (fVerbose > 1) Printf(" Bin %d, alpha is %.2f%%, fZ is %.3f, number of events with 0 mult.: %.2f (MC comparison: %.2f)", i, alpha * 100., fZ, events, |
b2f8342e |
398 | histVzMCTrg->GetBinContent(i,1)); |
9b0cb3c3 |
399 | correctedEvents->SetBinContent(i, 1, events); |
400 | } |
401 | |
402 | |
403 | |
404 | // Integrate correctedEvents over full z range |
7a0032a9 |
405 | fAllEvents = correctedEvents->Integral(0, correctedEvents->GetNbinsX()+1,0, correctedEvents->GetNbinsY()+1); |
9b0cb3c3 |
406 | // Integrate correctedEvents over needed z range |
7a0032a9 |
407 | fAllEventsZRange = correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange), |
9b0cb3c3 |
408 | 0, correctedEvents->GetNbinsY()+1); |
7a0032a9 |
409 | |
410 | fAllEventsZRangeMult1 = correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange), |
411 | 2,correctedEvents->GetNbinsX()+1); |
412 | |
413 | fAllEventsInBin0ZRange = correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange),1,1); |
414 | |
9b0cb3c3 |
415 | if(fVerbose > 1) AliInfo(Form("Results in |Vz| < %3.3f",fZRange)); |
416 | if(fVerbose > 1) AliInfo(Form(" Events in Bin0: %2.2f, With > 1 track: %2.2f, All corrected: %2.2f", |
7a0032a9 |
417 | fAllEventsInBin0ZRange, |
418 | fAllEventsZRangeMult1, |
419 | fAllEventsZRange)); |
420 | |
9b0cb3c3 |
421 | if(fVerbose > 1) { |
b2f8342e |
422 | Int_t nbin = histVzMCTrg->GetNbinsX(); |
7a0032a9 |
423 | fTrigEffBin0 = histVzMCTrg->Integral(0,nbin+1,1,1)/histVzMCGen->Integral(0,nbin+1,1,1); |
424 | AliInfo(Form("Trigger Efficiency in the zero bin: %3.3f", fTrigEffBin0 )); |
9b0cb3c3 |
425 | } |
426 | |
0d300b4f |
427 | |
7a0032a9 |
428 | AliInfo(Form("Number of collisions in full phase space: %2.2f", fAllEvents)); |
0d300b4f |
429 | // effVtxTrig->Draw(); |
430 | // new TCanvas(); |
431 | // correctedEvents->Draw(); // FIXME: debug |
9b0cb3c3 |
432 | |
7a0032a9 |
433 | return fAllEvents; |
9b0cb3c3 |
434 | } |
435 | |
436 | Long64_t AliCollisionNormalization::Merge(TCollection* list) |
437 | { |
438 | // Merge a list of AliCollisionNormalization objects with this |
439 | // (needed for PROOF). |
440 | // Returns the number of merged objects (including this). |
441 | |
442 | if (!list) |
443 | return 0; |
444 | |
445 | if (list->IsEmpty()) |
446 | return 1; |
447 | |
448 | TIterator* iter = list->MakeIterator(); |
449 | TObject* obj; |
450 | |
451 | // collections of all histograms |
0d300b4f |
452 | const Int_t nHists = kNProcs*3+5; |
9b0cb3c3 |
453 | TList collections[nHists]; |
454 | |
455 | Int_t count = 0; |
456 | while ((obj = iter->Next())) { |
457 | |
458 | AliCollisionNormalization* entry = dynamic_cast<AliCollisionNormalization*> (obj); |
459 | if (entry == 0) |
460 | continue; |
461 | Int_t ihist = -1; |
462 | |
463 | for(Int_t iproc = 0; iproc < kNProcs; iproc++){ |
464 | if (entry->fHistVzMCGen[iproc] ) collections[++ihist].Add(entry->fHistVzMCGen[iproc] ); |
465 | if (entry->fHistVzMCRec[iproc] ) collections[++ihist].Add(entry->fHistVzMCRec[iproc] ); |
466 | if (entry->fHistVzMCTrg[iproc] ) collections[++ihist].Add(entry->fHistVzMCTrg[iproc] ); |
467 | } |
468 | if (entry->fHistVzData ) collections[++ihist].Add(entry->fHistVzData ); |
469 | if (entry->fHistProcTypes ) collections[++ihist].Add(entry->fHistProcTypes ); |
470 | if (entry->fHistStatBin0 ) collections[++ihist].Add(entry->fHistStatBin0 ); |
0d300b4f |
471 | if (entry->fHistStat ) collections[++ihist].Add(entry->fHistStat ); |
9b0cb3c3 |
472 | |
473 | count++; |
474 | } |
475 | |
476 | Int_t ihist = -1; |
477 | for(Int_t iproc = 0; iproc < kNProcs; iproc++){ |
478 | if (fHistVzMCGen[iproc] ) fHistVzMCGen[iproc] ->Merge(&collections[++ihist]); |
479 | if (fHistVzMCRec[iproc] ) fHistVzMCRec[iproc] ->Merge(&collections[++ihist]); |
480 | if (fHistVzMCTrg[iproc] ) fHistVzMCTrg[iproc] ->Merge(&collections[++ihist]); |
481 | |
482 | } |
483 | if (fHistVzData ) fHistVzData ->Merge(&collections[++ihist]); |
484 | if (fHistProcTypes ) fHistProcTypes ->Merge(&collections[++ihist]); |
485 | if (fHistStatBin0 ) fHistStatBin0 ->Merge(&collections[++ihist]); |
0d300b4f |
486 | if (fHistStat ) fHistStat ->Merge(&collections[++ihist]); |
9b0cb3c3 |
487 | |
488 | |
489 | delete iter; |
490 | |
491 | return count+1; |
492 | } |
493 | |
494 | void AliCollisionNormalization::FillVzMCGen(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt) { |
495 | |
496 | // Fill MC gen histo and the process types statistics |
497 | |
498 | Int_t evtype = GetProcessType(mcEvt); |
499 | // When I fill gen histos, I also fill statistics of process types (used to detemine ratios afterwards). |
500 | fHistProcTypes->Fill(evtype); |
501 | fHistVzMCGen[evtype]->Fill(vz,ntrk); |
502 | } |
503 | |
504 | void AliCollisionNormalization::FillVzMCRec(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt){ |
505 | |
506 | // Fill MC rec histo |
507 | Int_t evtype = GetProcessType(mcEvt); |
508 | fHistVzMCRec[evtype]->Fill(vz,ntrk); |
509 | |
510 | } |
511 | void AliCollisionNormalization::FillVzMCTrg(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt) { |
512 | |
513 | // Fill MC trg histo |
514 | Int_t evtype = GetProcessType(mcEvt); |
515 | fHistVzMCTrg[evtype]->Fill(vz,ntrk); |
516 | } |
517 | |
518 | |
519 | Int_t AliCollisionNormalization::GetProcessType(AliMCEvent * mcEvt) { |
520 | |
521 | // Determine if the event was generated with pythia or phojet and return the process type |
522 | |
523 | // Check if mcEvt is fine |
524 | if (!mcEvt) { |
525 | AliFatal("NULL mc event"); |
526 | } |
527 | |
528 | // Determine if it was a pythia or phojet header, and return the correct process type |
529 | AliGenPythiaEventHeader * headPy = 0; |
530 | AliGenDPMjetEventHeader * headPho = 0; |
531 | AliGenEventHeader * htmp = mcEvt->GenEventHeader(); |
532 | if(!htmp) { |
533 | AliFatal("Cannot Get MC Header!!"); |
534 | } |
535 | if( TString(htmp->IsA()->GetName()) == "AliGenPythiaEventHeader") { |
536 | headPy = (AliGenPythiaEventHeader*) htmp; |
537 | } else if (TString(htmp->IsA()->GetName()) == "AliGenDPMjetEventHeader") { |
538 | headPho = (AliGenDPMjetEventHeader*) htmp; |
539 | } else { |
540 | AliError("Unknown header"); |
541 | } |
542 | |
543 | // Determine process type |
544 | if(headPy) { |
545 | if(headPy->ProcessType() == 92 || headPy->ProcessType() == 93) { |
546 | // single difractive |
547 | return kProcSD; |
548 | } else if (headPy->ProcessType() == 94) { |
549 | // double diffractive |
550 | return kProcDD; |
551 | } |
552 | else if(headPy->ProcessType() != 92 && headPy->ProcessType() != 93 && headPy->ProcessType() != 94) { |
553 | // non difractive |
554 | return kProcND; |
555 | } |
556 | } else if (headPho) { |
557 | if(headPho->ProcessType() == 5 || headPho->ProcessType() == 6 ) { |
558 | // single difractive |
559 | return kProcSD; |
560 | } else if (headPho->ProcessType() == 7) { |
561 | // double diffractive |
562 | return kProcDD; |
563 | } else if(headPho->ProcessType() != 5 && headPho->ProcessType() != 6 && headPho->ProcessType() != 7 ) { |
564 | // non difractive |
565 | return kProcND; |
566 | } |
567 | } |
568 | |
569 | |
570 | // no process type found? |
571 | AliError(Form("Unknown header: %s", htmp->IsA()->GetName())); |
572 | return kProcUnknown; |
573 | } |
574 | |
575 | Double_t AliCollisionNormalization::GetProcessWeight(Int_t proc) { |
576 | |
577 | // Return a weight to be used when combining the MC histos to |
0d300b4f |
578 | // compute efficiency, defined as the ratio XS in generator / XS |
579 | // measured |
9b0cb3c3 |
580 | |
581 | Float_t ref_SD, ref_DD, ref_ND, error_SD, error_DD, error_ND; |
582 | GetRelativeFractions(fReferenceXS,ref_SD, ref_DD, ref_ND, error_SD, error_DD, error_ND); |
583 | |
584 | static Double_t total = fHistProcTypes->Integral(); |
585 | if (fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcUnknown)) > 0) { |
586 | AliError("There were unknown processes!!!"); |
587 | } |
588 | static Double_t SD = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcSD))/total; |
589 | static Double_t DD = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcDD))/total; |
590 | static Double_t ND = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcND))/total; |
591 | |
592 | if (fVerbose > 2) { |
593 | AliInfo(Form("Total MC evts: %f",total)); |
594 | AliInfo(Form(" Frac SD %4.4f", SD)); |
595 | AliInfo(Form(" Frac DD %4.4f", DD)); |
596 | AliInfo(Form(" Frac ND %4.4f", ND)); |
597 | } |
598 | |
599 | switch(proc) { |
600 | case kProcSD: |
601 | return ref_SD/SD; |
602 | break; |
603 | case kProcDD: |
604 | return ref_DD/DD; |
605 | break; |
606 | case kProcND: |
607 | return ref_ND/ND; |
608 | break; |
609 | default: |
610 | AliError("Unknown process"); |
611 | } |
612 | |
613 | return 0; |
614 | |
615 | } |
616 | |
617 | void AliCollisionNormalization::GetRelativeFractions(Int_t origin, Float_t& ref_SD, Float_t& ref_DD, Float_t& ref_ND, Float_t& error_SD, Float_t& error_DD, Float_t& error_ND) |
618 | { |
619 | // Returns fraction of XS (SD, ND, DD) and corresponding error |
620 | // Stolen from Jan Fiete's drawSystematics macro |
621 | |
622 | // origin: |
623 | // -1 = Pythia (test) |
624 | // 0 = UA5 |
625 | // 1 = Data 1.8 TeV |
626 | // 2 = Tel-Aviv |
627 | // 3 = Durham |
628 | // |
629 | |
0d300b4f |
630 | Double_t epsilon = 0.0001; |
631 | if(TMath::Abs(fEnergy-900)<epsilon) { |
632 | |
633 | switch (origin) |
634 | { |
635 | case -10: // Pythia default at 7 GeV, 50% error |
636 | AliInfo("PYTHIA x-sections"); |
637 | ref_SD = 0.192637; error_SD = ref_SD * 0.5; |
638 | ref_DD = 0.129877; error_DD = ref_DD * 0.5; |
639 | ref_ND = 0.677486; error_ND = 0; |
640 | break; |
641 | |
642 | case -1: // Pythia default at 900 GeV, as test |
643 | AliInfo("PYTHIA x-sections"); |
9b0cb3c3 |
644 | ref_SD = 0.223788; |
645 | ref_DD = 0.123315; |
646 | ref_ND = 0.652897; |
647 | break; |
648 | |
0d300b4f |
649 | case 0: // UA5 |
650 | AliInfo("UA5 x-sections a la first paper"); |
651 | ref_SD = 0.153; error_SD = 0.05; |
652 | ref_DD = 0.080; error_DD = 0.05; |
653 | ref_ND = 0.767; error_ND = 0; |
654 | break; |
655 | |
656 | case 10: // UA5 |
657 | AliInfo("UA5 x-sections hadron level definition for Pythia"); |
658 | // Fractions in Pythia with UA5 cuts selection for SD |
659 | // ND: 0.688662 |
660 | // SD: 0.188588 --> this should be 15.3 |
661 | // DD: 0.122750 |
662 | ref_SD = 0.224 * 0.153 / 0.189; error_SD = 0.023 * 0.224 / 0.189; |
663 | ref_DD = 0.095; error_DD = 0.06; |
664 | ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0; |
665 | break; |
666 | |
667 | case 11: // UA5 |
668 | AliInfo("UA5 x-sections hadron level definition for Phojet"); |
669 | // Fractions in Phojet with UA5 cuts selection for SD |
670 | // ND: 0.783573 |
671 | // SD: 0.151601 --> this should be 15.3 |
672 | // DD: 0.064827 |
673 | ref_SD = 0.191 * 0.153 / 0.152; error_SD = 0.023 * 0.191 / 0.152; |
674 | ref_DD = 0.095; error_DD = 0.06; |
675 | ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0; |
676 | break; |
677 | case 2: // tel-aviv model |
678 | AliInfo("Tel-aviv model x-sections"); |
679 | ref_SD = 0.171; |
680 | ref_DD = 0.094; |
681 | ref_ND = 1 - ref_SD - ref_DD; |
682 | break; |
683 | |
684 | case 3: // durham model |
685 | AliInfo("Durham model x-sections"); |
686 | ref_SD = 0.190; |
687 | ref_DD = 0.125; |
688 | ref_ND = 1 - ref_SD - ref_DD; |
9b0cb3c3 |
689 | break; |
0d300b4f |
690 | default: |
691 | AliFatal(Form("Unknown origin %d, Energy %f", origin, fEnergy)); |
692 | } |
693 | } |
694 | else if(TMath::Abs(fEnergy-1800)<epsilon) { |
695 | switch (origin) |
696 | { |
697 | case 20: // E710, 1.8 TeV |
9b0cb3c3 |
698 | AliInfo("E710 x-sections hadron level definition for Pythia"); |
699 | // ND: 0.705709 |
700 | // SD: 0.166590 --> this should be 15.9 |
701 | // DD: 0.127701 |
702 | ref_SD = 0.217 * 0.159 / 0.167; error_SD = 0.024 * 0.217 / 0.167; |
703 | ref_DD = 0.075 * 1.43; error_DD = 0.02 * 1.43; |
704 | ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0; |
705 | break; |
706 | |
0d300b4f |
707 | case 21: // E710, 1.8 TeV |
708 | AliInfo("E710 x-sections hadron level definition for Phojet"); |
709 | // ND: 0.817462 |
710 | // SD: 0.125506 --> this should be 15.9 |
711 | // DD: 0.057032 |
712 | ref_SD = 0.161 * 0.159 / 0.126; error_SD = 0.024 * 0.161 / 0.126; |
713 | ref_DD = 0.075 * 1.43; error_DD = 0.02 * 1.43; |
714 | ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0; |
715 | break; |
716 | |
717 | case 1: // data 1.8 TeV |
718 | AliInfo("??? x-sections"); |
719 | ref_SD = 0.152; |
720 | ref_DD = 0.092; |
721 | ref_ND = 1 - ref_SD - ref_DD; |
722 | break; |
723 | default: |
724 | AliFatal(Form("Unknown origin %d, Energy %f", origin, fEnergy)); |
725 | } |
726 | } |
727 | else { |
e8b839ab |
728 | AliFatal(Form("Unknown energy %f", fEnergy)); |
9b0cb3c3 |
729 | } |
0d300b4f |
730 | |
9b0cb3c3 |
731 | } |
732 | |
0d300b4f |
733 | |