Added getters and members for some event counters (A. Marin)
[u/mrichter/AliRoot.git] / ANALYSIS / AliCollisionNormalization.cxx
CommitLineData
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
22ClassImp(AliCollisionNormalization)
23
24const char * AliCollisionNormalization::fProcLabel[] = {"SD","DD","ND", "Unknown"};
25
26AliCollisionNormalization::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
66AliCollisionNormalization::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
104AliCollisionNormalization::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
169AliCollisionNormalization::~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
185void 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
211TH1 * 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
239TH2F * 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}
254TH2F * 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
272TH2F * 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
289Double_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
436Long64_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
494void 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
504void 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}
511void 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
519Int_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
575Double_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
617void 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