]>
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 | ||
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 |