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   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")+ fProcLabel[iproc]      ,"Vz distribution of generated events vs rec multiplicity    ");
 196     fHistVzMCRec [iproc]      = (TH2F*) BookVzHisto(TString("fHistVzMCRec")+ fProcLabel[iproc]      ,"Vz distribution of reconstructed events vs rec multiplicity");
 197     fHistVzMCTrg [iproc]      = (TH2F*) BookVzHisto(TString("fHistVzMCTrg")+ fProcLabel[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* eTrigVtx_projx = eTrigVtx->ProjectionX("eTrigVtx_projx", 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("eTrigVtx_projx", 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 (!eTrigVtx_projx->GetBinContent(i) || !eTrig->Integral(0, eTrig->GetNbinsX()+1, 1, 1))
 386         continue;
 387
 388       Double_t fZ = eTrigVtx_projx->Integral(0, eTrigVtx_projx->GetNbinsX()+1) / eTrigVtx_projx->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),
 408                                                        0, correctedEvents->GetNbinsY()+1);
 409
 410   fAllEventsZRangeMult1 = correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange),
 411                                                     2,correctedEvents->GetNbinsX()+1);
 412
 413   fAllEventsInBin0ZRange = correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange),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
 421   if(fVerbose > 1) {
 422     Int_t nbin = histVzMCTrg->GetNbinsX();
 423     fTrigEffBin0 =  histVzMCTrg->Integral(0,nbin+1,1,1)/histVzMCGen->Integral(0,nbin+1,1,1);
 424     AliInfo(Form("Trigger Efficiency in the zero bin: %3.3f", fTrigEffBin0 ));
 425   }
 426
 427
 428   AliInfo(Form("Number of collisions in full phase space: %2.2f", fAllEvents));
 429 //   effVtxTrig->Draw();
 430 //   new TCanvas();
 431 //   correctedEvents->Draw(); // FIXME: debug
 432
 433   return fAllEvents;
 434 }
 435
 436 Long64_t AliCollisionNormalization::Merge(TCollection* list)
 437 {
 438   // Merge a list of AliCollisionNormalization objects with this
 439   // (needed for PROOF).
 440   // Returns the number of merged objects (including this).
 441
 442   if (!list)
 443     return 0;
 444
 445   if (list->IsEmpty())
 446     return 1;
 447
 448   TIterator* iter = list->MakeIterator();
 449   TObject* obj;
 450
 451   // collections of all histograms
 452   const Int_t nHists = kNProcs*3+5;
 453   TList collections[nHists];
 454
 455   Int_t count = 0;
 456   while ((obj = iter->Next())) {
 457
 458     AliCollisionNormalization* entry = dynamic_cast<AliCollisionNormalization*> (obj);
 459     if (entry == 0)
 460       continue;
 461     Int_t ihist = -1;
 462
 463     for(Int_t iproc = 0; iproc < kNProcs; iproc++){
 464       if (entry->fHistVzMCGen[iproc]      ) collections[++ihist].Add(entry->fHistVzMCGen[iproc]      );
 465       if (entry->fHistVzMCRec[iproc]      ) collections[++ihist].Add(entry->fHistVzMCRec[iproc]      );
 466       if (entry->fHistVzMCTrg[iproc]      ) collections[++ihist].Add(entry->fHistVzMCTrg[iproc]      );
 467     }
 468     if (entry->fHistVzData       ) collections[++ihist].Add(entry->fHistVzData       );
 469     if (entry->fHistProcTypes    ) collections[++ihist].Add(entry->fHistProcTypes    );
 470     if (entry->fHistStatBin0     ) collections[++ihist].Add(entry->fHistStatBin0     );
 471     if (entry->fHistStat         ) collections[++ihist].Add(entry->fHistStat         );
 472
 473     count++;
 474   }
 475
 476   Int_t ihist = -1;
 477   for(Int_t iproc = 0; iproc < kNProcs; iproc++){
 478     if (fHistVzMCGen[iproc]      ) fHistVzMCGen[iproc]      ->Merge(&collections[++ihist]);
 479     if (fHistVzMCRec[iproc]      ) fHistVzMCRec[iproc]      ->Merge(&collections[++ihist]);
 480     if (fHistVzMCTrg[iproc]      ) fHistVzMCTrg[iproc]      ->Merge(&collections[++ihist]);
 481
 482   }
 483   if (fHistVzData       ) fHistVzData       ->Merge(&collections[++ihist]);
 484   if (fHistProcTypes    ) fHistProcTypes    ->Merge(&collections[++ihist]);
 485   if (fHistStatBin0     ) fHistStatBin0     ->Merge(&collections[++ihist]);
 486   if (fHistStat         ) fHistStat         ->Merge(&collections[++ihist]);
 487
 488
 489   delete iter;
 490
 491   return count+1;
 492 }
 493
 494 void AliCollisionNormalization::FillVzMCGen(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt) {
 495
 496   // Fill MC gen histo and the process types statistics
 497
 498   Int_t evtype = GetProcessType(mcEvt);
 499   // When I fill gen histos, I also fill statistics of process types (used to detemine ratios afterwards).
 500   fHistProcTypes->Fill(evtype);
 501   fHistVzMCGen[evtype]->Fill(vz,ntrk);
 502 }
 503
 504 void AliCollisionNormalization::FillVzMCRec(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt){
 505
 506   // Fill MC rec histo
 507   Int_t evtype = GetProcessType(mcEvt);
 508   fHistVzMCRec[evtype]->Fill(vz,ntrk);
 509
 510 }
 511 void AliCollisionNormalization::FillVzMCTrg(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt) {
 512
 513   // Fill MC trg histo
 514   Int_t evtype = GetProcessType(mcEvt);
 515   fHistVzMCTrg[evtype]->Fill(vz,ntrk);
 516 }
 517
 518
 519 Int_t AliCollisionNormalization::GetProcessType(AliMCEvent * mcEvt) {
 520
 521   // Determine if the event was generated with pythia or phojet and return the process type
 522
 523   // Check if mcEvt is fine
 524   if (!mcEvt) {
 525     AliFatal("NULL mc event");
 526   }
 527
 528   // Determine if it was a pythia or phojet header, and return the correct process type
 529   AliGenPythiaEventHeader * headPy  = 0;
 530   AliGenDPMjetEventHeader * headPho = 0;
 531   AliGenEventHeader * htmp = mcEvt->GenEventHeader();
 532   if(!htmp) {
 533     AliFatal("Cannot Get MC Header!!");
 534   }
 535   if( TString(htmp->IsA()->GetName()) == "AliGenPythiaEventHeader") {
 536     headPy =  (AliGenPythiaEventHeader*) htmp;
 537   } else if (TString(htmp->IsA()->GetName()) == "AliGenDPMjetEventHeader") {
 538     headPho = (AliGenDPMjetEventHeader*) htmp;
 539   } else {
 540     AliError("Unknown header");
 541   }
 542
 543   // Determine process type
 544   if(headPy)   {
 545     if(headPy->ProcessType() == 92 || headPy->ProcessType() == 93) {
 546       // single difractive
 547       return kProcSD;
 548     } else if (headPy->ProcessType() == 94) {
 549       // double diffractive
 550       return kProcDD;
 551     }
 552     else if(headPy->ProcessType() != 92 && headPy->ProcessType() != 93 && headPy->ProcessType() != 94) {
 553       // non difractive
 554       return kProcND;
 555     }
 556   } else if (headPho) {
 557     if(headPho->ProcessType() == 5 || headPho->ProcessType() == 6 ) {
 558       // single difractive
 559       return kProcSD;
 560     } else if (headPho->ProcessType() == 7) {
 561       // double diffractive
 562       return kProcDD;
 563     } else if(headPho->ProcessType() != 5 && headPho->ProcessType() != 6  && headPho->ProcessType() != 7 ) {
 564       // non difractive
 565       return kProcND;
 566     }
 567   }
 568
 569
 570   // no process type found?
 571   AliError(Form("Unknown header: %s", htmp->IsA()->GetName()));
 572   return kProcUnknown;
 573 }
 574
 575 Double_t AliCollisionNormalization::GetProcessWeight(Int_t proc) {
 576
 577   // Return a weight to be used when combining the MC histos to
 578   // compute efficiency, defined as the ratio XS in generator / XS
 579   // measured
 580
 581   Float_t ref_SD,  ref_DD,  ref_ND,  error_SD,  error_DD,  error_ND;
 582   GetRelativeFractions(fReferenceXS,ref_SD,  ref_DD,  ref_ND,  error_SD,  error_DD,  error_ND);
 583
 584   static Double_t total = fHistProcTypes->Integral();
 585   if (fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcUnknown)) > 0) {
 586     AliError("There were unknown processes!!!");
 587   }
 588   static Double_t SD = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcSD))/total;
 589   static Double_t DD = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcDD))/total;
 590   static Double_t ND = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcND))/total;
 591
 592   if (fVerbose > 2) {
 593     AliInfo(Form("Total MC evts: %f",total));
 594     AliInfo(Form(" Frac SD %4.4f", SD));
 595     AliInfo(Form(" Frac DD %4.4f", DD));
 596     AliInfo(Form(" Frac ND %4.4f", ND));
 597   }
 598
 599   switch(proc) {
 600   case kProcSD:
 601     return ref_SD/SD;
 602     break;
 603   case kProcDD:
 604     return ref_DD/DD;
 605     break;
 606   case kProcND:
 607     return ref_ND/ND;
 608     break;
 609   default:
 610     AliError("Unknown process");
 611   }
 612
 613   return 0;
 614
 615 }
 616
 617 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)
 618 {
 619   // Returns fraction of XS (SD, ND, DD) and corresponding error
 620   // Stolen from Jan Fiete's drawSystematics macro
 621
 622   // origin:
 623   //   -1 = Pythia (test)
 624   //   0 = UA5
 625   //   1 = Data 1.8 TeV
 626   //   2 = Tel-Aviv
 627   //   3 = Durham
 628   //
 629
 630   Double_t epsilon = 0.0001;
 631   if(TMath::Abs(fEnergy-900)<epsilon) {
 632
 633     switch (origin)
 634       {
 635       case -10: // Pythia default at 7 GeV, 50% error
 636         AliInfo("PYTHIA x-sections");
 637         ref_SD = 0.192637; error_SD = ref_SD * 0.5;
 638         ref_DD = 0.129877; error_DD = ref_DD * 0.5;
 639         ref_ND = 0.677486; error_ND = 0;
 640         break;
 641
 642       case -1: // Pythia default at 900 GeV, as test
 643         AliInfo("PYTHIA x-sections");
 644       ref_SD = 0.223788;
 645       ref_DD = 0.123315;
 646       ref_ND = 0.652897;
 647       break;
 648
 649       case 0: // UA5
 650         AliInfo("UA5 x-sections a la first paper");
 651         ref_SD = 0.153; error_SD = 0.05;
 652         ref_DD = 0.080; error_DD = 0.05;
 653         ref_ND = 0.767; error_ND = 0;
 654         break;
 655
 656       case 10: // UA5
 657         AliInfo("UA5 x-sections hadron level definition for Pythia");
 658         // Fractions in Pythia with UA5 cuts selection for SD
 659         // ND: 0.688662
 660         // SD: 0.188588 --> this should be 15.3
 661         // DD: 0.122750
 662         ref_SD = 0.224 * 0.153 / 0.189; error_SD = 0.023 * 0.224 / 0.189;
 663         ref_DD = 0.095;                 error_DD = 0.06;
 664         ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
 665         break;
 666
 667       case 11: // UA5
 668         AliInfo("UA5 x-sections hadron level definition for Phojet");
 669         // Fractions in Phojet with UA5 cuts selection for SD
 670         // ND: 0.783573
 671         // SD: 0.151601 --> this should be 15.3
 672         // DD: 0.064827
 673         ref_SD = 0.191 * 0.153 / 0.152; error_SD = 0.023 * 0.191 / 0.152;
 674         ref_DD = 0.095;                 error_DD = 0.06;
 675         ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
 676         break;
 677       case 2: // tel-aviv model
 678         AliInfo("Tel-aviv model x-sections");
 679         ref_SD = 0.171;
 680         ref_DD = 0.094;
 681         ref_ND = 1 - ref_SD - ref_DD;
 682         break;
 683
 684       case 3: // durham model
 685         AliInfo("Durham model x-sections");
 686         ref_SD = 0.190;
 687         ref_DD = 0.125;
 688         ref_ND = 1 - ref_SD - ref_DD;
 689       break;
 690       default:
 691         AliFatal(Form("Unknown origin %d, Energy %f", origin, fEnergy));
 692       }
 693   }
 694   else if(TMath::Abs(fEnergy-1800)<epsilon) {
 695     switch (origin)
 696       {
 697       case 20: // E710, 1.8 TeV
 698       AliInfo("E710 x-sections hadron level definition for Pythia");
 699       // ND: 0.705709
 700       // SD: 0.166590 --> this should be 15.9
 701       // DD: 0.127701
 702       ref_SD = 0.217 * 0.159 / 0.167; error_SD = 0.024 * 0.217 / 0.167;
 703       ref_DD = 0.075 * 1.43;          error_DD = 0.02 * 1.43;
 704       ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
 705       break;
 706
 707       case 21: // E710, 1.8 TeV
 708         AliInfo("E710 x-sections hadron level definition for Phojet");
 709         // ND: 0.817462
 710         // SD: 0.125506 --> this should be 15.9
 711         // DD: 0.057032
 712         ref_SD = 0.161 * 0.159 / 0.126; error_SD = 0.024 * 0.161 / 0.126;
 713         ref_DD = 0.075 * 1.43;         error_DD = 0.02 * 1.43;
 714         ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
 715         break;
 716
 717       case 1: // data 1.8 TeV
 718         AliInfo("??? x-sections");
 719         ref_SD = 0.152;
 720         ref_DD = 0.092;
 721         ref_ND = 1 - ref_SD - ref_DD;
 722         break;
 723       default:
 724         AliFatal(Form("Unknown origin %d, Energy %f", origin, fEnergy));
 725       }
 726   }
 727   else {
 728     AliFatal(Form("Unknown energy %f", fEnergy));
 729   }
 730
 731 }
 732
 733