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