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