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