Added getters and members for some event counters (A. Marin)
[u/mrichter/AliRoot.git] / ANALYSIS / AliCollisionNormalization.cxx
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
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),
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
66 AliCollisionNormalization::AliCollisionNormalization(Int_t nbinz, Float_t minz, Float_t maxz):
67   TObject(),
68   fNbinsVz(0),
69   fMinVz(0),
70   fMaxVz(0),
71   fZRange(9.99),
72   fIsMC(0),
73   fReferenceXS(0),
74   fVerbose(0),
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
104 AliCollisionNormalization::AliCollisionNormalization(const char * dataFile, const char * dataListName, 
105                                                      const char * mcFile,   const char * mcListName,
106                                                      const char * eventStatFile) :
107   TObject(),
108   fNbinsVz(0),
109   fMinVz(0),
110   fMaxVz(0),
111   fZRange(9.99),
112   fIsMC(0),
113   fReferenceXS(0),
114   fVerbose(0),
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
169 AliCollisionNormalization::~AliCollisionNormalization(){
170
171   // dtor
172   for(Int_t iproc = 0; iproc < kNProcs; iproc++){
173     if(fHistVzMCGen[iproc]) { delete fHistVzMCGen[iproc]      ; fHistVzMCGen[iproc]      =0;}
174     if(fHistVzMCRec[iproc]) { delete fHistVzMCRec[iproc]      ; fHistVzMCRec[iproc]      =0;}
175     if(fHistVzMCTrg[iproc]) { delete fHistVzMCTrg[iproc]      ; fHistVzMCTrg[iproc]      =0;}    
176   }
177   
178   if(fHistVzData       ) { delete fHistVzData       ; fHistVzData       =0;}
179   if(fHistStatBin0     ) { delete fHistStatBin0     ; fHistStatBin0     =0;}
180   if(fHistStat         ) { delete fHistStat         ; fHistStat         =0;}
181   if(fHistProcTypes    ) { delete fHistProcTypes    ; fHistProcTypes    =0;}
182
183 }
184   
185 void AliCollisionNormalization::BookAllHistos(){
186   // books all histos
187   // Book histos of vz distributions vs multiplicity
188   // if vzOnly == kTRUE, it returns a 1d histo with vz dist only
189
190   // Do not attach histos to the current directory
191   Bool_t oldStatus = TH1::AddDirectoryStatus();
192   TH1::AddDirectory(kFALSE);
193
194   for(Int_t iproc = 0; iproc < kNProcs; iproc++){
195     fHistVzMCGen [iproc]      = (TH2F*) BookVzHisto(TString("fHistVzMCGen")+ fProcLabel[iproc]      ,"Vz distribution of generated events vs rec multiplicity    ");
196     fHistVzMCRec [iproc]      = (TH2F*) BookVzHisto(TString("fHistVzMCRec")+ fProcLabel[iproc]      ,"Vz distribution of reconstructed events vs rec multiplicity");
197     fHistVzMCTrg [iproc]      = (TH2F*) BookVzHisto(TString("fHistVzMCTrg")+ fProcLabel[iproc]      ,"Vz distribution of triggered events vs rec multiplicity    ");    
198   }
199   fHistVzData        = (TH2F*) BookVzHisto("fHistVzData"       ,"Vz distribution of triggered events vs rec multiplicity    ");
200   fHistProcTypes     = new TH1F   ("fHistProcTypes", "Number of events in the different process classes", kNProcs, -0.5 , kNProcs-0.5);
201
202   fHistProcTypes->GetXaxis()->SetBinLabel(kProcSD+1,"SD");
203   fHistProcTypes->GetXaxis()->SetBinLabel(kProcND+1,"ND");
204   fHistProcTypes->GetXaxis()->SetBinLabel(kProcDD+1,"DD");
205   fHistProcTypes->GetXaxis()->SetBinLabel(kProcUnknown+1,"Unknown");
206
207   TH1::AddDirectory(oldStatus);
208
209 }
210
211 TH1 * AliCollisionNormalization::BookVzHisto(const char * name , const char * title, Bool_t vzOnly) {
212
213   // Book histos of vz distributions vs multiplicity
214   // if vzOnly == kTRUE, it returns a 1d histo with vz dist only
215
216  
217   // Do not attach histos to the current directory
218   Bool_t oldStatus = TH1::AddDirectoryStatus();
219   TH1::AddDirectory(kFALSE);
220
221   TH1 * h;
222   Double_t binLimitsVtx[] = {-30,-25,-20,-15,-10,-7,-5.5,-4,-3,-2,-1,0,1,2,3,4,5.5,7,10,15,20,25,30};
223   if (vzOnly) {
224     h = new TH1F(name,title,22,binLimitsVtx);
225   } else {
226     h = new TH2F(name,title,22,binLimitsVtx,100,Double_t(-0.5),Double_t(99.5));
227   }
228  
229   h->SetXTitle("V_{z} (cm)");
230   h->SetYTitle("n_{trk}");
231   h->Sumw2();
232
233   TH1::AddDirectory(oldStatus);
234
235   return h;
236
237 }
238
239 TH2F * AliCollisionNormalization::GetVzMCGen (Int_t procType) { 
240
241   //  returns MC gen histo. If proc type < 0 sums over all proc types, reweighting the XS
242
243   if(procType>=0)  return fHistVzMCGen[procType] ;
244
245   TH2F * sum = (TH2F*) fHistVzMCGen[kProcSD]->Clone();
246   sum->Reset();
247
248   for(Int_t iproc = 0; iproc < kProcUnknown; iproc++){
249     sum->Add(fHistVzMCGen[iproc],GetProcessWeight(iproc));    
250   }
251   
252   return sum;
253 }
254 TH2F *   AliCollisionNormalization::GetVzMCRec       (Int_t procType) { 
255
256   //  returns MC rec histo. If proc type < 0 sums over all proc types, reweighting the XS
257
258   if(procType>=0)  return fHistVzMCRec[procType] ;
259
260   TH2F * sum = (TH2F*) fHistVzMCRec[kProcSD]->Clone();
261   sum->Reset();
262
263   for(Int_t iproc = 0; iproc < kProcUnknown; iproc++){
264     sum->Add(fHistVzMCRec[iproc],GetProcessWeight(iproc));    
265   }
266   
267   return sum;
268
269 }
270
271
272 TH2F *   AliCollisionNormalization::GetVzMCTrg       (Int_t procType) { 
273
274   //  returns MC trg histo. If proc type < 0 sums over all proc types, reweighting the XS
275
276   if(procType>=0)  return fHistVzMCTrg[procType] ;
277
278   TH2F * sum = (TH2F*) fHistVzMCTrg[kProcSD]->Clone();
279   sum->Reset();
280
281   for(Int_t iproc = 0; iproc < kProcUnknown; iproc++){
282     sum->Add(fHistVzMCTrg[iproc],GetProcessWeight(iproc));    
283   }
284   
285   return sum;
286
287 }
288
289 Double_t AliCollisionNormalization::ComputeNint() {
290
291   // Compute the number of collisions applying all corrections
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* 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       
383       // if (histVzMCRec->GetBinContent(i,1) == 0)
384       //   continue;
385       if (!eTrigVtx_projx->GetBinContent(i) || !eTrig->Integral(0, eTrig->GetNbinsX()+1, 1, 1))
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, 
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),
408                                                        0, correctedEvents->GetNbinsY()+1);
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
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
421   if(fVerbose > 1) {    
422     Int_t nbin = histVzMCTrg->GetNbinsX();
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 ));
425   }
426   
427
428   AliInfo(Form("Number of collisions in full phase space: %2.2f", fAllEvents));
429 //   effVtxTrig->Draw();
430 //   new TCanvas();
431 //   correctedEvents->Draw(); // FIXME: debug
432
433   return fAllEvents;
434 }
435
436 Long64_t AliCollisionNormalization::Merge(TCollection* list)
437 {
438   // Merge a list of AliCollisionNormalization objects with this
439   // (needed for PROOF).  
440   // Returns the number of merged objects (including this).
441
442   if (!list)
443     return 0;
444
445   if (list->IsEmpty())
446     return 1;
447
448   TIterator* iter = list->MakeIterator();
449   TObject* obj;
450   
451   // collections of all histograms
452   const Int_t nHists = kNProcs*3+5;
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     );
471     if (entry->fHistStat         ) collections[++ihist].Add(entry->fHistStat         );
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]);
486   if (fHistStat         ) fHistStat         ->Merge(&collections[++ihist]);
487     
488   
489   delete iter;
490
491   return count+1;
492 }
493
494 void AliCollisionNormalization::FillVzMCGen(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt) {
495
496   // Fill MC gen histo and the process types statistics
497
498   Int_t evtype = GetProcessType(mcEvt);
499   // When I fill gen histos, I also fill statistics of process types (used to detemine ratios afterwards).
500   fHistProcTypes->Fill(evtype);
501   fHistVzMCGen[evtype]->Fill(vz,ntrk);
502 }
503
504 void AliCollisionNormalization::FillVzMCRec(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt){
505
506   // Fill MC rec histo
507   Int_t evtype = GetProcessType(mcEvt);
508   fHistVzMCRec[evtype]->Fill(vz,ntrk);
509
510 }
511 void AliCollisionNormalization::FillVzMCTrg(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt) {      
512
513   // Fill MC trg histo
514   Int_t evtype = GetProcessType(mcEvt);
515   fHistVzMCTrg[evtype]->Fill(vz,ntrk);
516 }
517
518
519 Int_t AliCollisionNormalization::GetProcessType(AliMCEvent * mcEvt) {
520
521   // Determine if the event was generated with pythia or phojet and return the process type
522
523   // Check if mcEvt is fine
524   if (!mcEvt) {
525     AliFatal("NULL mc event");
526   } 
527
528   // Determine if it was a pythia or phojet header, and return the correct process type
529   AliGenPythiaEventHeader * headPy  = 0;
530   AliGenDPMjetEventHeader * headPho = 0;
531   AliGenEventHeader * htmp = mcEvt->GenEventHeader();
532   if(!htmp) {
533     AliFatal("Cannot Get MC Header!!");
534   }
535   if( TString(htmp->IsA()->GetName()) == "AliGenPythiaEventHeader") {
536     headPy =  (AliGenPythiaEventHeader*) htmp;
537   } else if (TString(htmp->IsA()->GetName()) == "AliGenDPMjetEventHeader") {
538     headPho = (AliGenDPMjetEventHeader*) htmp;
539   } else {
540     AliError("Unknown header");
541   }
542
543   // Determine process type
544   if(headPy)   {
545     if(headPy->ProcessType() == 92 || headPy->ProcessType() == 93) {
546       // single difractive
547       return kProcSD;
548     } else if (headPy->ProcessType() == 94) {
549       // double diffractive
550       return kProcDD;
551     }
552     else if(headPy->ProcessType() != 92 && headPy->ProcessType() != 93 && headPy->ProcessType() != 94) {    
553       // non difractive
554       return kProcND; 
555     }
556   } else if (headPho) {
557     if(headPho->ProcessType() == 5 || headPho->ProcessType() == 6 ) {
558       // single difractive
559       return kProcSD;
560     } else if (headPho->ProcessType() == 7) { 
561       // double diffractive
562       return kProcDD;      
563     } else if(headPho->ProcessType() != 5 && headPho->ProcessType() != 6  && headPho->ProcessType() != 7 ) {
564       // non difractive
565       return kProcND; 
566     }       
567   }
568   
569
570   // no process type found?
571   AliError(Form("Unknown header: %s", htmp->IsA()->GetName()));
572   return kProcUnknown;
573 }
574
575 Double_t AliCollisionNormalization::GetProcessWeight(Int_t proc) {
576
577   // Return a weight to be used when combining the MC histos to
578   // compute efficiency, defined as the ratio XS in generator / XS
579   // measured
580
581   Float_t ref_SD,  ref_DD,  ref_ND,  error_SD,  error_DD,  error_ND;
582   GetRelativeFractions(fReferenceXS,ref_SD,  ref_DD,  ref_ND,  error_SD,  error_DD,  error_ND);
583
584   static Double_t total = fHistProcTypes->Integral();
585   if (fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcUnknown)) > 0) {
586     AliError("There were unknown processes!!!");
587   }
588   static Double_t SD = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcSD))/total;
589   static Double_t DD = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcDD))/total;
590   static Double_t ND = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcND))/total;
591
592   if (fVerbose > 2) {
593     AliInfo(Form("Total MC evts: %f",total));
594     AliInfo(Form(" Frac SD %4.4f", SD));
595     AliInfo(Form(" Frac DD %4.4f", DD));
596     AliInfo(Form(" Frac ND %4.4f", ND));
597   }
598   
599   switch(proc) {
600   case kProcSD:
601     return ref_SD/SD;
602     break;
603   case kProcDD:
604     return ref_DD/DD;
605     break;
606   case kProcND:
607     return ref_ND/ND;
608     break;
609   default:
610     AliError("Unknown process");
611   }
612   
613   return 0;
614
615
616
617 void AliCollisionNormalization::GetRelativeFractions(Int_t origin, Float_t& ref_SD, Float_t& ref_DD, Float_t& ref_ND, Float_t& error_SD, Float_t& error_DD, Float_t& error_ND)
618 {
619   // Returns fraction of XS (SD, ND, DD) and corresponding error
620   // Stolen from Jan Fiete's drawSystematics macro
621
622   // origin: 
623   //   -1 = Pythia (test)
624   //   0 = UA5
625   //   1 = Data 1.8 TeV
626   //   2 = Tel-Aviv
627   //   3 = Durham
628   //
629
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");
644       ref_SD = 0.223788;
645       ref_DD = 0.123315;
646       ref_ND = 0.652897;
647       break;
648       
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;
689       break;
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
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     
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 {
728     AliFatal(Form("Unknown energy %f", fEnergy));
729   }
730     
731 }
732
733