PWG2/FORWARD/analysis2/AliForwardMCCorrectionsTask.cxx

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