]>
Commit | Line | Data |
---|---|---|
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 | ||
3aa0e5b3 | 24 | const char * AliCollisionNormalization::fgkProcLabel[] = {"SD","DD","ND", "Unknown"}; |
9b0cb3c3 | 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++){ | |
3aa0e5b3 | 195 | fHistVzMCGen [iproc] = (TH2F*) BookVzHisto(TString("fHistVzMCGen")+ fgkProcLabel[iproc] ,"Vz distribution of generated events vs rec multiplicity "); |
196 | fHistVzMCRec [iproc] = (TH2F*) BookVzHisto(TString("fHistVzMCRec")+ fgkProcLabel[iproc] ,"Vz distribution of reconstructed events vs rec multiplicity"); | |
197 | fHistVzMCTrg [iproc] = (TH2F*) BookVzHisto(TString("fHistVzMCTrg")+ fgkProcLabel[iproc] ,"Vz distribution of triggered events vs rec multiplicity "); | |
9b0cb3c3 | 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; | |
3aa0e5b3 | 365 | TH1* eTrigVtxProjx = eTrigVtx->ProjectionX("eTrigVtxProjx", 2, eTrigVtx->GetNbinsX()+1); |
9b0cb3c3 | 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 | ||
3aa0e5b3 | 375 | // TH1 * correctedEvents = fHistVzData->ProjectionX("eTrigVtxProjx", 2, eTrigVtx->GetNbinsX()+1); |
9b0cb3c3 | 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; | |
3aa0e5b3 | 385 | if (!eTrigVtxProjx->GetBinContent(i) || !eTrig->Integral(0, eTrig->GetNbinsX()+1, 1, 1)) |
9b0cb3c3 | 386 | continue; |
387 | ||
3aa0e5b3 | 388 | Double_t fZ = eTrigVtxProjx->Integral(0, eTrigVtxProjx->GetNbinsX()+1) / eTrigVtxProjx->GetBinContent(i) * |
9b0cb3c3 | 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 |
59a02c2e | 407 | fAllEventsZRange = correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange-0.0001), |
9b0cb3c3 | 408 | 0, correctedEvents->GetNbinsY()+1); |
7a0032a9 | 409 | |
59a02c2e | 410 | fAllEventsZRangeMult1 = correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange-0.0001), |
411 | 2,correctedEvents->GetNbinsY()+1); | |
7a0032a9 | 412 | |
59a02c2e | 413 | fAllEventsInBin0ZRange = correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange-0.0001),1,1); |
7a0032a9 | 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)); | |
9c4a15a7 | 420 | Int_t nbin = histVzMCTrg->GetNbinsX(); |
421 | fTrigEffBin0 = histVzMCTrg->Integral(0,nbin+1,1,1)/histVzMCGen->Integral(0,nbin+1,1,1); | |
9b0cb3c3 | 422 | if(fVerbose > 1) { |
7a0032a9 | 423 | AliInfo(Form("Trigger Efficiency in the zero bin: %3.3f", fTrigEffBin0 )); |
9b0cb3c3 | 424 | } |
425 | ||
0d300b4f | 426 | |
7a0032a9 | 427 | AliInfo(Form("Number of collisions in full phase space: %2.2f", fAllEvents)); |
0d300b4f | 428 | // effVtxTrig->Draw(); |
429 | // new TCanvas(); | |
430 | // correctedEvents->Draw(); // FIXME: debug | |
9b0cb3c3 | 431 | |
7a0032a9 | 432 | return fAllEvents; |
9b0cb3c3 | 433 | } |
434 | ||
435 | Long64_t AliCollisionNormalization::Merge(TCollection* list) | |
436 | { | |
437 | // Merge a list of AliCollisionNormalization objects with this | |
438 | // (needed for PROOF). | |
439 | // Returns the number of merged objects (including this). | |
440 | ||
441 | if (!list) | |
442 | return 0; | |
443 | ||
444 | if (list->IsEmpty()) | |
445 | return 1; | |
446 | ||
447 | TIterator* iter = list->MakeIterator(); | |
448 | TObject* obj; | |
449 | ||
450 | // collections of all histograms | |
0d300b4f | 451 | const Int_t nHists = kNProcs*3+5; |
9b0cb3c3 | 452 | TList collections[nHists]; |
453 | ||
454 | Int_t count = 0; | |
455 | while ((obj = iter->Next())) { | |
456 | ||
457 | AliCollisionNormalization* entry = dynamic_cast<AliCollisionNormalization*> (obj); | |
458 | if (entry == 0) | |
459 | continue; | |
460 | Int_t ihist = -1; | |
461 | ||
462 | for(Int_t iproc = 0; iproc < kNProcs; iproc++){ | |
463 | if (entry->fHistVzMCGen[iproc] ) collections[++ihist].Add(entry->fHistVzMCGen[iproc] ); | |
464 | if (entry->fHistVzMCRec[iproc] ) collections[++ihist].Add(entry->fHistVzMCRec[iproc] ); | |
465 | if (entry->fHistVzMCTrg[iproc] ) collections[++ihist].Add(entry->fHistVzMCTrg[iproc] ); | |
466 | } | |
467 | if (entry->fHistVzData ) collections[++ihist].Add(entry->fHistVzData ); | |
468 | if (entry->fHistProcTypes ) collections[++ihist].Add(entry->fHistProcTypes ); | |
469 | if (entry->fHistStatBin0 ) collections[++ihist].Add(entry->fHistStatBin0 ); | |
0d300b4f | 470 | if (entry->fHistStat ) collections[++ihist].Add(entry->fHistStat ); |
9b0cb3c3 | 471 | |
472 | count++; | |
473 | } | |
474 | ||
475 | Int_t ihist = -1; | |
476 | for(Int_t iproc = 0; iproc < kNProcs; iproc++){ | |
477 | if (fHistVzMCGen[iproc] ) fHistVzMCGen[iproc] ->Merge(&collections[++ihist]); | |
478 | if (fHistVzMCRec[iproc] ) fHistVzMCRec[iproc] ->Merge(&collections[++ihist]); | |
479 | if (fHistVzMCTrg[iproc] ) fHistVzMCTrg[iproc] ->Merge(&collections[++ihist]); | |
480 | ||
481 | } | |
482 | if (fHistVzData ) fHistVzData ->Merge(&collections[++ihist]); | |
483 | if (fHistProcTypes ) fHistProcTypes ->Merge(&collections[++ihist]); | |
484 | if (fHistStatBin0 ) fHistStatBin0 ->Merge(&collections[++ihist]); | |
0d300b4f | 485 | if (fHistStat ) fHistStat ->Merge(&collections[++ihist]); |
9b0cb3c3 | 486 | |
487 | ||
488 | delete iter; | |
489 | ||
490 | return count+1; | |
491 | } | |
492 | ||
493 | void AliCollisionNormalization::FillVzMCGen(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt) { | |
494 | ||
495 | // Fill MC gen histo and the process types statistics | |
496 | ||
497 | Int_t evtype = GetProcessType(mcEvt); | |
498 | // When I fill gen histos, I also fill statistics of process types (used to detemine ratios afterwards). | |
499 | fHistProcTypes->Fill(evtype); | |
500 | fHistVzMCGen[evtype]->Fill(vz,ntrk); | |
501 | } | |
502 | ||
503 | void AliCollisionNormalization::FillVzMCRec(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt){ | |
504 | ||
505 | // Fill MC rec histo | |
506 | Int_t evtype = GetProcessType(mcEvt); | |
507 | fHistVzMCRec[evtype]->Fill(vz,ntrk); | |
508 | ||
509 | } | |
510 | void AliCollisionNormalization::FillVzMCTrg(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt) { | |
511 | ||
512 | // Fill MC trg histo | |
513 | Int_t evtype = GetProcessType(mcEvt); | |
514 | fHistVzMCTrg[evtype]->Fill(vz,ntrk); | |
515 | } | |
516 | ||
517 | ||
3aa0e5b3 | 518 | Int_t AliCollisionNormalization::GetProcessType(const AliMCEvent * mcEvt) { |
9b0cb3c3 | 519 | |
520 | // Determine if the event was generated with pythia or phojet and return the process type | |
521 | ||
522 | // Check if mcEvt is fine | |
523 | if (!mcEvt) { | |
524 | AliFatal("NULL mc event"); | |
525 | } | |
526 | ||
527 | // Determine if it was a pythia or phojet header, and return the correct process type | |
528 | AliGenPythiaEventHeader * headPy = 0; | |
529 | AliGenDPMjetEventHeader * headPho = 0; | |
530 | AliGenEventHeader * htmp = mcEvt->GenEventHeader(); | |
531 | if(!htmp) { | |
532 | AliFatal("Cannot Get MC Header!!"); | |
533 | } | |
534 | if( TString(htmp->IsA()->GetName()) == "AliGenPythiaEventHeader") { | |
535 | headPy = (AliGenPythiaEventHeader*) htmp; | |
536 | } else if (TString(htmp->IsA()->GetName()) == "AliGenDPMjetEventHeader") { | |
537 | headPho = (AliGenDPMjetEventHeader*) htmp; | |
538 | } else { | |
539 | AliError("Unknown header"); | |
540 | } | |
541 | ||
542 | // Determine process type | |
543 | if(headPy) { | |
544 | if(headPy->ProcessType() == 92 || headPy->ProcessType() == 93) { | |
545 | // single difractive | |
546 | return kProcSD; | |
547 | } else if (headPy->ProcessType() == 94) { | |
548 | // double diffractive | |
549 | return kProcDD; | |
550 | } | |
551 | else if(headPy->ProcessType() != 92 && headPy->ProcessType() != 93 && headPy->ProcessType() != 94) { | |
552 | // non difractive | |
553 | return kProcND; | |
554 | } | |
555 | } else if (headPho) { | |
556 | if(headPho->ProcessType() == 5 || headPho->ProcessType() == 6 ) { | |
557 | // single difractive | |
558 | return kProcSD; | |
559 | } else if (headPho->ProcessType() == 7) { | |
560 | // double diffractive | |
561 | return kProcDD; | |
562 | } else if(headPho->ProcessType() != 5 && headPho->ProcessType() != 6 && headPho->ProcessType() != 7 ) { | |
563 | // non difractive | |
564 | return kProcND; | |
565 | } | |
566 | } | |
567 | ||
568 | ||
569 | // no process type found? | |
570 | AliError(Form("Unknown header: %s", htmp->IsA()->GetName())); | |
571 | return kProcUnknown; | |
572 | } | |
573 | ||
574 | Double_t AliCollisionNormalization::GetProcessWeight(Int_t proc) { | |
575 | ||
576 | // Return a weight to be used when combining the MC histos to | |
0d300b4f | 577 | // compute efficiency, defined as the ratio XS in generator / XS |
578 | // measured | |
9b0cb3c3 | 579 | |
3aa0e5b3 | 580 | Float_t refSD, refDD, refND, errorSD, errorDD, errorND; |
581 | GetRelativeFractions(fReferenceXS,refSD, refDD, refND, errorSD, errorDD, errorND); | |
9b0cb3c3 | 582 | |
583 | static Double_t total = fHistProcTypes->Integral(); | |
584 | if (fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcUnknown)) > 0) { | |
585 | AliError("There were unknown processes!!!"); | |
586 | } | |
3aa0e5b3 | 587 | static Double_t lSD = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcSD))/total; |
588 | static Double_t lDD = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcDD))/total; | |
589 | static Double_t lND = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcND))/total; | |
9b0cb3c3 | 590 | |
591 | if (fVerbose > 2) { | |
592 | AliInfo(Form("Total MC evts: %f",total)); | |
3aa0e5b3 | 593 | AliInfo(Form(" Frac SD %4.4f", lSD)); |
594 | AliInfo(Form(" Frac DD %4.4f", lDD)); | |
595 | AliInfo(Form(" Frac ND %4.4f", lND)); | |
9b0cb3c3 | 596 | } |
597 | ||
598 | switch(proc) { | |
599 | case kProcSD: | |
3aa0e5b3 | 600 | return refSD/lSD; |
9b0cb3c3 | 601 | break; |
602 | case kProcDD: | |
3aa0e5b3 | 603 | return refDD/lDD; |
9b0cb3c3 | 604 | break; |
605 | case kProcND: | |
3aa0e5b3 | 606 | return refND/lND; |
9b0cb3c3 | 607 | break; |
608 | default: | |
609 | AliError("Unknown process"); | |
610 | } | |
611 | ||
612 | return 0; | |
613 | ||
614 | } | |
615 | ||
3aa0e5b3 | 616 | void AliCollisionNormalization::GetRelativeFractions(Int_t origin, Float_t& refSD, Float_t& refDD, Float_t& refND, Float_t& errorSD, Float_t& errorDD, Float_t& errorND) |
9b0cb3c3 | 617 | { |
618 | // Returns fraction of XS (SD, ND, DD) and corresponding error | |
619 | // Stolen from Jan Fiete's drawSystematics macro | |
620 | ||
621 | // origin: | |
622 | // -1 = Pythia (test) | |
623 | // 0 = UA5 | |
624 | // 1 = Data 1.8 TeV | |
625 | // 2 = Tel-Aviv | |
626 | // 3 = Durham | |
627 | // | |
628 | ||
0d300b4f | 629 | Double_t epsilon = 0.0001; |
630 | if(TMath::Abs(fEnergy-900)<epsilon) { | |
631 | ||
632 | switch (origin) | |
633 | { | |
634 | case -10: // Pythia default at 7 GeV, 50% error | |
635 | AliInfo("PYTHIA x-sections"); | |
3aa0e5b3 | 636 | refSD = 0.192637; errorSD = refSD * 0.5; |
637 | refDD = 0.129877; errorDD = refDD * 0.5; | |
638 | refND = 0.677486; errorND = 0; | |
0d300b4f | 639 | break; |
640 | ||
641 | case -1: // Pythia default at 900 GeV, as test | |
642 | AliInfo("PYTHIA x-sections"); | |
3aa0e5b3 | 643 | refSD = 0.223788; |
644 | refDD = 0.123315; | |
645 | refND = 0.652897; | |
9b0cb3c3 | 646 | break; |
647 | ||
0d300b4f | 648 | case 0: // UA5 |
649 | AliInfo("UA5 x-sections a la first paper"); | |
3aa0e5b3 | 650 | refSD = 0.153; errorSD = 0.05; |
651 | refDD = 0.080; errorDD = 0.05; | |
652 | refND = 0.767; errorND = 0; | |
0d300b4f | 653 | break; |
654 | ||
655 | case 10: // UA5 | |
656 | AliInfo("UA5 x-sections hadron level definition for Pythia"); | |
657 | // Fractions in Pythia with UA5 cuts selection for SD | |
658 | // ND: 0.688662 | |
659 | // SD: 0.188588 --> this should be 15.3 | |
660 | // DD: 0.122750 | |
3aa0e5b3 | 661 | refSD = 0.224 * 0.153 / 0.189; errorSD = 0.023 * 0.224 / 0.189; |
662 | refDD = 0.095; errorDD = 0.06; | |
663 | refND = 1.0 - refSD - refDD; errorND = 0; | |
0d300b4f | 664 | break; |
665 | ||
666 | case 11: // UA5 | |
667 | AliInfo("UA5 x-sections hadron level definition for Phojet"); | |
668 | // Fractions in Phojet with UA5 cuts selection for SD | |
669 | // ND: 0.783573 | |
670 | // SD: 0.151601 --> this should be 15.3 | |
671 | // DD: 0.064827 | |
3aa0e5b3 | 672 | refSD = 0.191 * 0.153 / 0.152; errorSD = 0.023 * 0.191 / 0.152; |
673 | refDD = 0.095; errorDD = 0.06; | |
674 | refND = 1.0 - refSD - refDD; errorND = 0; | |
0d300b4f | 675 | break; |
676 | case 2: // tel-aviv model | |
677 | AliInfo("Tel-aviv model x-sections"); | |
3aa0e5b3 | 678 | refSD = 0.171; |
679 | refDD = 0.094; | |
680 | refND = 1 - refSD - refDD; | |
0d300b4f | 681 | break; |
682 | ||
683 | case 3: // durham model | |
684 | AliInfo("Durham model x-sections"); | |
3aa0e5b3 | 685 | refSD = 0.190; |
686 | refDD = 0.125; | |
687 | refND = 1 - refSD - refDD; | |
9b0cb3c3 | 688 | break; |
0d300b4f | 689 | default: |
690 | AliFatal(Form("Unknown origin %d, Energy %f", origin, fEnergy)); | |
691 | } | |
692 | } | |
693 | else if(TMath::Abs(fEnergy-1800)<epsilon) { | |
694 | switch (origin) | |
695 | { | |
696 | case 20: // E710, 1.8 TeV | |
9b0cb3c3 | 697 | AliInfo("E710 x-sections hadron level definition for Pythia"); |
698 | // ND: 0.705709 | |
699 | // SD: 0.166590 --> this should be 15.9 | |
700 | // DD: 0.127701 | |
3aa0e5b3 | 701 | refSD = 0.217 * 0.159 / 0.167; errorSD = 0.024 * 0.217 / 0.167; |
702 | refDD = 0.075 * 1.43; errorDD = 0.02 * 1.43; | |
703 | refND = 1.0 - refSD - refDD; errorND = 0; | |
9b0cb3c3 | 704 | break; |
705 | ||
0d300b4f | 706 | case 21: // E710, 1.8 TeV |
707 | AliInfo("E710 x-sections hadron level definition for Phojet"); | |
708 | // ND: 0.817462 | |
709 | // SD: 0.125506 --> this should be 15.9 | |
710 | // DD: 0.057032 | |
3aa0e5b3 | 711 | refSD = 0.161 * 0.159 / 0.126; errorSD = 0.024 * 0.161 / 0.126; |
712 | refDD = 0.075 * 1.43; errorDD = 0.02 * 1.43; | |
713 | refND = 1.0 - refSD - refDD; errorND = 0; | |
0d300b4f | 714 | break; |
715 | ||
716 | case 1: // data 1.8 TeV | |
717 | AliInfo("??? x-sections"); | |
3aa0e5b3 | 718 | refSD = 0.152; |
719 | refDD = 0.092; | |
720 | refND = 1 - refSD - refDD; | |
0d300b4f | 721 | break; |
722 | default: | |
723 | AliFatal(Form("Unknown origin %d, Energy %f", origin, fEnergy)); | |
724 | } | |
725 | } | |
726 | else { | |
e8b839ab | 727 | AliFatal(Form("Unknown energy %f", fEnergy)); |
9b0cb3c3 | 728 | } |
0d300b4f | 729 | |
9b0cb3c3 | 730 | } |
731 | ||
0d300b4f | 732 |