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