ANALYSIS/AliCollisionNormalization.cxx

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