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