Added loop for extraction of clusters really attached to its track.
[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
3aa0e5b3 24const char * AliCollisionNormalization::fgkProcLabel[] = {"SD","DD","ND", "Unknown"};
9b0cb3c3 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++){
3aa0e5b3 195 fHistVzMCGen [iproc] = (TH2F*) BookVzHisto(TString("fHistVzMCGen")+ fgkProcLabel[iproc] ,"Vz distribution of generated events vs rec multiplicity ");
196 fHistVzMCRec [iproc] = (TH2F*) BookVzHisto(TString("fHistVzMCRec")+ fgkProcLabel[iproc] ,"Vz distribution of reconstructed events vs rec multiplicity");
197 fHistVzMCTrg [iproc] = (TH2F*) BookVzHisto(TString("fHistVzMCTrg")+ fgkProcLabel[iproc] ,"Vz distribution of triggered events vs rec multiplicity ");
9b0cb3c3 198 }
199 fHistVzData = (TH2F*) BookVzHisto("fHistVzData" ,"Vz distribution of triggered events vs rec multiplicity ");
200 fHistProcTypes = new TH1F ("fHistProcTypes", "Number of events in the different process classes", kNProcs, -0.5 , kNProcs-0.5);
201
202 fHistProcTypes->GetXaxis()->SetBinLabel(kProcSD+1,"SD");
203 fHistProcTypes->GetXaxis()->SetBinLabel(kProcND+1,"ND");
204 fHistProcTypes->GetXaxis()->SetBinLabel(kProcDD+1,"DD");
205 fHistProcTypes->GetXaxis()->SetBinLabel(kProcUnknown+1,"Unknown");
206
207 TH1::AddDirectory(oldStatus);
208
209}
210
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;
3aa0e5b3 365 TH1* eTrigVtxProjx = eTrigVtx->ProjectionX("eTrigVtxProjx", 2, eTrigVtx->GetNbinsX()+1);
9b0cb3c3 366
367 // compute trigger and vertex efficiency
368 TH2 * effVtxTrig = (TH2*) histVzMCRec->Clone("effVtxTrig");
369 effVtxTrig->Reset();
370 effVtxTrig->Divide(histVzMCRec,histVzMCGen,1,1,"B");
371 // Apply correction to data to get corrected events
372 TH2 * correctedEvents = (TH2*) fHistVzData->Clone("correctedEvents");
373 correctedEvents->Divide(effVtxTrig);
374
3aa0e5b3 375 // TH1 * correctedEvents = fHistVzData->ProjectionX("eTrigVtxProjx", 2, eTrigVtx->GetNbinsX()+1);
9b0cb3c3 376
377 // loop over vertex bins
378 for (Int_t i = 1; i <= fHistVzData->GetNbinsX(); i++)
379 {
380 Double_t alpha = (Double_t) hVzData->GetBinContent(i) / allEventsWithVertex;
381 Double_t events = alpha * triggeredEventsWith0Mult;
382
b2f8342e 383 // if (histVzMCRec->GetBinContent(i,1) == 0)
384 // continue;
3aa0e5b3 385 if (!eTrigVtxProjx->GetBinContent(i) || !eTrig->Integral(0, eTrig->GetNbinsX()+1, 1, 1))
9b0cb3c3 386 continue;
387
3aa0e5b3 388 Double_t fZ = eTrigVtxProjx->Integral(0, eTrigVtxProjx->GetNbinsX()+1) / eTrigVtxProjx->GetBinContent(i) *
9b0cb3c3 389 eTrig->GetBinContent(i, 1) / eTrig->Integral(0, eTrig->GetNbinsX()+1, 1, 1);
390
391 events *= fZ;
392
393 // multiply with trigger correction
394 Double_t correction = histVzMCGen->GetBinContent(i,1)/histVzMCTrg->GetBinContent(i,1);
395 events *= correction;
396
397 if (fVerbose > 1) Printf(" Bin %d, alpha is %.2f%%, fZ is %.3f, number of events with 0 mult.: %.2f (MC comparison: %.2f)", i, alpha * 100., fZ, events,
b2f8342e 398 histVzMCTrg->GetBinContent(i,1));
9b0cb3c3 399 correctedEvents->SetBinContent(i, 1, events);
400 }
401
402
403
404 // Integrate correctedEvents over full z range
7a0032a9 405 fAllEvents = correctedEvents->Integral(0, correctedEvents->GetNbinsX()+1,0, correctedEvents->GetNbinsY()+1);
9b0cb3c3 406 // Integrate correctedEvents over needed z range
59a02c2e 407 fAllEventsZRange = correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange-0.0001),
9b0cb3c3 408 0, correctedEvents->GetNbinsY()+1);
7a0032a9 409
59a02c2e 410 fAllEventsZRangeMult1 = correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange-0.0001),
411 2,correctedEvents->GetNbinsY()+1);
7a0032a9 412
59a02c2e 413 fAllEventsInBin0ZRange = correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange-0.0001),1,1);
7a0032a9 414
9b0cb3c3 415 if(fVerbose > 1) AliInfo(Form("Results in |Vz| < %3.3f",fZRange));
416 if(fVerbose > 1) AliInfo(Form(" Events in Bin0: %2.2f, With > 1 track: %2.2f, All corrected: %2.2f",
7a0032a9 417 fAllEventsInBin0ZRange,
418 fAllEventsZRangeMult1,
419 fAllEventsZRange));
9c4a15a7 420 Int_t nbin = histVzMCTrg->GetNbinsX();
421 fTrigEffBin0 = histVzMCTrg->Integral(0,nbin+1,1,1)/histVzMCGen->Integral(0,nbin+1,1,1);
9b0cb3c3 422 if(fVerbose > 1) {
7a0032a9 423 AliInfo(Form("Trigger Efficiency in the zero bin: %3.3f", fTrigEffBin0 ));
9b0cb3c3 424 }
425
0d300b4f 426
7a0032a9 427 AliInfo(Form("Number of collisions in full phase space: %2.2f", fAllEvents));
0d300b4f 428// effVtxTrig->Draw();
429// new TCanvas();
430// correctedEvents->Draw(); // FIXME: debug
9b0cb3c3 431
7a0032a9 432 return fAllEvents;
9b0cb3c3 433}
434
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
0d300b4f 451 const Int_t nHists = kNProcs*3+5;
9b0cb3c3 452 TList collections[nHists];
453
454 Int_t count = 0;
455 while ((obj = iter->Next())) {
456
457 AliCollisionNormalization* entry = dynamic_cast<AliCollisionNormalization*> (obj);
458 if (entry == 0)
459 continue;
460 Int_t ihist = -1;
461
462 for(Int_t iproc = 0; iproc < kNProcs; iproc++){
463 if (entry->fHistVzMCGen[iproc] ) collections[++ihist].Add(entry->fHistVzMCGen[iproc] );
464 if (entry->fHistVzMCRec[iproc] ) collections[++ihist].Add(entry->fHistVzMCRec[iproc] );
465 if (entry->fHistVzMCTrg[iproc] ) collections[++ihist].Add(entry->fHistVzMCTrg[iproc] );
466 }
467 if (entry->fHistVzData ) collections[++ihist].Add(entry->fHistVzData );
468 if (entry->fHistProcTypes ) collections[++ihist].Add(entry->fHistProcTypes );
469 if (entry->fHistStatBin0 ) collections[++ihist].Add(entry->fHistStatBin0 );
0d300b4f 470 if (entry->fHistStat ) collections[++ihist].Add(entry->fHistStat );
9b0cb3c3 471
472 count++;
473 }
474
475 Int_t ihist = -1;
476 for(Int_t iproc = 0; iproc < kNProcs; iproc++){
477 if (fHistVzMCGen[iproc] ) fHistVzMCGen[iproc] ->Merge(&collections[++ihist]);
478 if (fHistVzMCRec[iproc] ) fHistVzMCRec[iproc] ->Merge(&collections[++ihist]);
479 if (fHistVzMCTrg[iproc] ) fHistVzMCTrg[iproc] ->Merge(&collections[++ihist]);
480
481 }
482 if (fHistVzData ) fHistVzData ->Merge(&collections[++ihist]);
483 if (fHistProcTypes ) fHistProcTypes ->Merge(&collections[++ihist]);
484 if (fHistStatBin0 ) fHistStatBin0 ->Merge(&collections[++ihist]);
0d300b4f 485 if (fHistStat ) fHistStat ->Merge(&collections[++ihist]);
9b0cb3c3 486
487
488 delete iter;
489
490 return count+1;
491}
492
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
3aa0e5b3 518Int_t AliCollisionNormalization::GetProcessType(const AliMCEvent * mcEvt) {
9b0cb3c3 519
520 // Determine if the event was generated with pythia or phojet and return the process type
521
522 // Check if mcEvt is fine
523 if (!mcEvt) {
524 AliFatal("NULL mc event");
525 }
526
527 // Determine if it was a pythia or phojet header, and return the correct process type
528 AliGenPythiaEventHeader * headPy = 0;
529 AliGenDPMjetEventHeader * headPho = 0;
530 AliGenEventHeader * htmp = mcEvt->GenEventHeader();
531 if(!htmp) {
532 AliFatal("Cannot Get MC Header!!");
533 }
534 if( TString(htmp->IsA()->GetName()) == "AliGenPythiaEventHeader") {
535 headPy = (AliGenPythiaEventHeader*) htmp;
536 } else if (TString(htmp->IsA()->GetName()) == "AliGenDPMjetEventHeader") {
537 headPho = (AliGenDPMjetEventHeader*) htmp;
538 } else {
539 AliError("Unknown header");
540 }
541
542 // Determine process type
543 if(headPy) {
544 if(headPy->ProcessType() == 92 || headPy->ProcessType() == 93) {
545 // single difractive
546 return kProcSD;
547 } else if (headPy->ProcessType() == 94) {
548 // double diffractive
549 return kProcDD;
550 }
551 else if(headPy->ProcessType() != 92 && headPy->ProcessType() != 93 && headPy->ProcessType() != 94) {
552 // non difractive
553 return kProcND;
554 }
555 } else if (headPho) {
556 if(headPho->ProcessType() == 5 || headPho->ProcessType() == 6 ) {
557 // single difractive
558 return kProcSD;
559 } else if (headPho->ProcessType() == 7) {
560 // double diffractive
561 return kProcDD;
562 } else if(headPho->ProcessType() != 5 && headPho->ProcessType() != 6 && headPho->ProcessType() != 7 ) {
563 // non difractive
564 return kProcND;
565 }
566 }
567
568
569 // no process type found?
570 AliError(Form("Unknown header: %s", htmp->IsA()->GetName()));
571 return kProcUnknown;
572}
573
574Double_t AliCollisionNormalization::GetProcessWeight(Int_t proc) {
575
576 // Return a weight to be used when combining the MC histos to
0d300b4f 577 // compute efficiency, defined as the ratio XS in generator / XS
578 // measured
9b0cb3c3 579
3aa0e5b3 580 Float_t refSD, refDD, refND, errorSD, errorDD, errorND;
581 GetRelativeFractions(fReferenceXS,refSD, refDD, refND, errorSD, errorDD, errorND);
9b0cb3c3 582
583 static Double_t total = fHistProcTypes->Integral();
584 if (fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcUnknown)) > 0) {
585 AliError("There were unknown processes!!!");
586 }
3aa0e5b3 587 static Double_t lSD = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcSD))/total;
588 static Double_t lDD = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcDD))/total;
589 static Double_t lND = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcND))/total;
9b0cb3c3 590
591 if (fVerbose > 2) {
592 AliInfo(Form("Total MC evts: %f",total));
3aa0e5b3 593 AliInfo(Form(" Frac SD %4.4f", lSD));
594 AliInfo(Form(" Frac DD %4.4f", lDD));
595 AliInfo(Form(" Frac ND %4.4f", lND));
9b0cb3c3 596 }
597
598 switch(proc) {
599 case kProcSD:
3aa0e5b3 600 return refSD/lSD;
9b0cb3c3 601 break;
602 case kProcDD:
3aa0e5b3 603 return refDD/lDD;
9b0cb3c3 604 break;
605 case kProcND:
3aa0e5b3 606 return refND/lND;
9b0cb3c3 607 break;
608 default:
609 AliError("Unknown process");
610 }
611
612 return 0;
613
614}
615
3aa0e5b3 616void AliCollisionNormalization::GetRelativeFractions(Int_t origin, Float_t& refSD, Float_t& refDD, Float_t& refND, Float_t& errorSD, Float_t& errorDD, Float_t& errorND)
9b0cb3c3 617{
618 // Returns fraction of XS (SD, ND, DD) and corresponding error
619 // Stolen from Jan Fiete's drawSystematics macro
620
621 // origin:
622 // -1 = Pythia (test)
623 // 0 = UA5
624 // 1 = Data 1.8 TeV
625 // 2 = Tel-Aviv
626 // 3 = Durham
627 //
628
0d300b4f 629 Double_t epsilon = 0.0001;
630 if(TMath::Abs(fEnergy-900)<epsilon) {
631
632 switch (origin)
633 {
634 case -10: // Pythia default at 7 GeV, 50% error
635 AliInfo("PYTHIA x-sections");
3aa0e5b3 636 refSD = 0.192637; errorSD = refSD * 0.5;
637 refDD = 0.129877; errorDD = refDD * 0.5;
638 refND = 0.677486; errorND = 0;
0d300b4f 639 break;
640
641 case -1: // Pythia default at 900 GeV, as test
642 AliInfo("PYTHIA x-sections");
3aa0e5b3 643 refSD = 0.223788;
644 refDD = 0.123315;
645 refND = 0.652897;
9b0cb3c3 646 break;
647
0d300b4f 648 case 0: // UA5
649 AliInfo("UA5 x-sections a la first paper");
3aa0e5b3 650 refSD = 0.153; errorSD = 0.05;
651 refDD = 0.080; errorDD = 0.05;
652 refND = 0.767; errorND = 0;
0d300b4f 653 break;
654
655 case 10: // UA5
656 AliInfo("UA5 x-sections hadron level definition for Pythia");
657 // Fractions in Pythia with UA5 cuts selection for SD
658 // ND: 0.688662
659 // SD: 0.188588 --> this should be 15.3
660 // DD: 0.122750
3aa0e5b3 661 refSD = 0.224 * 0.153 / 0.189; errorSD = 0.023 * 0.224 / 0.189;
662 refDD = 0.095; errorDD = 0.06;
663 refND = 1.0 - refSD - refDD; errorND = 0;
0d300b4f 664 break;
665
666 case 11: // UA5
667 AliInfo("UA5 x-sections hadron level definition for Phojet");
668 // Fractions in Phojet with UA5 cuts selection for SD
669 // ND: 0.783573
670 // SD: 0.151601 --> this should be 15.3
671 // DD: 0.064827
3aa0e5b3 672 refSD = 0.191 * 0.153 / 0.152; errorSD = 0.023 * 0.191 / 0.152;
673 refDD = 0.095; errorDD = 0.06;
674 refND = 1.0 - refSD - refDD; errorND = 0;
0d300b4f 675 break;
676 case 2: // tel-aviv model
677 AliInfo("Tel-aviv model x-sections");
3aa0e5b3 678 refSD = 0.171;
679 refDD = 0.094;
680 refND = 1 - refSD - refDD;
0d300b4f 681 break;
682
683 case 3: // durham model
684 AliInfo("Durham model x-sections");
3aa0e5b3 685 refSD = 0.190;
686 refDD = 0.125;
687 refND = 1 - refSD - refDD;
9b0cb3c3 688 break;
0d300b4f 689 default:
690 AliFatal(Form("Unknown origin %d, Energy %f", origin, fEnergy));
691 }
692 }
693 else if(TMath::Abs(fEnergy-1800)<epsilon) {
694 switch (origin)
695 {
696 case 20: // E710, 1.8 TeV
9b0cb3c3 697 AliInfo("E710 x-sections hadron level definition for Pythia");
698 // ND: 0.705709
699 // SD: 0.166590 --> this should be 15.9
700 // DD: 0.127701
3aa0e5b3 701 refSD = 0.217 * 0.159 / 0.167; errorSD = 0.024 * 0.217 / 0.167;
702 refDD = 0.075 * 1.43; errorDD = 0.02 * 1.43;
703 refND = 1.0 - refSD - refDD; errorND = 0;
9b0cb3c3 704 break;
705
0d300b4f 706 case 21: // E710, 1.8 TeV
707 AliInfo("E710 x-sections hadron level definition for Phojet");
708 // ND: 0.817462
709 // SD: 0.125506 --> this should be 15.9
710 // DD: 0.057032
3aa0e5b3 711 refSD = 0.161 * 0.159 / 0.126; errorSD = 0.024 * 0.161 / 0.126;
712 refDD = 0.075 * 1.43; errorDD = 0.02 * 1.43;
713 refND = 1.0 - refSD - refDD; errorND = 0;
0d300b4f 714 break;
715
716 case 1: // data 1.8 TeV
717 AliInfo("??? x-sections");
3aa0e5b3 718 refSD = 0.152;
719 refDD = 0.092;
720 refND = 1 - refSD - refDD;
0d300b4f 721 break;
722 default:
723 AliFatal(Form("Unknown origin %d, Energy %f", origin, fEnergy));
724 }
725 }
726 else {
e8b839ab 727 AliFatal(Form("Unknown energy %f", fEnergy));
9b0cb3c3 728 }
0d300b4f 729
9b0cb3c3 730}
731
0d300b4f 732