AliCollisionNormalization:
[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 {
41
42   // ctor
43
44   fNbinsVz =  30;
45   fMinVz   = -15;
46   fMaxVz   = +15;
47
48   for(Int_t iproc = 0; iproc < kNProcs; iproc++){
49       fHistVzMCGen[iproc] = 0;
50       fHistVzMCRec[iproc] = 0;
51       fHistVzMCTrg[iproc] = 0;
52   }
53   
54   
55   BookAllHistos();
56 }
57
58 AliCollisionNormalization::AliCollisionNormalization(Int_t nbinz, Float_t minz, Float_t maxz):
59   TObject(),
60   fNbinsVz(0),
61   fMinVz(0),
62   fMaxVz(0),
63   fZRange(9.99),
64   fIsMC(0),
65   fReferenceXS(0),
66   fVerbose(0),
67   fEnergy(900),
68   fHistVzData       (0),
69   fHistProcTypes    (0),
70   fHistStatBin0     (0),
71   fHistStat     (0)
72
73 {
74
75   // ctor, allows setting binning
76   fNbinsVz = nbinz;
77   fMinVz   = minz;
78   fMaxVz   = maxz;
79
80   for(Int_t iproc = 0; iproc < kNProcs; iproc++){
81       fHistVzMCGen[iproc] = 0;
82       fHistVzMCRec[iproc] = 0;
83       fHistVzMCTrg[iproc] = 0;
84   }
85
86   BookAllHistos();
87 }
88
89 AliCollisionNormalization::AliCollisionNormalization(const char * dataFile, const char * dataListName, 
90                                                      const char * mcFile,   const char * mcListName,
91                                                      const char * eventStatFile) :
92   TObject(),
93   fNbinsVz(0),
94   fMinVz(0),
95   fMaxVz(0),
96   fZRange(9.99),
97   fIsMC(0),
98   fReferenceXS(0),
99   fVerbose(0),
100   fEnergy(900),
101   fHistVzData       (0),
102   fHistProcTypes    (0),
103   fHistStatBin0     (0),
104   fHistStat     (0)
105
106 {
107
108   // ctor, loads histograms from file
109   for(Int_t iproc = 0; iproc < kNProcs; iproc++){
110       fHistVzMCGen[iproc] = 0;
111       fHistVzMCRec[iproc] = 0;
112       fHistVzMCTrg[iproc] = 0;
113   }
114
115
116   TFile * fdata = new TFile (dataFile);
117   TFile * fmc   = new TFile (mcFile  );
118   TFile * fstat = new TFile(eventStatFile);
119   
120   if (!fdata->IsOpen() || !fmc->IsOpen() || !fstat->IsOpen()) {
121     AliFatal("Cannot open input file(s)");
122   }
123   
124   TList * ldata = (TList*) fdata->Get(dataListName);
125   TList * lmc   = (TList*) fmc  ->Get(mcListName  );
126
127   AliCollisionNormalization * cndata = (AliCollisionNormalization*) ldata->FindObject("AliCollisionNormalization");
128   AliCollisionNormalization * cnmc   = (AliCollisionNormalization*) lmc  ->FindObject("AliCollisionNormalization");
129
130
131   // Assign or book all histos
132   for(Int_t iproc = 0; iproc < kNProcs; iproc++){
133     fHistVzMCGen[iproc]= cnmc->GetVzMCGen(iproc);
134     fHistVzMCRec[iproc]= cnmc->GetVzMCRec(iproc);
135     fHistVzMCTrg[iproc]= cnmc->GetVzMCTrg(iproc);
136   }
137   fHistVzData        = cndata->GetVzData();
138   fHistProcTypes     = cnmc->GetHistProcTypes();
139     
140   fHistStatBin0      =  (TH1F*) fstat->Get("fHistStatistics_Bin0");
141   fHistStat          =  (TH1F*) fstat->Get("fHistStatistics");
142   
143 }
144
145
146 AliCollisionNormalization::~AliCollisionNormalization(){
147
148   // dtor
149   for(Int_t iproc = 0; iproc < kNProcs; iproc++){
150     if(fHistVzMCGen[iproc]) { delete fHistVzMCGen[iproc]      ; fHistVzMCGen[iproc]      =0;}
151     if(fHistVzMCRec[iproc]) { delete fHistVzMCRec[iproc]      ; fHistVzMCRec[iproc]      =0;}
152     if(fHistVzMCTrg[iproc]) { delete fHistVzMCTrg[iproc]      ; fHistVzMCTrg[iproc]      =0;}    
153   }
154   
155   if(fHistVzData       ) { delete fHistVzData       ; fHistVzData       =0;}
156   if(fHistStatBin0     ) { delete fHistStatBin0     ; fHistStatBin0     =0;}
157   if(fHistStat         ) { delete fHistStat         ; fHistStat         =0;}
158   if(fHistProcTypes    ) { delete fHistProcTypes    ; fHistProcTypes    =0;}
159
160 }
161   
162 void AliCollisionNormalization::BookAllHistos(){
163   // books all histos
164   // Book histos of vz distributions vs multiplicity
165   // if vzOnly == kTRUE, it returns a 1d histo with vz dist only
166
167   // Do not attach histos to the current directory
168   Bool_t oldStatus = TH1::AddDirectoryStatus();
169   TH1::AddDirectory(kFALSE);
170
171   for(Int_t iproc = 0; iproc < kNProcs; iproc++){
172     fHistVzMCGen [iproc]      = (TH2F*) BookVzHisto(TString("fHistVzMCGen")+ fProcLabel[iproc]      ,"Vz distribution of generated events vs rec multiplicity    ");
173     fHistVzMCRec [iproc]      = (TH2F*) BookVzHisto(TString("fHistVzMCRec")+ fProcLabel[iproc]      ,"Vz distribution of reconstructed events vs rec multiplicity");
174     fHistVzMCTrg [iproc]      = (TH2F*) BookVzHisto(TString("fHistVzMCTrg")+ fProcLabel[iproc]      ,"Vz distribution of triggered events vs rec multiplicity    ");    
175   }
176   fHistVzData        = (TH2F*) BookVzHisto("fHistVzData"       ,"Vz distribution of triggered events vs rec multiplicity    ");
177   fHistProcTypes     = new TH1F   ("fHistProcTypes", "Number of events in the different process classes", kNProcs, -0.5 , kNProcs-0.5);
178
179   fHistProcTypes->GetXaxis()->SetBinLabel(kProcSD+1,"SD");
180   fHistProcTypes->GetXaxis()->SetBinLabel(kProcND+1,"ND");
181   fHistProcTypes->GetXaxis()->SetBinLabel(kProcDD+1,"DD");
182   fHistProcTypes->GetXaxis()->SetBinLabel(kProcUnknown+1,"Unknown");
183
184   TH1::AddDirectory(oldStatus);
185
186 }
187
188 TH1 * AliCollisionNormalization::BookVzHisto(const char * name , const char * title, Bool_t vzOnly) {
189
190   // Book histos of vz distributions vs multiplicity
191   // if vzOnly == kTRUE, it returns a 1d histo with vz dist only
192
193  
194   // Do not attach histos to the current directory
195   Bool_t oldStatus = TH1::AddDirectoryStatus();
196   TH1::AddDirectory(kFALSE);
197
198   TH1 * h;
199   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};
200   if (vzOnly) {
201     h = new TH1F(name,title,22,binLimitsVtx);
202   } else {
203     h = new TH2F(name,title,22,binLimitsVtx,100,Double_t(-0.5),Double_t(99.5));
204   }
205  
206   h->SetXTitle("V_{z} (cm)");
207   h->SetYTitle("n_{trk}");
208   h->Sumw2();
209
210   TH1::AddDirectory(oldStatus);
211
212   return h;
213
214 }
215
216 TH2F * AliCollisionNormalization::GetVzMCGen (Int_t procType) { 
217
218   //  returns MC gen histo. If proc type < 0 sums over all proc types, reweighting the XS
219
220   if(procType>=0)  return fHistVzMCGen[procType] ;
221
222   TH2F * sum = (TH2F*) fHistVzMCGen[kProcSD]->Clone();
223   sum->Reset();
224
225   for(Int_t iproc = 0; iproc < kProcUnknown; iproc++){
226     sum->Add(fHistVzMCGen[iproc],GetProcessWeight(iproc));    
227   }
228   
229   return sum;
230 }
231 TH2F *   AliCollisionNormalization::GetVzMCRec       (Int_t procType) { 
232
233   //  returns MC rec histo. If proc type < 0 sums over all proc types, reweighting the XS
234
235   if(procType>=0)  return fHistVzMCRec[procType] ;
236
237   TH2F * sum = (TH2F*) fHistVzMCRec[kProcSD]->Clone();
238   sum->Reset();
239
240   for(Int_t iproc = 0; iproc < kProcUnknown; iproc++){
241     sum->Add(fHistVzMCRec[iproc],GetProcessWeight(iproc));    
242   }
243   
244   return sum;
245
246 }
247
248
249 TH2F *   AliCollisionNormalization::GetVzMCTrg       (Int_t procType) { 
250
251   //  returns MC trg histo. If proc type < 0 sums over all proc types, reweighting the XS
252
253   if(procType>=0)  return fHistVzMCTrg[procType] ;
254
255   TH2F * sum = (TH2F*) fHistVzMCTrg[kProcSD]->Clone();
256   sum->Reset();
257
258   for(Int_t iproc = 0; iproc < kProcUnknown; iproc++){
259     sum->Add(fHistVzMCTrg[iproc],GetProcessWeight(iproc));    
260   }
261   
262   return sum;
263
264 }
265
266 Double_t AliCollisionNormalization::ComputeNint() {
267
268   // Compute the number of collisions applying all corrections
269   // TODO: check error propagation  
270
271   TH1 * hVzData = fHistVzData->ProjectionX("_px",2,-1,"E"); // Skip zero bin 
272   Int_t allEventsWithVertex = (Int_t) fHistVzData->Integral(0, fHistVzData->GetNbinsX()+1,
273                                                             2, fHistVzData->GetNbinsY()+1); // include under/overflow!, skip 0 bin
274
275   // Assign histos reweighted with measured XS
276   TH2F * histVzMCGen = GetVzMCGen(-1);
277   TH2F * histVzMCTrg = GetVzMCTrg(-1);
278   TH2F * histVzMCRec = GetVzMCRec(-1);
279
280   // Before we start: print number of input events to the physics
281   // selection: this allows to crosscheck that all runs were
282   // successfully processed (useful if running on grid: you may have a
283   // crash without noticing it).
284   AliInfo(Form("Input Events (No cuts: %d, After Phys. Sel.:%d)",
285                Int_t(fHistStat->GetBinContent(1,1)),
286                Int_t(fHistStat->GetBinContent(fHistStat->GetNbinsX(),1)))); // Fixme: will this change with a different trigger scheme?
287
288   // Get or compute BG. This assumes the CINT1B suite
289   Double_t triggeredEventsWith0MultWithBG = fHistVzData->Integral(0, fHistVzData->GetNbinsX()+1,1, 1);
290   Double_t bg = 0; // This will include beam gas + accidentals
291   if (fHistStatBin0->GetNbinsY() > 4) { // FIXME: we need a better criterion to decide...
292     AliInfo("Using BG computed by Physics Selection");
293     bg = fHistStatBin0->GetBinContent(fHistStatBin0->GetNbinsX(),AliPhysicsSelection::kStatRowBG);
294     bg += fHistStatBin0->GetBinContent(fHistStatBin0->GetNbinsX(),AliPhysicsSelection::kStatRowAcc);
295     Int_t cint1B = (Int_t) fHistStatBin0->GetBinContent(fHistStatBin0->GetNbinsX(),1);  
296     if (cint1B != Int_t(triggeredEventsWith0MultWithBG)) {
297       AliWarning(Form("Events in bin0 from physics selection and local counter not consistent: %d - %d", cint1B, Int_t(triggeredEventsWith0MultWithBG)));
298     }
299   } else {
300     AliInfo("Computing BG using CINT1A/B/C/E, ignoring intensities");
301     Int_t icol = fHistStatBin0->GetNbinsX();
302     Int_t cint1B = (Int_t) fHistStatBin0->GetBinContent(icol,1);        
303     Int_t cint1A = (Int_t) fHistStatBin0->GetBinContent(icol,2);        
304     Int_t cint1C = (Int_t) fHistStatBin0->GetBinContent(icol,3);        
305     Int_t cint1E = (Int_t) fHistStatBin0->GetBinContent(icol,4);      
306     bg   = cint1A + cint1C-2*cint1E ;
307     if (cint1B != triggeredEventsWith0MultWithBG) {
308       AliWarning(Form("Events in bin0 from physics selection and local counter not consistent: %d - %d", cint1B, triggeredEventsWith0MultWithBG));
309     }
310   }
311
312   Double_t triggeredEventsWith0Mult = triggeredEventsWith0MultWithBG - bg; 
313   if(fVerbose > 0) AliInfo(Form("Measured events with vertex: %d",allEventsWithVertex));
314   Double_t bin0 = fHistVzData->Integral(0, fHistVzData->GetNbinsX()+1,
315                                         1, 1);
316   if(fVerbose > 0) AliInfo(Form("Zero Bin, Meas: %2.2f, BG: %2.2f, Meas - BG: %2.2f", bin0, bg, triggeredEventsWith0Mult)); 
317
318   // This pointers are here so that I use the same names used in Jan Fiete's code (allows for faster/easier comparison)
319
320   TH2* eTrig =    histVzMCTrg; 
321   TH2* eTrigVtx = histVzMCRec; 
322   TH1* eTrigVtx_projx = eTrigVtx->ProjectionX("eTrigVtx_projx", 2, eTrigVtx->GetNbinsX()+1);
323
324   // compute trigger and vertex efficiency
325   TH2 * effVtxTrig = (TH2*) histVzMCRec->Clone("effVtxTrig");
326   effVtxTrig->Reset();
327   effVtxTrig->Divide(histVzMCRec,histVzMCGen,1,1,"B");
328   // Apply correction to data to get corrected events
329   TH2 * correctedEvents = (TH2*) fHistVzData->Clone("correctedEvents");
330   correctedEvents->Divide(effVtxTrig);
331
332   //  TH1 * correctedEvents = fHistVzData->ProjectionX("eTrigVtx_projx", 2, eTrigVtx->GetNbinsX()+1); 
333
334   // loop over vertex bins
335   for (Int_t i = 1; i <= fHistVzData->GetNbinsX(); i++)
336     {
337       Double_t alpha = (Double_t) hVzData->GetBinContent(i) / allEventsWithVertex;
338       Double_t events = alpha * triggeredEventsWith0Mult;
339       
340       if (histVzMCRec->GetBinContent(i,1) == 0)
341         continue;
342
343       Double_t fZ = eTrigVtx_projx->Integral(0, eTrigVtx_projx->GetNbinsX()+1) / eTrigVtx_projx->GetBinContent(i) *
344         eTrig->GetBinContent(i, 1) / eTrig->Integral(0, eTrig->GetNbinsX()+1, 1, 1);
345
346       events *= fZ;
347
348       // multiply with trigger correction
349       Double_t correction = histVzMCGen->GetBinContent(i,1)/histVzMCTrg->GetBinContent(i,1);
350       events *= correction;
351
352       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, 
353                                histVzMCRec->GetBinContent(i,1));
354       correctedEvents->SetBinContent(i, 1, events);
355     }
356
357   
358   
359   // Integrate correctedEvents over full z range
360   Double_t allEvents       = correctedEvents->Integral(0, correctedEvents->GetNbinsX()+1,0, correctedEvents->GetNbinsY()+1);
361   // Integrate correctedEvents over needed z range
362   Double_t allEventsZrange = correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange),
363                                                        0, correctedEvents->GetNbinsY()+1);
364   
365   if(fVerbose > 1) AliInfo(Form("Results in |Vz| < %3.3f",fZRange));
366   if(fVerbose > 1) AliInfo(Form(" Events in Bin0: %2.2f, With > 1 track: %2.2f, All corrected: %2.2f",
367                                 correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange),1,1),
368                                 correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange), 
369                                                           2,correctedEvents->GetNbinsX()+1),
370                                 allEventsZrange));
371   if(fVerbose > 1) {    
372     Int_t nbin = histVzMCRec->GetNbinsX();
373     AliInfo(Form("Efficiency in the zero bin: %3.3f", histVzMCRec->Integral(0,nbin+1,1,1)/histVzMCGen->Integral(0,nbin+1,1,1) ));
374   }
375   
376
377   AliInfo(Form("Number of collisions in full phase space: %2.2f", allEvents));
378 //   effVtxTrig->Draw();
379 //   new TCanvas();
380 //   correctedEvents->Draw(); // FIXME: debug
381
382   return allEvents;
383 }
384
385 Long64_t AliCollisionNormalization::Merge(TCollection* list)
386 {
387   // Merge a list of AliCollisionNormalization objects with this
388   // (needed for PROOF).  
389   // Returns the number of merged objects (including this).
390
391   if (!list)
392     return 0;
393
394   if (list->IsEmpty())
395     return 1;
396
397   TIterator* iter = list->MakeIterator();
398   TObject* obj;
399   
400   // collections of all histograms
401   const Int_t nHists = kNProcs*3+5;
402   TList collections[nHists];
403
404   Int_t count = 0;
405   while ((obj = iter->Next())) {
406
407     AliCollisionNormalization* entry = dynamic_cast<AliCollisionNormalization*> (obj);
408     if (entry == 0) 
409       continue;
410     Int_t ihist = -1;
411
412     for(Int_t iproc = 0; iproc < kNProcs; iproc++){
413       if (entry->fHistVzMCGen[iproc]      ) collections[++ihist].Add(entry->fHistVzMCGen[iproc]      );
414       if (entry->fHistVzMCRec[iproc]      ) collections[++ihist].Add(entry->fHistVzMCRec[iproc]      );
415       if (entry->fHistVzMCTrg[iproc]      ) collections[++ihist].Add(entry->fHistVzMCTrg[iproc]      );
416     }
417     if (entry->fHistVzData       ) collections[++ihist].Add(entry->fHistVzData       );
418     if (entry->fHistProcTypes    ) collections[++ihist].Add(entry->fHistProcTypes    );
419     if (entry->fHistStatBin0     ) collections[++ihist].Add(entry->fHistStatBin0     );
420     if (entry->fHistStat         ) collections[++ihist].Add(entry->fHistStat         );
421
422     count++;
423   }
424
425   Int_t ihist = -1;
426   for(Int_t iproc = 0; iproc < kNProcs; iproc++){
427     if (fHistVzMCGen[iproc]      ) fHistVzMCGen[iproc]      ->Merge(&collections[++ihist]);
428     if (fHistVzMCRec[iproc]      ) fHistVzMCRec[iproc]      ->Merge(&collections[++ihist]);
429     if (fHistVzMCTrg[iproc]      ) fHistVzMCTrg[iproc]      ->Merge(&collections[++ihist]);
430
431   }  
432   if (fHistVzData       ) fHistVzData       ->Merge(&collections[++ihist]);
433   if (fHistProcTypes    ) fHistProcTypes    ->Merge(&collections[++ihist]);
434   if (fHistStatBin0     ) fHistStatBin0     ->Merge(&collections[++ihist]);
435   if (fHistStat         ) fHistStat         ->Merge(&collections[++ihist]);
436     
437   
438   delete iter;
439
440   return count+1;
441 }
442
443 void AliCollisionNormalization::FillVzMCGen(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt) {
444
445   // Fill MC gen histo and the process types statistics
446
447   Int_t evtype = GetProcessType(mcEvt);
448   // When I fill gen histos, I also fill statistics of process types (used to detemine ratios afterwards).
449   fHistProcTypes->Fill(evtype);
450   fHistVzMCGen[evtype]->Fill(vz,ntrk);
451 }
452
453 void AliCollisionNormalization::FillVzMCRec(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt){
454
455   // Fill MC rec histo
456   Int_t evtype = GetProcessType(mcEvt);
457   fHistVzMCRec[evtype]->Fill(vz,ntrk);
458
459 }
460 void AliCollisionNormalization::FillVzMCTrg(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt) {      
461
462   // Fill MC trg histo
463   Int_t evtype = GetProcessType(mcEvt);
464   fHistVzMCTrg[evtype]->Fill(vz,ntrk);
465 }
466
467
468 Int_t AliCollisionNormalization::GetProcessType(AliMCEvent * mcEvt) {
469
470   // Determine if the event was generated with pythia or phojet and return the process type
471
472   // Check if mcEvt is fine
473   if (!mcEvt) {
474     AliFatal("NULL mc event");
475   } 
476
477   // Determine if it was a pythia or phojet header, and return the correct process type
478   AliGenPythiaEventHeader * headPy  = 0;
479   AliGenDPMjetEventHeader * headPho = 0;
480   AliGenEventHeader * htmp = mcEvt->GenEventHeader();
481   if(!htmp) {
482     AliFatal("Cannot Get MC Header!!");
483   }
484   if( TString(htmp->IsA()->GetName()) == "AliGenPythiaEventHeader") {
485     headPy =  (AliGenPythiaEventHeader*) htmp;
486   } else if (TString(htmp->IsA()->GetName()) == "AliGenDPMjetEventHeader") {
487     headPho = (AliGenDPMjetEventHeader*) htmp;
488   } else {
489     AliError("Unknown header");
490   }
491
492   // Determine process type
493   if(headPy)   {
494     if(headPy->ProcessType() == 92 || headPy->ProcessType() == 93) {
495       // single difractive
496       return kProcSD;
497     } else if (headPy->ProcessType() == 94) {
498       // double diffractive
499       return kProcDD;
500     }
501     else if(headPy->ProcessType() != 92 && headPy->ProcessType() != 93 && headPy->ProcessType() != 94) {    
502       // non difractive
503       return kProcND; 
504     }
505   } else if (headPho) {
506     if(headPho->ProcessType() == 5 || headPho->ProcessType() == 6 ) {
507       // single difractive
508       return kProcSD;
509     } else if (headPho->ProcessType() == 7) { 
510       // double diffractive
511       return kProcDD;      
512     } else if(headPho->ProcessType() != 5 && headPho->ProcessType() != 6  && headPho->ProcessType() != 7 ) {
513       // non difractive
514       return kProcND; 
515     }       
516   }
517   
518
519   // no process type found?
520   AliError(Form("Unknown header: %s", htmp->IsA()->GetName()));
521   return kProcUnknown;
522 }
523
524 Double_t AliCollisionNormalization::GetProcessWeight(Int_t proc) {
525
526   // Return a weight to be used when combining the MC histos to
527   // compute efficiency, defined as the ratio XS in generator / XS
528   // measured
529
530   Float_t ref_SD,  ref_DD,  ref_ND,  error_SD,  error_DD,  error_ND;
531   GetRelativeFractions(fReferenceXS,ref_SD,  ref_DD,  ref_ND,  error_SD,  error_DD,  error_ND);
532
533   static Double_t total = fHistProcTypes->Integral();
534   if (fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcUnknown)) > 0) {
535     AliError("There were unknown processes!!!");
536   }
537   static Double_t SD = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcSD))/total;
538   static Double_t DD = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcDD))/total;
539   static Double_t ND = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcND))/total;
540
541   if (fVerbose > 2) {
542     AliInfo(Form("Total MC evts: %f",total));
543     AliInfo(Form(" Frac SD %4.4f", SD));
544     AliInfo(Form(" Frac DD %4.4f", DD));
545     AliInfo(Form(" Frac ND %4.4f", ND));
546   }
547   
548   switch(proc) {
549   case kProcSD:
550     return ref_SD/SD;
551     break;
552   case kProcDD:
553     return ref_DD/DD;
554     break;
555   case kProcND:
556     return ref_ND/ND;
557     break;
558   default:
559     AliError("Unknown process");
560   }
561   
562   return 0;
563
564
565
566 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)
567 {
568   // Returns fraction of XS (SD, ND, DD) and corresponding error
569   // Stolen from Jan Fiete's drawSystematics macro
570
571   // origin: 
572   //   -1 = Pythia (test)
573   //   0 = UA5
574   //   1 = Data 1.8 TeV
575   //   2 = Tel-Aviv
576   //   3 = Durham
577   //
578
579   Double_t epsilon = 0.0001;
580   if(TMath::Abs(fEnergy-900)<epsilon) {
581
582     switch (origin)
583       {
584       case -10: // Pythia default at 7 GeV, 50% error
585         AliInfo("PYTHIA x-sections");
586         ref_SD = 0.192637; error_SD = ref_SD * 0.5;
587         ref_DD = 0.129877; error_DD = ref_DD * 0.5;
588         ref_ND = 0.677486; error_ND = 0;
589         break;
590         
591       case -1: // Pythia default at 900 GeV, as test
592         AliInfo("PYTHIA x-sections");
593       ref_SD = 0.223788;
594       ref_DD = 0.123315;
595       ref_ND = 0.652897;
596       break;
597       
598       case 0: // UA5
599         AliInfo("UA5 x-sections a la first paper");
600         ref_SD = 0.153; error_SD = 0.05;
601         ref_DD = 0.080; error_DD = 0.05;
602         ref_ND = 0.767; error_ND = 0;
603         break;
604         
605       case 10: // UA5
606         AliInfo("UA5 x-sections hadron level definition for Pythia"); 
607         // Fractions in Pythia with UA5 cuts selection for SD
608         // ND: 0.688662
609         // SD: 0.188588 --> this should be 15.3
610         // DD: 0.122750
611         ref_SD = 0.224 * 0.153 / 0.189; error_SD = 0.023 * 0.224 / 0.189;
612         ref_DD = 0.095;                 error_DD = 0.06; 
613         ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
614         break;
615         
616       case 11: // UA5
617         AliInfo("UA5 x-sections hadron level definition for Phojet"); 
618         // Fractions in Phojet with UA5 cuts selection for SD
619         // ND: 0.783573
620         // SD: 0.151601 --> this should be 15.3
621         // DD: 0.064827
622         ref_SD = 0.191 * 0.153 / 0.152; error_SD = 0.023 * 0.191 / 0.152;
623         ref_DD = 0.095;                 error_DD = 0.06; 
624         ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
625         break;
626       case 2: // tel-aviv model
627         AliInfo("Tel-aviv model x-sections");
628         ref_SD = 0.171;
629         ref_DD = 0.094;
630         ref_ND = 1 - ref_SD - ref_DD;
631         break;
632         
633       case 3: // durham model
634         AliInfo("Durham model x-sections");
635         ref_SD = 0.190;
636         ref_DD = 0.125;
637         ref_ND = 1 - ref_SD - ref_DD;
638       break;
639       default:
640         AliFatal(Form("Unknown origin %d, Energy %f", origin, fEnergy));
641       }
642   }
643   else if(TMath::Abs(fEnergy-1800)<epsilon) {     
644     switch (origin)
645       {
646       case 20: // E710, 1.8 TeV
647       AliInfo("E710 x-sections hadron level definition for Pythia");
648       // ND: 0.705709
649       // SD: 0.166590 --> this should be 15.9
650       // DD: 0.127701
651       ref_SD = 0.217 * 0.159 / 0.167; error_SD = 0.024 * 0.217 / 0.167;
652       ref_DD = 0.075 * 1.43;          error_DD = 0.02 * 1.43; 
653       ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
654       break;
655     
656       case 21: // E710, 1.8 TeV
657         AliInfo("E710 x-sections hadron level definition for Phojet"); 
658         // ND: 0.817462
659         // SD: 0.125506 --> this should be 15.9
660         // DD: 0.057032
661         ref_SD = 0.161 * 0.159 / 0.126; error_SD = 0.024 * 0.161 / 0.126;
662         ref_DD = 0.075 * 1.43;         error_DD = 0.02 * 1.43;
663         ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
664         break;
665         
666       case 1: // data 1.8 TeV
667         AliInfo("??? x-sections");
668         ref_SD = 0.152;
669         ref_DD = 0.092;
670         ref_ND = 1 - ref_SD - ref_DD;
671         break;
672       default:
673         AliFatal(Form("Unknown origin %d, Energy %f", origin, fEnergy));
674       }    
675   } 
676   else {
677     AliFatal(Form("Unknown energy %f", origin, fEnergy));
678   }
679     
680 }
681
682