]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - ANALYSIS/AliCollisionNormalization.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Fix for productions in which the ProductionInfo exists, but does not contain the...
[u/mrichter/AliRoot.git] / ANALYSIS / AliCollisionNormalization.cxx
... / ...
CommitLineData
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
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::fgkProcLabel[] = {"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),
35 fEnergy(900),
36 fHistVzData (0),
37 fHistProcTypes (0),
38 fHistStatBin0 (0),
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)
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),
75 fEnergy(900),
76 fHistVzData (0),
77 fHistProcTypes (0),
78 fHistStatBin0 (0),
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)
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),
115 fEnergy(900),
116 fHistVzData (0),
117 fHistProcTypes (0),
118 fHistStatBin0 (0),
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)
128
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 );
141 TFile * fstat = new TFile(eventStatFile);
142
143 if (!fdata->IsOpen() || !fmc->IsOpen() || !fstat->IsOpen()) {
144 AliFatal("Cannot open input file(s)");
145 }
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
163 fHistStatBin0 = (TH1F*) fstat->Get("fHistStatistics_Bin0");
164 fHistStat = (TH1F*) fstat->Get("fHistStatistics");
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;}
180 if(fHistStat ) { delete fHistStat ; fHistStat =0;}
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")+ 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 ");
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
292 // TODO: check error propagation
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
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).
307
308 fInputEvents = fHistStat->GetBinContent(1,1);
309 fPhysSelEvents = fHistStat->GetBinContent( fHistStat->GetNbinsX(),1);
310
311 AliInfo(Form("Input Events (No cuts: %d, After Phys. Sel.:%d)",
312 Int_t(fInputEvents),
313 Int_t(fPhysSelEvents))); // Fixme: will this change with a different trigger scheme?
314
315 // Get or compute BG. This assumes the CINT1B suite
316 Double_t triggeredEventsWith0MultWithBG = fHistVzData->Integral(0, fHistVzData->GetNbinsX()+1,1, 1);
317 // Double_t bg = 0; // This will include beam gas + accidentals
318 if (fHistStatBin0->GetNbinsY() > 4) { // FIXME: we need a better criterion to decide...
319 AliInfo("Using BG computed by Physics Selection");
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
335 fBgEvents = fHistStatBin0->GetBinContent(fHistStatBin0->GetNbinsX(), bgOffset+AliPhysicsSelection::kStatRowBG);
336 fBgEvents += fHistStatBin0->GetBinContent(fHistStatBin0->GetNbinsX(),bgOffset+AliPhysicsSelection::kStatRowAcc);
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 }
341 } else {
342 AliInfo("Computing BG using CINT1A/B/C/E, ignoring intensities");
343 AliError("This will only work for early runs!!! USE BG FROM PHYSICS SELECTION INSTEAD");
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);
349 fBgEvents = cint1A + cint1C-2*cint1E ; // FIXME: to be changed to take into account ratios of events
350 if (cint1B != triggeredEventsWith0MultWithBG) {
351 AliWarning(Form("Events in bin0 from physics selection and local counter not consistent: %d - %d", cint1B, (Int_t)triggeredEventsWith0MultWithBG));
352 }
353 }
354
355 Double_t triggeredEventsWith0Mult = triggeredEventsWith0MultWithBG - fBgEvents;
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);
359 if(fVerbose > 0) AliInfo(Form("Zero Bin, Meas: %2.2f, BG: %2.2f, Meas - BG: %2.2f", bin0, fBgEvents, triggeredEventsWith0Mult));
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* eTrigVtxProjx = eTrigVtx->ProjectionX("eTrigVtxProjx", 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("eTrigVtxProjx", 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
383 // if (histVzMCRec->GetBinContent(i,1) == 0)
384 // continue;
385 if (!eTrigVtxProjx->GetBinContent(i) || !eTrig->Integral(0, eTrig->GetNbinsX()+1, 1, 1))
386 continue;
387
388 Double_t fZ = eTrigVtxProjx->Integral(0, eTrigVtxProjx->GetNbinsX()+1) / eTrigVtxProjx->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,
398 histVzMCTrg->GetBinContent(i,1));
399 correctedEvents->SetBinContent(i, 1, events);
400 }
401
402
403
404 // Integrate correctedEvents over full z range
405 fAllEvents = correctedEvents->Integral(0, correctedEvents->GetNbinsX()+1,0, correctedEvents->GetNbinsY()+1);
406 // Integrate correctedEvents over needed z range
407 fAllEventsZRange = correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange-0.0001),
408 0, correctedEvents->GetNbinsY()+1);
409
410 fAllEventsZRangeMult1 = correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange-0.0001),
411 2,correctedEvents->GetNbinsY()+1);
412
413 fAllEventsInBin0ZRange = correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange-0.0001),1,1);
414
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",
417 fAllEventsInBin0ZRange,
418 fAllEventsZRangeMult1,
419 fAllEventsZRange));
420 Int_t nbin = histVzMCTrg->GetNbinsX();
421 fTrigEffBin0 = histVzMCTrg->Integral(0,nbin+1,1,1)/histVzMCGen->Integral(0,nbin+1,1,1);
422 if(fVerbose > 1) {
423 AliInfo(Form("Trigger Efficiency in the zero bin: %3.3f", fTrigEffBin0 ));
424 }
425
426
427 AliInfo(Form("Number of collisions in full phase space: %2.2f", fAllEvents));
428// effVtxTrig->Draw();
429// new TCanvas();
430// correctedEvents->Draw(); // FIXME: debug
431
432 return fAllEvents;
433}
434
435Long64_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
451 const Int_t nHists = kNProcs*3+5;
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 );
470 if (entry->fHistStat ) collections[++ihist].Add(entry->fHistStat );
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]);
485 if (fHistStat ) fHistStat ->Merge(&collections[++ihist]);
486
487
488 delete iter;
489
490 return count+1;
491}
492
493void 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
503void 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}
510void 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
518Int_t AliCollisionNormalization::GetProcessType(const AliMCEvent * mcEvt) {
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
574Double_t AliCollisionNormalization::GetProcessWeight(Int_t proc) {
575
576 // Return a weight to be used when combining the MC histos to
577 // compute efficiency, defined as the ratio XS in generator / XS
578 // measured
579
580 Float_t refSD, refDD, refND, errorSD, errorDD, errorND;
581 GetRelativeFractions(fReferenceXS,refSD, refDD, refND, errorSD, errorDD, errorND);
582
583 static Double_t total = fHistProcTypes->Integral();
584 if (fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcUnknown)) > 0) {
585 AliError("There were unknown processes!!!");
586 }
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;
590
591 if (fVerbose > 2) {
592 AliInfo(Form("Total MC evts: %f",total));
593 AliInfo(Form(" Frac SD %4.4f", lSD));
594 AliInfo(Form(" Frac DD %4.4f", lDD));
595 AliInfo(Form(" Frac ND %4.4f", lND));
596 }
597
598 switch(proc) {
599 case kProcSD:
600 return refSD/lSD;
601 break;
602 case kProcDD:
603 return refDD/lDD;
604 break;
605 case kProcND:
606 return refND/lND;
607 break;
608 default:
609 AliError("Unknown process");
610 }
611
612 return 0;
613
614}
615
616void AliCollisionNormalization::GetRelativeFractions(Int_t origin, Float_t& refSD, Float_t& refDD, Float_t& refND, Float_t& errorSD, Float_t& errorDD, Float_t& errorND)
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
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");
636 refSD = 0.192637; errorSD = refSD * 0.5;
637 refDD = 0.129877; errorDD = refDD * 0.5;
638 refND = 0.677486; errorND = 0;
639 break;
640
641 case -1: // Pythia default at 900 GeV, as test
642 AliInfo("PYTHIA x-sections");
643 refSD = 0.223788;
644 refDD = 0.123315;
645 refND = 0.652897;
646 break;
647
648 case 0: // UA5
649 AliInfo("UA5 x-sections a la first paper");
650 refSD = 0.153; errorSD = 0.05;
651 refDD = 0.080; errorDD = 0.05;
652 refND = 0.767; errorND = 0;
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
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;
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
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;
675 break;
676 case 2: // tel-aviv model
677 AliInfo("Tel-aviv model x-sections");
678 refSD = 0.171;
679 refDD = 0.094;
680 refND = 1 - refSD - refDD;
681 break;
682
683 case 3: // durham model
684 AliInfo("Durham model x-sections");
685 refSD = 0.190;
686 refDD = 0.125;
687 refND = 1 - refSD - refDD;
688 break;
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
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
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;
704 break;
705
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
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;
714 break;
715
716 case 1: // data 1.8 TeV
717 AliInfo("??? x-sections");
718 refSD = 0.152;
719 refDD = 0.092;
720 refND = 1 - refSD - refDD;
721 break;
722 default:
723 AliFatal(Form("Unknown origin %d, Energy %f", origin, fEnergy));
724 }
725 }
726 else {
727 AliFatal(Form("Unknown energy %f", fEnergy));
728 }
729
730}
731
732