DQM configure file
[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::fgkProcLabel[] = {"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")+ 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
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* 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
435 Long64_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
493 void 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
503 void 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 }
510 void 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
518 Int_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
574 Double_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
616 void 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