PWG2/FORWARD/analysis2/AliForwardMCCorrectionsTask.cxx

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