PWG2/FORWARD/analysis2/AliForwardMCCorrections.cxx

   1 #include "AliForwardMCCorrections.h"
   2 #include "AliTriggerAnalysis.h"
   3 #include "AliPhysicsSelection.h"
   4 #include "AliLog.h"
   5 #include "AliFMDAnaParameters.h"
   6 #include "AliESDEvent.h"
   7 #include "AliAODHandler.h"
   8 #include "AliMultiplicity.h"
   9 #include "AliInputEventHandler.h"
  10 #include <TH1.h>
  11 #include <TDirectory.h>
  12 #include <TTree.h>
  13
  14 //====================================================================
  15 namespace {
  16   const char* GetEventName(Bool_t tr, Bool_t vtx)
  17   {
  18     return Form("nEvents%s%s", (tr ? "Tr" : ""), (vtx ? "Vtx", ""));
  19   }
  20   const char* GetHitsName(UShort_t d, Char_t r)
  21   {
  22     return Form("hitsFMD%d%c", d, r);
  23   }
  24   const char* GetStripsName(UShort_t d, Char_t r)
  25   {
  26     return Form("stripsFMD%d%c", d, r);
  27   }
  28   const char* GetPrimaryName(Char_t r, Bool_t trVtx)
  29   {
  30     return Form("primaries%s%s",
  31                 ((r == 'I' || r == 'i') ? "Inner" : "Outer"),
  32                 (trVtx ? "TrVtx" : "All"));
  33   }
  34 }
  35
  36 //====================================================================
  37 AliForwardMCCorrections::AliForwardMCCorrections()
  38   : AliAnalysisTaskSE(),
  39     fHEventsTr(0),
  40     fHEventsTrVtx(0),
  41     fHTriggers(0),
  42     fVtxAxis(),
  43     fEtaAxis()
  44 {
  45 }
  46
  47 //____________________________________________________________________
  48 AliForwardMCCorrections::AliForwardMCCorrections(const char* name)
  49   : AliAnalysisTaskSE(name),
  50     fHEventsTr(0),
  51     fHEventsTrVtx(0),
  52     fHTriggers(0),
  53     fVtxAxis(10,-10,10),
  54     fEtaAxis(200,-4,6)
  55 {
  56   DefineOutput(1, TList::Class());
  57 }
  58
  59 //____________________________________________________________________
  60 AliForwardMCCorrections::AliForwardMCCorrections(const AliForwardMCCorrections& o)
  61   : AliAnalysisTaskSE(o),
  62     fHEventsTr(o.fHEventsTr),
  63     fHEventsTrVtx(o.fHEventsTrVtx),
  64     fHTriggers(o.fHTriggers),
  65     fVtxAxis(o.fVtxAxis),
  66     fEtaAxis(o.fEtaAxis)
  67 {
  68 }
  69
  70 //____________________________________________________________________
  71 AliForwardMCCorrections&
  72 AliForwardMCCorrections::operator=(const AliForwardMCCorrections& o)
  73 {
  74   fHEventsTr         = o.fHEventsTr;
  75   fHEventsTrVtx      = o.fHEventsTrVtx;
  76   fHTriggers         = o.fHTriggers;
  77   fHData             = o.fHData;
  78   fVtxAxis           = o.fVtxAxis;
  79   fEtaAxis           = o.fEtaAxis;
  80
  81   return *this;
  82 }
  83
  84 //____________________________________________________________________
  85 void
  86 AliForwardMCCorrections::Init()
  87 {
  88 }
  89
  90 //____________________________________________________________________
  91 void
  92 AliForwardMCCorrections::SetVertexAxis(Int_t nBin, Double_t min, Double_t max)
  93 {
  94   if (max < min) max = -min;
  95   if (min < max) {
  96     Double_t tmp = min;
  97     min          = max;
  98     max          = tmp;
  99   }
 100   if (min < -10)
 101     AliWarning(Form("Minimum vertex %f < -10, make sure you want this",min));
 102   if (max > +10)
 103     AliWarning(Form("Minimum vertex %f > +10, make sure you want this",max));
 104   fVtxAxis.Set(nBin, min, max);
 105 }
 106 //____________________________________________________________________
 107 void
 108 AliForwardMCCorrections::SetVertexAxis(const TAxis& axis)
 109 {
 110   SetVertexAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
 111 }
 112
 113 //____________________________________________________________________
 114 void
 115 AliForwardMCCorrections::SetEtaAxis(Int_t nBin, Double_t min, Double_t max)
 116 {
 117   if (max < min) max = -min;
 118   if (min < max) {
 119     Double_t tmp = min;
 120     min          = max;
 121     max          = tmp;
 122   }
 123   if (min < -4)
 124     AliWarning(Form("Minimum eta %f < -4, make sure you want this",min));
 125   if (max > +6)
 126     AliWarning(Form("Minimum eta %f > +6, make sure you want this",max));
 127   fVtxAxis.Set(nBin, min, max);
 128 }
 129 //____________________________________________________________________
 130 void
 131 AliForwardMCCorrections::SetVertexAxis(const TAxis& axis)
 132 {
 133   SetEtaAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
 134 }
 135
 136 //____________________________________________________________________
 137 void
 138 AliForwardMCCorrections::InitializeSubs()
 139 {
 140 }
 141
 142 //____________________________________________________________________
 143 TH3D*
 144 AliForwardMCCorrections::Make3D(const char* name, const char* title,
 145                                Int_t nPhi) const
 146 {
 147   TH3D* ret = new TH3D(name, title,
 148                        fVtxAxis.GetNbins(),
 149                        fVtxAxis.GetXmin(),
 150                        fVtxAxis.GetXmax(),
 151                        fEtaAxis.GetNbins(),
 152                        fEtaAxis.GetXmin(),
 153                        fEtaAxis.GetXmax(),
 154                        nPhi, 0, 2*TMath::Pi());
 155   ret->SetXTitle("v_{z} [cm]");
 156   ret->SetYTitle("#eta");
 157   ret->SetZTitle("#varphi [radians]");
 158   ret->SetDirectory(0);
 159   ret->SetStats(0);
 160   ret->Sumw2();
 161
 162   return ret;
 163 }
 164 //____________________________________________________________________
 165 TH1D*
 166 AliForwardMCCorrections::Make1D(const char* name, const char* title) const
 167 {
 168   TH1D* ret = new TH1D(name, title,
 169                        fEtaAxis.GetNbins(),
 170                        fEtaAxis.GetXmin(),
 171                        fEtaAxis.GetXmax());
 172   ret->SetXTitle("#eta");
 173   ret->SetFillColor(kRed+1);
 174   ret->SetFillStyle(3001);
 175   ret->SetDirectory(0);
 176   ret->SetStats(0);
 177   ret->Sumw2();
 178
 179   return ret;
 180 }
 181
 182 //____________________________________________________________________
 183 void
 184 AliForwardMCCorrections::UserCreateOutputObjects()
 185 {
 186   fList = new TList;
 187   fList->SetName(GetName());
 188
 189   fHEvents = new TH1I(GetEventsName(false,false),
 190                       "Number of all events",
 191                       fAxis.GetNbins(),
 192                       fAxis.GetXmin(),
 193                       fAxis.GetXmax());
 194   fHEvents->SetXTitle("v_{z} [cm]");
 195   fHEvents->SetYTitle("# of events");
 196   fHEvents->SetFillColor(kBlue+1);
 197   fHEvents->SetFillStyle(3001);
 198   fHEvents->SetDirectory(0);
 199   fList->Add(fHEvents);
 200
 201   fHEventsTr = new TH1I(GetEventsName(true, false),
 202                         fVtxAxis.GetNbins(),
 203                         fVtxAxis.GetXmin(),
 204                         fVtxAxis.GetXmax());
 205   fHEventsTr->SetXTitle("v_{z} [cm]");
 206   fHEventsTr->SetYTitle("# of events");
 207   fHEventsTr->SetFillColor(kRed+1);
 208   fHEventsTr->SetFillStyle(3001);
 209   fHEventsTr->SetDirectory(0);
 210   fList->Add(fHEventsTr);
 211
 212   fHEventsTrVtx = new TH1I(GetEventsName(true,true),
 213                            "Number of events w/trigger and vertex",
 214                            fVtxAxis.GetNbins(),
 215                            fVtxAxis.GetXmin(),
 216                            fVtxAxis.GetXmax());
 217   fHEventsTrVtx->SetXTitle("v_{z} [cm]");
 218   fHEventsTrVtx->SetYTitle("# of events");
 219   fHEventsTrVtx->SetFillColor(kBlue+1);
 220   fHEventsTrVtx->SetFillStyle(3001);
 221   fHEventsTrVtx->SetDirectory(0);
 222   fList->Add(fHEventsTrVtx);
 223
 224   fHEventsVtx = new TH1I(GetEventsName(false,true),
 225                          "Number of events w/vertex",
 226                          fVtxAxis.GetNbins(),
 227                          fVtxAxis.GetXmin(),
 228                          fVtxAxis.GetXmax());
 229   fHEventsVtx->SetXTitle("v_{z} [cm]");
 230   fHEventsVtx->SetYTitle("# of events");
 231   fHEventsVtx->SetFillColor(kBlue+1);
 232   fHEventsVtx->SetFillStyle(3001);
 233   fHEventsVtx->SetDirectory(0);
 234   fList->Add(fHEventsVtx);
 235
 236
 237   fHTriggers = new TH1I("triggers", "Triggers", 10, 0, 10);
 238   fHTriggers->SetFillColor(kRed+1);
 239   fHTriggers->SetFillStyle(3001);
 240   fHTriggers->SetStats(0);
 241   fHTriggers->SetDirectory(0);
 242   fHTriggers->GetXaxis()->SetBinLabel(1,"INEL");
 243   fHTriggers->GetXaxis()->SetBinLabel(2,"INEL>0");
 244   fHTriggers->GetXaxis()->SetBinLabel(3,"NSD");
 245   fHTriggers->GetXaxis()->SetBinLabel(4,"Empty");
 246   fHTriggers->GetXaxis()->SetBinLabel(5,"A");
 247   fHTriggers->GetXaxis()->SetBinLabel(6,"B");
 248   fHTriggers->GetXaxis()->SetBinLabel(7,"C");
 249   fHTriggers->GetXaxis()->SetBinLabel(8,"E");
 250   fList->Add(fHTriggers);
 251
 252   fPrimaryInnerAll   = Make3D(GetPrimaryName('I',false), "Primary particles, "
 253                               "20 #varphi bins, all events", 20);
 254   fPrimaryOuterAll   = Make3D(GetPrimaryName('O',false), "Primary particles, "
 255                               "40 #varphi bins, all events", 40);
 256   fPrimaryInnerTrVtx = Make3D(GetPrimaryName('I',true), "Primary particles, "
 257                               "20 #varphi bins, Tr+Vtx events", 20);
 258   fPrimaryOuterTrVtx = Make3D(GetPrimaryName('O',true), "Primary particles, "
 259                               "40 #varphi bins, Tr+Vtx events", 40);
 260   TList* primaries = new TList;
 261   primaries->SetName("primaries");
 262   primaries->Add(fPrimaryInnerAll);
 263   primaries->Add(fPrimaryOuterAll);
 264   primaries->Add(fPrimaryInnerTrVtx);
 265   primaries->Add(fPrimaryOuterTrVtx);
 266   fList->Add(primaries);
 267
 268
 269   fHitsFMD1i   = Make3D(GetHitsName(1,'i'),   "All hits in FMD1i, Tr+Vtx", 20);
 270   fHitsFMD2i   = Make3D(GetHitsName(2,'i'),   "All hits in FMD2i, Tr+Vtx", 20);
 271   fHitsFMD2o   = Make3D(GetHitsName(2,'o'),   "All hits in FMD2o, Tr+Vtx", 40);
 272   fHitsFMD3i   = Make3D(GetHitsName(3,'i'),   "All hits in FMD3i, Tr+Vtx", 20);
 273   fHitsFMD3o   = Make3D(GetHitsName(3,'o'),   "All hits in FMD3o, Tr+Vtx", 40);
 274   TList* hits = new TList;
 275   hits->SetName("hits");
 276   hits->Add(fHitsFMD1i);
 277   hits->Add(fHitsFMD2i);
 278   hits->Add(fHitsFMD2o);
 279   hits->Add(fHitsFMD3i);
 280   hits->Add(fHitsFMD3o);
 281   fList->Add(hits);
 282
 283   fStripsFMD1i = Make1D(GetStripsName(1,'i'), "# hit strips in FMD1i, Tr+Vtx");
 284   fStripsFMD2i = Make1D(GetStripsName(2,'i'), "# hit strips in FMD2i, Tr+Vtx");
 285   fStripsFMD2o = Make1D(GetStripsName(2,'o'), "# hit strips in FMD2o, Tr+Vtx");
 286   fStripsFMD3i = Make1D(GetStripsName(3,'i'), "# hit strips in FMD3i, Tr+Vtx");
 287   fStripsFMD3o = Make1D(GetStripsName(3,'o'), "# hit strips in FMD3o, Tr+Vtx");
 288   TList* strips = new TList;
 289   strips->SetName("strips");
 290   strips->Add(fStripsFMD1i);
 291   strips->Add(fStripsFMD2i);
 292   strips->Add(fStripsFMD2o);
 293   strips->Add(fStripsFMD3i);
 294   strips->Add(fStripsFMD3o);
 295   fList->Add(strips);
 296
 297   PostData(1, fList);
 298 }
 299 //____________________________________________________________________
 300 void
 301 AliForwardMCCorrections::UserExec(Option_t*)
 302 {
 303   // Get the input data - MC event
 304   AliMCEvent*  mcEvent = MCEvent();
 305   if (mcevent) {
 306     AliWarning("No MC event found");
 307     return;
 308   }
 309
 310   // Get the input data - ESD event
 311   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 312   if (!esd) {
 313     AliWarning("No ESD event found for input event");
 314     return;
 315   }
 316
 317   // Get the particle stack
 318   AliStack* stack = mcEvent->Stack();
 319
 320   // Get the event generate header
 321   AliHeader*          mcHeader  = mcEvent->Header();
 322   AliGenEventHeader*  genHeader = mcHeader->GetCenEventHeader();
 323
 324   // Get the generator vertex
 325   TArrayF mcVertex;
 326   genHeader->PrimaryVertex();
 327   Double_t mcVtxZ = mcVertex.At(2);
 328
 329   // Check the MC vertex is in range
 330   Int_t mcVtxBin = fVtxAxis.FindBin(vZ);
 331   if (vZ < 1 || vZ > fVtxAxis.GetNbins()) {
 332 #ifdef VERBOSE
 333     AliWarning(Form("MC v_z=%f is out of rante [%,%f]",
 334                     vZ, fVtxAxis.GetXmin(), fVtxAxis.GetXmax()));
 335 #endif
 336     return;
 337   }
 338
 339   // Get the triggers
 340   UInt_t triggers = 0;
 341   Bool_t gotTriggers = AliForwardUtil::ReadTriggers(esd,triggers,fHTriggers);
 342   Bool_t gotInel     = triggers & AliAODForwardMult::kInel;
 343
 344   // Get the ESD vertex
 345   Double_t vZ = -1000000;
 346   Bool_t gotVertex = AliForwardUtil::ReadVertex(esd,vZ);
 347
 348
 349   // Fill the event count histograms
 350   if (gotInel)              fHEventsTr->Fill(mcVtxZ);
 351   if (gotInel && gotVertex) fHEventsTrVtx->Fill(mcVtxZ);
 352   if (gotVertex)            fHEventsVtx->Fill(mcVtxZ);
 353   fHEvents->Fill(mcVtxZ);
 354
 355   // Cache of the hits
 356   AliFMDFloatMap hitsByStrip;
 357   AliFMDFloatMap lastByStrip;
 358
 359   // Loop over all tracks
 360   Int_t nTracks = mcEvent->GetNumberOfTracks();
 361   for (Int_t iTr = 0; iTr < nTracks; iTr++) {
 362     AliMCParticle* particle =
 363       static_cast<AliMCParticle*>(mcEvent->GetTrack(iTr));
 364
 365     // Check the returned particle
 366     if (!particle) continue;
 367
 368     // Check if this charged and a primary
 369     Bool_t isCharged = particle->Charge() != 0;
 370     Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);
 371
 372     // Fill (eta,phi) of the particle into histograsm for b
 373     Double_t eta = particle->Eta();
 374     Double_t phi = partcile->Phi();
 375
 376     if (isCharged && isPrimary)
 377       FillPrimary(gotInel, gotVertex, mcVtxZ, eta, phi);
 378
 379     // For the rest, ignore non-collisions, and events out of vtx range
 380     if (!gotInel || gotVtx) continue;
 381
 382     Int_t nTrRef = particle->GetNumberOfTrackReferences();
 383     for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) {
 384       AliTrackReference* trackRef = particle->GetTrackReference(iTrRef);
 385
 386       // Check existence
 387       if (!ret) continue;
 388
 389       // Check that we hit an FMD element
 390       if (ref->DetectorId() != AliTrackReference::kFMD)
 391         continue;
 392
 393       // Get the detector coordinates
 394       UShort_t d, s, t;
 395       Char_t r;
 396       AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
 397
 398       // Check if mother (?) is charged and that we haven't dealt with it
 399       // already
 400       Int_t lastTrack = Int_t(lastByStrip(d,r,s,t));
 401       if (!isCharged || iTr == lastTrack) continue;
 402
 403       // Increase counter of hits per strip
 404       hitsByStrip(d,r,s,t) += 1;
 405
 406       Double_t trRefEta = esd->GetFMDData()->Eta(d,r,s,t);
 407       Double_t trRefPhi = esd->GetFMDData()->Phi(d,r,s,t);
 408
 409       // Fill strip information into histograms
 410       FillStrip(mcVtxZ, eta, phi, hitsByStrip(d,r,s,t) == 1);
 411
 412       // Set the last processed track number - marking it as done for
 413       // this strip
 414       lastByStrip(d,r,s,t) = Float_t(iTr);
 415
 416       // Flag neighboring strips too
 417       Int_t nStrip = (r == 'I' || r == 'i' ? 512 : 256);
 418       if (t > 0)        lastByStrip(d,r,s,t-1) = Float_t(iTr);
 419       if (t < nStrip-1) lastByStrip(d,r,s,t+1) = Float_t(iTr);
 420     }
 421   }
 422 }
 423
 424 //____________________________________________________________________
 425 void
 426 AliForwardMCCorrections::FillPrimary(Bool_t gotInel, Bool_t gotVtx,
 427                                     Double_t vz, Double_t eta, Double_t phi)
 428 {
 429   if (gotInel && gotVtx) {
 430     // This takes the place of hPrimary_FMD_<r>_vtx<v> and
 431     // Analysed_FMD<r>_vtx<v>
 432     fPrimaryInnerTrVtx->Fill(vz,eta,phi);
 433     fPrimaryOuterTrVtx->Fill(vz,eta,phi);
 434   }
 435   // This takes the place of Inel_FMD<r>_vtx<v> and
 436   // Analysed_FMD<r>_vtx<v>
 437   fPrimaryInnerAll->Fill(vz,eta,phi);
 438   fPrimaryOuterAll->Fill(vz,eta,phi);
 439 }
 440
 441 //____________________________________________________________________
 442 void
 443 AliForwardMCCorrections::FillStrip(UShort_t d, Char_t r,
 444                                   Double_t vz, Double_t eta, Double_t phi,
 445                                   Bool_t first)
 446 {
 447   // Number of hit strips per eta bin
 448   TH1D* hits   = 0;
 449   // All hits per ring, vertex in eta,phi bins.  This takes the place
 450   // of hHits_FMD<d><r>_vtx<v> and DoubleHits_FMD<d><r> (the later
 451   // being just the 2D projection over the X bins)
 452   TH3D* strips = 0;
 453   switch (d) {
 454   case 1:
 455     hits   = fHitsFMD1i;
 456     strips = fStripsFMD1i;
 457     break;
 458   case 2:
 459     hits   = (r == 'I' || r == 'i' ? fHitsFMD2i   : fHitsFMD2o);
 460     strips = (r == 'I' || r == 'i' ? fStripsFMD2i : fStripsFMD2o);
 461     break;
 462   case 3:
 463     hits   = (r == 'I' || r == 'i' ? fHitsFMD3i   : fHitsFMD3o);
 464     strips = (r == 'I' || r == 'i' ? fStripsFMD3i : fStripsFMD3o);
 465     break;
 466   }
 467   if (!hits || !strips) return;
 468
 469   if (first) strips->Fill(eta);
 470
 471   hits->Fill(vz, eta, phi);
 472 }
 473
 474 //____________________________________________________________________
 475 TH2D*
 476 AliForwardMCCorrections::GetVertexProj(Int_t v, TH3D* src) const
 477 {
 478   Int_t nX = src->GetNbinsX();
 479   if (v > nX || v < 1) return 0;
 480
 481   src->GetXaxis()->SetRange(v,v+1);
 482
 483   TH2D* ret   = static_cast<TH2D*>(src->Project3D("zy"));
 484   ret->SetName(Form("%s_vtx%02d", src->GetName(), v));
 485   ret->SetDirectory(0);
 486
 487   src->GetXaxis()->SetRange(0,nX+1);
 488
 489   return ret;
 490 }
 491
 492
 493 //____________________________________________________________________
 494 void
 495 AliForwardMCCorrections::Terminate(Option_t*)
 496 {
 497   TList* list = dynamic_cast<TList*>(GetOutputData(1));
 498   if (!list) {
 499     AliError("No output list defined");
 500     return;
 501   }
 502
 503   TList* primaries = static_cast<TList*>(list->FindObject("primaries"));
 504   TList* hits      = static_cast<TList*>(list->FindObject("hits"));
 505   TList* strips    = static_cast<TList*>(list->FindObject("strips"));
 506   if (!primaries || !hits || !strips) {
 507     AliError(Form("Could not find all sub-lists in %s (p=%p,h=%p,s=%p)",
 508                   list->GetName(), primaries, hits, strips));
 509     return;
 510   }
 511
 512   TH1I* eventsAll =
 513     static_cast<TH1I*>(list->FindObject(GetEventsName(false,false)));
 514   TH1I* eventsTr =
 515     static_cast<TH1I*>(list->FindObject(GetEventsName(true,false)));
 516   TH1I* eventsVtx =
 517     static_cast<TH1I*>(list->FindObject(GetEventsName(false,true)));
 518   TH1I* eventsTrVtx =
 519     static_cast<TH1I*>(list->FindObject(GetEventsName(true,true)));
 520   if (!eventsAll || !eventsTr || !eventsVtx || !eventsTrVtx) {
 521     AliError(Form("Could not find all event histograms in %s "
 522                   "(a=%p,t=%p,v=%p,tv=%p)", list->GetName(),
 523                   eventsAll, eventsTr, eventsVtx, eventsTrVtx));
 524     return;
 525   }
 526
 527   TH3D* primIAll =
 528     static_cast<TH3D*>(primaries->FindObject(GetPrimaryName('I',false)));
 529   TH3D* primOAll =
 530     static_cast<TH3D*>(primaries->FindObject(GetPrimaryName('O',false)));
 531   TH3D* primITrVtx =
 532     static_cast<TH3D*>(primaries->FindObject(GetPrimaryName('I',true)));
 533   TH3D* primOTrVtx =
 534     static_cast<TH3D*>(primaries->FindObject(GetPrimaryName('O',true)));
 535   if (!primIAll || !primOAll || !primITrVtx || !primOTrVtx) {
 536     AliError(Form("Could not find all primary particle histograms in %s "
 537                   "(ai=%p,ao=%p,tvi=%p,tvo=%p)", list->GetName(),
 538                   primIAll, primOAll, primITrVtx, primOTrVtx));
 539     return;
 540   }
 541
 542   TH3D* hits1i = static_cast<TH3D*>(hits->FindObject(GetHitsName(1,'i')));
 543   TH3D* hits2i = static_cast<TH3D*>(hits->FindObject(GetHitsName(2,'i')));
 544   TH3D* hits2o = static_cast<TH3D*>(hits->FindObject(GetHitsName(2,'o')));
 545   TH3D* hits3i = static_cast<TH3D*>(hits->FindObject(GetHitsName(3,'i')));
 546   TH3D* hits3o = static_cast<TH3D*>(hits->FindObject(GetHitsName(3,'o')));
 547   if (!hits1i || !hits2i || !hits2o || hits3i || hits3o) {
 548     AliError(Form("Could not find all ring hits histograms in %s "
 549                   "(1i=%p,2i=%p,2o=%p,3i=%p,3o=%p)", hits->GetName(),
 550                   hits1i, hits2i, hits2o, hits3i, hits3o));
 551     return;
 552   }
 553
 554   TH1D* strips1i = static_cast<TH1D*>(strips->FindObject(GetStripsName(1,'i')));
 555   TH1D* strips2i = static_cast<TH1D*>(strips->FindObject(GetStripsName(2,'i')));
 556   TH1D* strips2o = static_cast<TH1D*>(strips->FindObject(GetStripsName(2,'o')));
 557   TH1D* strips3i = static_cast<TH1D*>(strips->FindObject(GetStripsName(3,'i')));
 558   TH1D* strips3o = static_cast<TH1D*>(strips->FindObject(GetStripsName(3,'o')));
 559   if (!strips1i || !strips2i || !strips2o || strips3i || strips3o) {
 560     AliError(Form("Could not find all ring strips histograms in %s "
 561                   "(1i=%p,2i=%p,2o=%p,3i=%p,3o=%p)", strips->GetName(),
 562                   strips1i, strips2i, strips2o, strips3i, strips3o));
 563     return;
 564   }
 565
 566   // Calculate the over-all event selection efficiency
 567   TH1D* selEff = new TH1D("selEff", "Event selection efficiency",
 568                           fVtxArray.GetNbins(),
 569                           fVtxArray.GetXmin(),
 570                           fVtxArray.GetXmax());
 571   selEff->Sumw2();
 572   selEff->SetDirectory(0);
 573   selEff->SetFillColor(kMagenta+1);
 574   selEff->SetFillStyle(3001);
 575   selEff->Add(eventsAll);
 576   selEff->Divide(eventsTrVtx);
 577   list->Add(selEff);
 578
 579   // Loop over vertex bins and do vertex dependendt stuff
 580   for (Int_t v = 1; v <= fVtxAxis.GetNbins(); v++) {
 581     // Make a sub-list in the output
 582     TList* vl = new TList;
 583     vl->SetName("vtx%02d", v);
 584     list->Add(vl);
 585
 586     // Get event counts
 587     Int_t nEventsAll   = eventsAll->GetEntries(v);
 588     Int_t nEventsTr    = eventsTr->GetEntries(v);
 589     Int_t nEventsVtx   = eventsVtx->GetEntries(v);
 590     Int_t nEventsTrVtx = eventsTrVtx->GetEntries(v);
 591
 592     // Project event histograms, set names, and store
 593     TH2D* primIAllV   = GetVertexProj(v, primIAll);
 594     TH2D* primOAllV   = GetVertexProj(v, primOAll);
 595     TH2D* primITrVtxV = GetVertexProj(v, primITrVtx);
 596     TH2D* primOTrVtxV = GetVertexProj(v, primOTrVtx);
 597     vl->Add(primIAllV);
 598     vl->Add(primOAllV);
 599     vl->Add(primITrVtxV);
 600     vl->Add(primOTrVtxV);
 601
 602     primIAllV->Scale(1. / nEventAll);
 603     primOAllV->Scale(1. / nEventAll);
 604     primITrVtxV->Scale(1. / nEventTr);
 605     primOTrVtxV->Scale(1. / nEventTr);
 606
 607     // Calculate the vertex bias on the d^2N/detadphi
 608     TH2D* selBiasI = static_cast<TH2D*>(primITrVtxV->Clone("selBiasI%02d",v));
 609     TH2D* selBiasO = static_cast<TH2D*>(primOTrVtxV->Clone("selBiasO%02d",v));
 610     selBiasI->SetTitle(Form("Event selection bias in vertex bin %d", v));
 611     selBiasO->SetTitle(Form("Event selection bias in vertex bin %d", v));
 612     selBiasI->Divide(primIAllV);
 613     selBiasO->Divide(primOAllV);
 614     vl->Add(selBiasI);
 615     vl->Add(selBiasO);
 616
 617     for (UShort_t d = 1; d <= 3; d++) {
 618       UShort_t nQ = (d == 1 ? 1 : 2);
 619       for (UShort_t q = 0; q < nQ; q++) {
 620         Char_t   r  = (q == 0 ? 'I' : 'O');
 621         TH3D* hits3D = 0;
 622         TH2D* prim2D = (q == 0 ? primITrVtxV : primOTrVtxV);
 623         switch (d) {
 624         case 1: hits3D = hits1i; break;
 625         case 2: hits3D = (q == 0 ? hits2i : hits2o); break;
 626         case 3: hits3D = (q == 0 ? hits3i : hits3o); break;
 627         }
 628
 629         TH2D* sec = GetVertexProj(v, hits3D);
 630         sec->SetName(Form("secondaryFMD%d%c_vtx%02d", d, r, v));
 631         sec->SetTitle(Form("Secondary correction map for FMD%d%c "
 632                            "in vertex bin %d", d, r, v));
 633         sec->Divide(prim2D);
 634         vl->Add(sec);
 635
 636         if (v > 1) continue;
 637
 638         // Do the double hit correction (only done once per ring in
 639         // the vertex loop)
 640         TH1D* strips = 0;
 641         switch (d) {
 642         case 1: strips = strips1i; break;
 643         case 2: strips = (q == 0 ? strips2i : strips2o); break;
 644         case 3: strips = (q == 0 ? strips3i : strips3o); break;
 645         }
 646
 647         TH2D* hits2D    = GetVertexProj(v, hits3D);
 648         TH1D* doubleHit = hits2D->ProjectionX(Form("doubleHitFMD%d%c",d,r));
 649         doubleHit->SetTitle(Form("Double hit correction for FMD%d%c",d,r));
 650         doubleHit->SetDirectory(0);
 651         doubleHit->SetFillColor(kGreen+1);
 652         doubleHit->SetFillStyle(3001);
 653         doubleHit->Sumw2();
 654         doubleHit->Divide(strips);
 655         list->Add(doubleHit);
 656       }
 657     }
 658   }
 659 }
 660
 661 //____________________________________________________________________
 662 void
 663 AliForwardMCCorrections::Print(Option_t*) const
 664 {
 665 }
 666
 667 //
 668 // EOF
 669 //