PWG2/FORWARD/analysis2/MakeEvaluateTriggers.C

   1 #ifdef BUILD
   2 #include "AliAnalysisTaskSE.h"
   3 #include "AliAnalysisManager.h"
   4 #include "AliAnalysisDataContainer.h"
   5 #include "AliMCEvent.h"
   6 #include "AliESDEvent.h"
   7 #include "AliStack.h"
   8 #include "AliMultiplicity.h"
   9 #include "AliFMDMCEventInspector.h"
  10 #include "AliAODForwardMult.h"
  11 #include "AliLog.h"
  12 #include <TH1I.h>
  13 #include <TH2D.h>
  14 #include <TAxis.h>
  15 #include <TList.h>
  16 #include <TObjArray.h>
  17 #include <TParameter.h>
  18 #include <TStopwatch.h>
  19 #include <TROOT.h>
  20 #include <THStack.h>
  21 #include <TStyle.h>
  22
  23 //====================================================================
  24 /**
  25  * Task to evaluate trigger bias in pp
  26  *
  27  */
  28 class EvaluateTrigger : public AliAnalysisTaskSE
  29 {
  30 public:
  31   enum {
  32     kNone = 0x0,
  33     kESD  = 0x1,
  34     kMC   = 0x2
  35   };
  36   /**
  37    * Constructor
  38    */
  39   EvaluateTrigger()
  40     : AliAnalysisTaskSE(),
  41       fInel(),
  42       fInelGt0(),
  43       fNSD(),
  44       fInspector(),
  45       fFirstEvent(true),
  46       fData(0),
  47       fTriggers(0),
  48       fTrackletRequirement(kESD),
  49       fVertexRequirement(kESD),
  50       fVertexAxis(0, 0, 0),
  51       fVertexESD(0),
  52       fVertexMC(0),
  53       fM(0)
  54   {}
  55   /**
  56    * Constructor
  57    */
  58   EvaluateTrigger(const char* /*name*/)
  59     : AliAnalysisTaskSE("evaluateTrigger"),
  60       fInel(AliAODForwardMult::kInel),
  61       fInelGt0(AliAODForwardMult::kInelGt0),
  62       fNSD(AliAODForwardMult::kNSD),
  63       fInspector("eventInspector"),
  64       fFirstEvent(true),
  65       fData(0),
  66       fTriggers(0),
  67       fTrackletRequirement(kESD),
  68       fVertexRequirement(kESD),
  69       fVertexAxis(10, -10, 10),
  70       fVertexESD(0),
  71       fVertexMC(0),
  72       fM(0)
  73   {
  74     DefineOutput(1, TList::Class());
  75     DefineOutput(2, TList::Class());
  76   }
  77   void SetVertexRequirement(UShort_t m) { fVertexRequirement = m; }
  78   void SetTrackletRequirement(UShort_t m) { fTrackletRequirement = m; }
  79   void SetVertexAxis(Int_t nBins, Double_t low, Double_t high)
  80   {
  81     fVertexAxis.Set(nBins, low, high);
  82   }
  83   /**
  84    * Intialize
  85    */
  86   void Init() {}
  87   /**
  88    * Create output objects
  89    */
  90   void UserCreateOutputObjects()
  91   {
  92     fList = new TList;
  93     fList->SetOwner();
  94     fList->SetName(GetName());
  95
  96     Double_t mb[] = { 0, 1, 2, 3, 4, 5, 6, 8, 9, 10, 1000 };
  97     Int_t    nM   = 10;
  98     TAxis mAxis(nM, mb);
  99     TAxis eAxis(200, -4, 6);
 100     TAxis pAxis(40, 0, 2*TMath::Pi());
 101
 102     fData = new TH2D("data", "Cache",
 103                      eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax(),
 104                      pAxis.GetNbins(), pAxis.GetXmin(), pAxis.GetXmax());
 105     fData->SetDirectory(0);
 106     fData->SetXTitle("#eta");
 107     fData->SetYTitle("#varphi [radians]");
 108     fData->SetZTitle("N_{ch}(#eta,#varphi)");
 109     fData->Sumw2();
 110
 111     fM = new TH1D("m", "Distribution of N_{ch}|_{|#eta|<1}",nM+1, -0.5, nM+.5);
 112     fM->SetXTitle("N_{ch}|_{|#eta|<1}");
 113     fM->SetYTitle("Events");
 114     fM->SetFillColor(kRed+1);
 115     fM->SetFillStyle(3001);
 116     fM->SetDirectory(0);
 117     fList->Add(fM);
 118
 119     for (Int_t i = 0; i <= nM; i++) {
 120       TString lbl;
 121       if (i == 0)       lbl = "all";
 122       else if (i == nM) lbl = Form("%d+",int(mAxis.GetBinLowEdge(i)+.5));
 123       else              lbl = Form("<%d",int(mAxis.GetBinUpEdge(i)+.5));
 124       fM->GetXaxis()->SetBinLabel(i+1, lbl);
 125     }
 126
 127     fTriggers = new TH1I("triggers", "Triggers", 6, -.5, 5.5);
 128     fTriggers->SetDirectory(0);
 129     fTriggers->GetXaxis()->SetBinLabel(1, "INEL (MC)");
 130     fTriggers->GetXaxis()->SetBinLabel(2, "INEL (ESD)");
 131     fTriggers->GetXaxis()->SetBinLabel(3, "INEL>0 (MC)");
 132     fTriggers->GetXaxis()->SetBinLabel(4, "INEL>0 (ESD)");
 133     fTriggers->GetXaxis()->SetBinLabel(5, "NSD (MC)");
 134     fTriggers->GetXaxis()->SetBinLabel(6, "NSD (ESD)");
 135     fTriggers->SetFillColor(kYellow+1);
 136     fTriggers->SetFillStyle(3001);
 137     fList->Add(fTriggers);
 138
 139     fVertexESD = new TH1D("vertexESD", "ESD vertex distribution",
 140                           fVertexAxis.GetNbins(),
 141                           fVertexAxis.GetXmin(),
 142                           fVertexAxis.GetXmax());
 143     fVertexESD->SetDirectory(0);
 144     fVertexESD->SetFillColor(kRed+1);
 145     fVertexESD->SetFillStyle(3001);
 146     fVertexESD->SetXTitle("v_{z} [cm]");
 147     fVertexESD->SetYTitle("P(v_{z})");
 148     fList->Add(fVertexESD);
 149
 150     fVertexMC = new TH1D("vertexMC", "MC vertex distribution",
 151                           fVertexAxis.GetNbins(),
 152                           fVertexAxis.GetXmin(),
 153                           fVertexAxis.GetXmax());
 154     fVertexMC->SetDirectory(0);
 155     fVertexMC->SetFillColor(kBlue+1);
 156     fVertexMC->SetFillStyle(3001);
 157     fVertexMC->SetXTitle("v_{z} [cm]");
 158     fVertexMC->SetYTitle("P(v_{z})");
 159     fList->Add(fVertexMC);
 160
 161     fInel.CreateObjects(fList, fM, fData);
 162     fInelGt0.CreateObjects(fList, fM, fData);
 163     fNSD.CreateObjects(fList, fM, fData);
 164
 165
 166     fInspector.DefineOutput(fList);
 167     fInspector.Init(fVertexAxis);
 168
 169     PostData(1, fList);
 170   }
 171   /**
 172    * Event processing
 173    */
 174   void UserExec(Option_t*)
 175   {
 176     // Get the input data - MC event
 177     AliMCEvent*  mcEvent = MCEvent();
 178     if (!mcEvent) {
 179       AliWarning("No MC event found");
 180       return;
 181     }
 182
 183     // Get the input data - ESD event
 184     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 185     if (!esd) {
 186       AliWarning("No ESD event found for input event");
 187       return;
 188     }
 189
 190     if (fFirstEvent && esd->GetESDRun()) {
 191       fInspector.ReadRunDetails(esd);
 192
 193       AliInfo(Form("Initializing with parameters from the ESD:\n"
 194                    "         AliESDEvent::GetBeamEnergy()   ->%f\n"
 195                    "         AliESDEvent::GetBeamType()     ->%s\n"
 196                    "         AliESDEvent::GetCurrentL3()    ->%f\n"
 197                    "         AliESDEvent::GetMagneticField()->%f\n"
 198                    "         AliESDEvent::GetRunNumber()    ->%d\n",
 199                    esd->GetBeamEnergy(),
 200                    esd->GetBeamType(),
 201                    esd->GetCurrentL3(),
 202                    esd->GetMagneticField(),
 203                    esd->GetRunNumber()));
 204
 205       fFirstEvent = false;
 206     }
 207
 208     // Get the particle stack
 209     AliStack* stack = mcEvent->Stack();
 210
 211     // Some variables
 212     UInt_t   triggers; // Trigger bits
 213     Bool_t   lowFlux;  // Low flux flag
 214     UShort_t iVz;      // Vertex bin from ESD
 215     Double_t vZ;       // Z coordinate from ESD
 216     Double_t cent;     // Centrality
 217     UShort_t iVzMc;    // Vertex bin from MC
 218     Double_t vZMc;     // Z coordinate of IP vertex from MC
 219     Double_t b;        // Impact parameter
 220     Int_t    nPart;    // Number of participants
 221     Int_t    nBin;     // Number of binary collisions
 222     Double_t phiR;     // Reaction plane from MC
 223
 224     // Process the data
 225     Int_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ, cent);
 226     Int_t retMC  = fInspector.ProcessMC(mcEvent, triggers, iVzMc,
 227                                         vZMc, b, nPart, nBin, phiR);
 228
 229     Bool_t hasESDVtx = retESD == AliFMDEventInspector::kOk;
 230     Bool_t hasMCVtx  = retMC  == AliFMDEventInspector::kOk;
 231     if (hasESDVtx) fVertexESD->Fill(vZ);
 232     if (hasMCVtx)  fVertexMC->Fill(vZMc);
 233
 234     Bool_t isMcInel = true; // (triggers & AliAODForwardMult::kB);
 235     Bool_t isMcNSD  = (triggers & AliAODForwardMult::kMCNSD);
 236
 237     Int_t mESD = 0;
 238     const AliMultiplicity* spdmult = esd->GetMultiplicity();
 239     if (!spdmult) {
 240       AliWarning("No SPD multiplicity");
 241     }
 242     else {
 243       // Check if we have one or more tracklets
 244       // in the range -1 < eta < 1 to set the INEL>0
 245       // trigger flag.
 246       Int_t n = spdmult->GetNumberOfTracklets();
 247       for (Int_t j = 0; j < n; j++)
 248         if(TMath::Abs(spdmult->GetEta(j)) < 1) mESD++;
 249     }
 250
 251     // Reset cache
 252     fData->Reset();
 253     Int_t mMC = 0; // Number of particles in |eta|<1
 254
 255     // Loop over all tracks
 256     Int_t nTracks = mcEvent->GetNumberOfTracks();
 257     for (Int_t iTr = 0; iTr < nTracks; iTr++) {
 258       AliMCParticle* particle =
 259         static_cast<AliMCParticle*>(mcEvent->GetTrack(iTr));
 260
 261       // Check the returned particle
 262       if (!particle) continue;
 263
 264       // Check if this charged and a primary
 265       Bool_t isCharged = particle->Charge() != 0;
 266       Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);
 267
 268       if (!isCharged || !isPrimary) continue;
 269
 270
 271       // Fill (eta,phi) of the particle into histograsm for b
 272       Double_t eta = particle->Eta();
 273       Double_t phi = particle->Phi();
 274
 275       fData->Fill(eta, phi);
 276       if (TMath::Abs(eta) <= 1) mMC++;
 277     }
 278     Int_t m = mESD;
 279     if (fTrackletRequirement == kMC) m = mMC;
 280     fM->Fill(m);
 281
 282     bool isMcInelGt0 = isMcInel && (mMC > 0);
 283
 284     bool hasVertex   = true;
 285     if (fVertexRequirement & kMC)  hasVertex = hasVertex && hasMCVtx;
 286     if (fVertexRequirement & kESD) hasVertex = hasVertex && hasESDVtx;
 287
 288     if (isMcInel) {
 289       fTriggers->Fill(0);
 290       bool triggered = (triggers & AliAODForwardMult::kInel);
 291       if (triggered) fTriggers->Fill(1);
 292       fInel.AddEvent(triggered, hasVertex, m, fData);
 293     }
 294     if (isMcInelGt0) {
 295       fTriggers->Fill(2);
 296       bool triggered = (triggers & AliAODForwardMult::kInelGt0);
 297       if (triggered) fTriggers->Fill(3);
 298       fInelGt0.AddEvent(triggered, hasVertex, m, fData);
 299     }
 300     if (isMcNSD) {
 301       fTriggers->Fill(4);
 302       bool triggered = (triggers & AliAODForwardMult::kNSD);
 303       if (triggered) fTriggers->Fill(5);
 304       fNSD.AddEvent(triggered, hasVertex, m, fData);
 305     }
 306     PostData(1, fList);
 307   }
 308   /**
 309    * End of job processing
 310    */
 311   void Terminate(Option_t*)
 312   {
 313     fList = dynamic_cast<TList*>(GetOutputData(1));
 314     if (!fList) {
 315       AliError(Form("No output list defined (%p)", GetOutputData(1)));
 316       if (GetOutputData(1)) GetOutputData(1)->Print();
 317       return;
 318     }
 319
 320
 321     TList* output = new TList;
 322     output->SetName(GetName());
 323     output->SetOwner();
 324
 325     fVertexMC = static_cast<TH1D*>(fList->FindObject("vertexMC"));
 326     fVertexESD = static_cast<TH1D*>(fList->FindObject("vertexESD"));
 327     fM         = static_cast<TH1D*>(fList->FindObject("m"));
 328     if (fVertexMC) {
 329       TH1D* vtxMC = static_cast<TH1D*>(fVertexMC->Clone("vertexMC"));
 330       vtxMC->SetDirectory(0);
 331       vtxMC->Scale(1. / vtxMC->GetEntries());
 332       output->Add(vtxMC);
 333     }
 334     if (fVertexESD) {
 335       TH1D* vtxESD = static_cast<TH1D*>(fVertexESD->Clone("vertexESD"));
 336       vtxESD->SetDirectory(0);
 337       vtxESD->Scale(1. / vtxESD->GetEntries());
 338       output->Add(vtxESD);
 339     }
 340     if (fM) {
 341       TH1D* m = static_cast<TH1D*>(fM->Clone("m"));
 342       m->SetDirectory(0);
 343       m->SetYTitle("P(N_{ch}|_{|#eta|<1} < X)");
 344       m->Scale(1. / m->GetBinContent(1));
 345       output->Add(m);
 346     }
 347
 348     fInel.Finish(fList, output);
 349     fInelGt0.Finish(fList, output);
 350     fNSD.Finish(fList, output);
 351
 352     PostData(2, output);
 353   }
 354
 355 protected:
 356   //__________________________________________________________________
 357   /**
 358    * Structure to hold per trigger type information
 359    */
 360   struct TriggerType : public TNamed
 361   {
 362     //________________________________________________________________
 363     /**
 364      * Structure to hold per multiplicity bin information
 365      */
 366     struct MBin : public TNamed
 367     {
 368       TH2D*     fTruth;
 369       TH2D*     fTriggered;
 370       TH2D*     fAccepted;
 371       TH1I*     fCounts;
 372       UShort_t  fLow;
 373       UShort_t  fHigh;
 374       /**
 375        * Constructor
 376        */
 377       MBin()
 378         : fTruth(0), fTriggered(0), fAccepted(0),
 379           fCounts(0), fLow(0), fHigh(1000) {}
 380       /**
 381        * Constructor
 382        *
 383        * @param p      Parent list
 384        * @param low    Low cut
 385        * @param high   High cut
 386        * @param eAxis  Eta axis
 387        * @param pAxis  Phi axis
 388        */
 389       MBin(TList* p, UShort_t low, UShort_t high, const TH2D* dHist)
 390         : fTruth(0),
 391           fTriggered(0),
 392           fAccepted(0),
 393           fCounts(0),
 394           fLow(low),
 395           fHigh(high)
 396       {
 397         if (low >= high) SetName("all");
 398         else             SetName(Form("m%03d_%03d", fLow, fHigh));
 399         TList* l = new TList;
 400         l->SetOwner();
 401         l->SetName(GetName());
 402         p->Add(l);
 403
 404         fTruth = static_cast<TH2D*>(dHist->Clone(("truth")));
 405         fTruth->SetTitle("MC truth");
 406         fTruth->SetDirectory(0);
 407         fTruth->SetZTitle("#sum_i^{N_X} N_{ch}(#eta,#varphi)");
 408         fTruth->Sumw2();
 409         fTruth->Reset();
 410         l->Add(fTruth);
 411
 412         fTriggered = static_cast<TH2D*>(fTruth->Clone(("triggered")));
 413         fTriggered->SetTitle("Triggered");
 414         fTriggered->SetDirectory(0);
 415         fTriggered->SetZTitle("#sum_i^{N_T} N_{ch}(#eta,#varphi)");
 416         fTriggered->Sumw2();
 417         fTriggered->Reset();
 418         l->Add(fTriggered);
 419
 420         fAccepted = static_cast<TH2D*>(fTruth->Clone(("accepted")));
 421         fAccepted->SetTitle("Accepted");
 422         fAccepted->SetDirectory(0);
 423         fAccepted->SetZTitle("#sum_i^{N_T} N_{ch}(#eta,#varphi)");
 424         fAccepted->Sumw2();
 425         fAccepted->Reset();
 426         l->Add(fAccepted);
 427
 428         fCounts = new TH1I("counts", "Event counts", 3, -.5, 2.5);
 429         fCounts->SetDirectory(0);
 430         fCounts->GetXaxis()->SetBinLabel(1, "Truth");
 431         fCounts->GetXaxis()->SetBinLabel(2, "Triggered");
 432         fCounts->GetXaxis()->SetBinLabel(3, "Accepted");
 433         fCounts->SetYTitle("# events");
 434         l->Add(fCounts);
 435       }
 436       /**
 437        * Add event observation
 438        *
 439        * @param triggered Whether the event was triggered
 440        * @param event     Data for this event
 441        */
 442       void AddEvent(Bool_t triggered, Bool_t hasVtx, const TH2D* event)
 443       {
 444         fCounts->Fill(0);
 445         fTruth->Add(event);
 446         if (triggered) {
 447           fCounts->Fill(1);
 448           fTriggered->Add(event);
 449           if (hasVtx) {
 450             fCounts->Fill(2);
 451             fAccepted->Add(event);
 452           }
 453         }
 454       }
 455       /**
 456        * End of job processing
 457        *
 458        * @param p      Parent list
 459        * @param o      Output parent list
 460        * @param stack  Stack of histograms
 461        *
 462        * @return Trigger efficiency
 463        */
 464      Double_t Finish(const TList* p, TList* o, THStack* stack)
 465       {
 466         TList* l = dynamic_cast<TList*>(p->FindObject(GetName()));
 467         if (!l) {
 468           Warning("Finish", "Cannot find %s in %s", GetName(), p->GetName());
 469           return 0;
 470         }
 471         fTruth     = static_cast<TH2D*>(l->FindObject("truth"));
 472         fTriggered = static_cast<TH2D*>(l->FindObject("triggered"));
 473         fAccepted  = static_cast<TH2D*>(l->FindObject("accepted"));
 474         fCounts    = static_cast<TH1I*>(l->FindObject("counts"));
 475
 476         Int_t    nTruth     = fCounts->GetBinContent(1);
 477         Int_t    nTriggered = fCounts->GetBinContent(2);
 478         Int_t    nAccepted  = fCounts->GetBinContent(3);
 479         Double_t eff        = 0;
 480         if (nTruth > 0) eff = double(nTriggered) / nTruth;
 481         else if (nTriggered == nTruth) eff = 1;
 482
 483         if (nTruth > 0)     fTruth->Scale(1. / nTruth);
 484         if (nTriggered > 0) fTriggered->Scale(1. / nTriggered);
 485         if (nAccepted > 0)  fAccepted->Scale(1. / nAccepted);
 486
 487         if (fLow >= fHigh)
 488           Info("Finish", "%-6s  [all]  E_X=N_T/N_X=%9d/%-9d=%f "
 489                "E_V=N_A/N_T=%9d/%-9d=%f",
 490                p->GetName(), nTriggered, nTruth, eff, nAccepted, nTriggered,
 491                (nTriggered > 0 ? double(nAccepted) / nTriggered: 0));
 492         else
 493           Info("Finish", "%-6s [%2d-%2d] E_X=N_T/N_X=%9d/%-9d=%f "
 494                "E_V=N_A/N_T=%9d/%-9d=%f",
 495                p->GetName(), fLow, fHigh, nTriggered, nTruth, eff,
 496                nAccepted, nTriggered,
 497                (nTriggered > 0 ? double(nAccepted) / nTriggered: 0));
 498
 499         TList* out = new TList;
 500         out->SetName(GetName());
 501         out->SetOwner();
 502         o->Add(out);
 503
 504         out->Add(fTruth);
 505         out->Add(fTriggered);
 506         out->Add(new TParameter<double>("eff", eff));
 507
 508         TH2D* bias = static_cast<TH2D*>(fAccepted->Clone("bias"));
 509         bias->Divide(fTruth);
 510         bias->SetDirectory(0);
 511         bias->SetZTitle("Trigger bias (accepted/truth)");
 512         out->Add(bias);
 513
 514         TH1D* truth_px = static_cast<TH1D*>(fTruth->ProjectionX("truth_eta"));
 515         truth_px->SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1} < %d",
 516                                fLow, fHigh));
 517         truth_px->Scale(1. / fTruth->GetNbinsY());
 518         truth_px->SetDirectory(0);
 519         truth_px->SetLineColor(kRed+1);
 520         truth_px->SetMarkerColor(kRed+1);
 521         truth_px->SetFillColor(kRed+1);
 522         truth_px->SetFillStyle(0);
 523         out->Add(truth_px);
 524
 525         TH1D* triggered_px =
 526           static_cast<TH1D*>(fTriggered->ProjectionX("triggered_eta"));
 527         triggered_px->SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1} < %d",
 528                                fLow, fHigh));
 529         triggered_px->Scale(1. / fTriggered->GetNbinsY());
 530         triggered_px->SetDirectory(0);
 531         triggered_px->SetLineColor(kGreen+1);
 532         triggered_px->SetMarkerColor(kGreen+1);
 533         triggered_px->SetFillColor(kGreen+1);
 534         triggered_px->SetFillStyle(0);
 535         out->Add(triggered_px);
 536
 537         TH1D* accepted_px =
 538           static_cast<TH1D*>(fAccepted->ProjectionX("accepted_eta"));
 539         accepted_px->SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1} < %d",
 540                                    fLow, fHigh));
 541         accepted_px->Scale(1. / fAccepted->GetNbinsY());
 542         accepted_px->SetLineColor(kBlue+1);
 543         accepted_px->SetMarkerColor(kBlue+1);
 544         accepted_px->SetFillColor(kBlue+1);
 545         accepted_px->SetDirectory(0);
 546         out->Add(accepted_px);
 547
 548         THStack* data = new THStack("data", "Data distributions");
 549         data->Add(truth_px);
 550         data->Add(triggered_px);
 551         data->Add(accepted_px);
 552         out->Add(data);
 553
 554 #if 0
 555         TH1D* bias_px = static_cast<TH1D*>(bias->ProjectionX("bias_eta"));
 556         bias_px->SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1} < %d",
 557                                fLow, fHigh));
 558         bias_px->Scale(1. / bias->GetNbinsY());
 559 #else
 560         TH1D* bias_px = static_cast<TH1D*>(accepted_px->Clone("bias_px"));
 561         bias_px->Divide(truth_px);
 562         bias_px->SetYTitle("Trigger bias (triggered/truth)");
 563 #endif
 564         bias_px->SetDirectory(0);
 565         bias_px->SetMarkerStyle(20);
 566         bias_px->SetFillStyle(0);
 567         bias_px->SetMinimum(0);
 568         out->Add(bias_px);
 569
 570         stack->Add(bias_px);
 571
 572         return eff;
 573       }
 574       ClassDef(MBin,1);
 575     };
 576
 577     /**
 578      * Constructor
 579      */
 580     TriggerType()
 581       : TNamed(),
 582         fMask(0),
 583         fM(0),
 584         fBins(0)
 585     {}
 586     //--- Constructor ------------------------------------------------
 587     /**
 588      * Constructor
 589      *
 590      * @param mask  Trigger mask
 591      */
 592     TriggerType(UShort_t mask)
 593       : TNamed(AliAODForwardMult::GetTriggerString(mask), ""),
 594         fMask(mask),
 595         fM(0),
 596         fBins(0)
 597     {
 598     }
 599     /**
 600      * Create objects
 601      *
 602      * @param list   PArent list
 603      * @param mAxis  Multiplicity axis
 604      * @param eAxis  Eta axis
 605      * @param pAxis  Phi axis
 606      */
 607     void CreateObjects(TList* list,
 608                        const TH1D* mHist,
 609                        const TH2D* dHist)
 610     {
 611       TList* ours  = new TList;
 612       ours->SetName(GetName());
 613       ours->SetOwner();
 614       list->Add(ours);
 615
 616       fM = static_cast<TH1D*>(mHist->Clone("m"));
 617       fM->SetDirectory(0);
 618       ours->Add(fM);
 619
 620       fBins = new TObjArray;
 621       fBins->AddAt(new MBin(ours, 0, 0, dHist), 0);
 622       for (Int_t i = 1; i <= fM->GetNbinsX(); i++) {
 623         Double_t low  = fM->GetXaxis()->GetBinLowEdge(i);
 624         Double_t high = fM->GetXaxis()->GetBinUpEdge(i);
 625
 626         fBins->AddAt(new MBin(ours, low, high, dHist), i);
 627       }
 628     }
 629     /**
 630      * Find bin corresponding to m
 631      *
 632      * @param m Multiplicity
 633      *
 634      * @return Bin.
 635      */
 636     MBin* FindBin(UShort_t m)
 637     {
 638       Int_t b = fM->GetXaxis()->FindBin(m);
 639       if (b <= 0) return 0;
 640       if (b >= fM->GetNbinsX()+1) b = fM->GetNbinsX();
 641       return static_cast<MBin*>(fBins->At(b));
 642     }
 643     /**
 644      * Add event observation
 645      *
 646      * @param triggered  IF this is triggered
 647      * @param m          Multiplicity
 648      * @param data       Observation
 649      */
 650     void AddEvent(Bool_t triggered, Bool_t hasVtx, UShort_t m, const TH2D* data)
 651     {
 652       fM->AddBinContent(1);
 653       fM->AddBinContent(TMath::Min(fM->GetNbinsX(), m+2));
 654
 655       MBin* all = static_cast<MBin*>(fBins->At(0));
 656       all->AddEvent(triggered, hasVtx, data);
 657
 658       MBin* bin = FindBin(m);
 659       bin->AddEvent(triggered, hasVtx, data);
 660     }
 661     /**
 662      * End of job processing
 663      *
 664      * @param p Parent list
 665      * @param o Parent output list
 666      */
 667     void Finish(const TList* p, TList* o)
 668     {
 669       TList* l = dynamic_cast<TList*>(p->FindObject(GetName()));
 670       if (!l) {
 671         Warning("Finish", "Cannot find %s in %s", GetName(), p->GetName());
 672         return;
 673       }
 674
 675       TList* ours  = new TList;
 676       ours->SetName(GetName());
 677       ours->SetOwner();
 678       o->Add(ours);
 679
 680       fM = static_cast<TH1D*>(l->FindObject("m"));
 681       if (!fM) {
 682         Warning("Finish", "Didn't find object 'm' in %s", l->GetName());
 683         return;
 684       }
 685       TH1D* m = static_cast<TH1D*>(fM->Clone("m"));
 686       m->SetDirectory(0);
 687       m->SetYTitle("P(N_{ch}|_{|#eta|<1} < X)");
 688       if (m->GetBinContent(1) > 0)
 689         m->Scale(1. / m->GetBinContent(1));
 690       ours->Add(m);
 691
 692       Int_t nBin = fM->GetNbinsX();
 693       TH1D* effs = static_cast<TH1D*>(fM->Clone("effs"));
 694       effs->SetYTitle("#epsilon_{X}");
 695       effs->SetFillColor(kRed+1);
 696       effs->SetDirectory(0);
 697       effs->SetMinimum(0);
 698
 699       gStyle->SetPalette(1);
 700       Int_t   nCol = gStyle->GetNumberOfColors();
 701       THStack* stack = new THStack("biases", "Trigger biases");
 702       for (Int_t i = 0; i <= nBin; i++) {
 703         MBin* bin = static_cast<MBin*>(fBins->At(i));
 704         effs->SetBinContent(i+1, bin->Finish(l, ours, stack));
 705         TH1* h = static_cast<TH1*>(stack->GetHists()->At(i));
 706         Int_t col = kBlack;
 707         if (i != 0) {
 708           Int_t icol = TMath::Min(nCol-1,int(double(i)/nBin * nCol + .5));
 709           col        = gStyle->GetColorPalette(icol);
 710         }
 711         h->SetMarkerColor(col);
 712         h->SetFillColor(col);
 713         h->SetLineColor(col);
 714       }
 715
 716       ours->Add(stack);
 717       ours->Add(effs);
 718     }
 719     UShort_t   fMask;
 720     TH1D*      fM;
 721     TObjArray* fBins;
 722     ClassDef(TriggerType,1);
 723   };
 724   TriggerType            fInel;
 725   TriggerType            fInelGt0;
 726   TriggerType            fNSD;
 727   AliFMDMCEventInspector fInspector;
 728   TList*                 fList;
 729   Bool_t                 fFirstEvent;
 730   TH2D*                  fData;
 731   TH1I*                  fTriggers;
 732   UInt_t                 fTrackletRequirement;
 733   UInt_t                 fVertexRequirement;
 734   TAxis                  fVertexAxis;
 735   TH1D*                  fVertexESD;
 736   TH1D*                  fVertexMC;
 737   TH1D*                  fM;
 738   ClassDef(EvaluateTrigger,1);
 739 };
 740 #else
 741 //====================================================================
 742 void MakeEvaluateTriggers(const char* esddir,
 743                           Int_t       nEvents    = -1,
 744                           UInt_t      vtx        = 0x1,
 745                           UInt_t      trk        = 0x1,
 746                           UInt_t      vz         = 10,
 747                           Int_t       proof      = 0)
 748 {
 749   // --- Libraries to load -------------------------------------------
 750   gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
 751
 752   // --- Check for proof mode, and possibly upload pars --------------
 753   if (proof> 0) {
 754     gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadPars.C");
 755     if (!LoadPars(proof)) {
 756       Error("MakeAOD", "Failed to load PARs");
 757       return;
 758     }
 759   }
 760
 761   // --- Our data chain ----------------------------------------------
 762   gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/MakeChain.C");
 763   TChain* chain = MakeChain("ESD", esddir,true);
 764   // If 0 or less events is select, choose all
 765   if (nEvents <= 0) nEvents = chain->GetEntries();
 766
 767   // --- Set the macro path ------------------------------------------
 768   gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2:"
 769                            "$ALICE_ROOT/ANALYSIS/macros",
 770                            gROOT->GetMacroPath()));
 771
 772   // --- Creating the manager and handlers ---------------------------
 773   AliAnalysisManager *mgr  = new AliAnalysisManager("Triggers",
 774                                                     "Forward multiplicity");
 775   AliAnalysisManager::SetCommonFileName("triggers.root");
 776
 777   // --- ESD input handler -------------------------------------------
 778   AliESDInputHandler *esdHandler = new AliESDInputHandler();
 779   mgr->SetInputEventHandler(esdHandler);
 780
 781   // --- Monte Carlo handler -----------------------------------------
 782   AliMCEventHandler* mcHandler = new AliMCEventHandler();
 783   mgr->SetMCtruthEventHandler(mcHandler);
 784   mcHandler->SetReadTR(true);
 785
 786   // --- Add tasks ---------------------------------------------------
 787   // Physics selection
 788   gROOT->Macro(Form("AddTaskPhysicsSelection.C(1,1,0)"));
 789
 790 #if 0
 791   // --- Fix up physics selection to give proper A,C, and E triggers -
 792   AliInputEventHandler* ih =
 793     static_cast<AliInputEventHandler*>(mgr->GetInputEventHandler());
 794   AliPhysicsSelection* ps =
 795     static_cast<AliPhysicsSelection*>(ih->GetEventSelection());
 796   // Ignore trigger class when selecting events.  This mean that we
 797   // get offline+(A,C,E) events too
 798   ps->SetSkipTriggerClassSelection(true);
 799 #endif
 800
 801   // --- compile our code --------------------------------------------
 802   gSystem->AddIncludePath("-I${ALICE_ROOT}/PWG2/FORWARD/analysis2 "
 803                           "-I${ALICE_ROOT}/ANALYSIS "
 804                           "-I${ALICE_ROOT}/include -DBUILD=1");
 805   gROOT->LoadMacro("${ALICE_ROOT}/PWG2/FORWARD/analysis2/MakeEvaluateTriggers.C++g");
 806
 807   // --- Make our object ---------------------------------------------
 808   EvaluateTrigger* task = new EvaluateTrigger("triggers");
 809   mgr->AddTask(task);
 810   task->SetVertexRequirement(vtx);
 811   task->SetTrackletRequirement(trk);
 812   task->SetVertexAxis(10, -vz, vz);
 813
 814   // --- create containers for input/output --------------------------
 815   AliAnalysisDataContainer *sums =
 816     mgr->CreateContainer("triggerSums", TList::Class(),
 817                          AliAnalysisManager::kOutputContainer,
 818                          AliAnalysisManager::GetCommonFileName());
 819   AliAnalysisDataContainer *output =
 820     mgr->CreateContainer("triggerResults", TList::Class(),
 821                          AliAnalysisManager::kParamContainer,
 822                          AliAnalysisManager::GetCommonFileName());
 823
 824   // --- connect input/output ----------------------------------------
 825   mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
 826   mgr->ConnectOutput(task, 1, sums);
 827   mgr->ConnectOutput(task, 2, output);
 828
 829   // --- Run the analysis --------------------------------------------
 830   TStopwatch t;
 831   if (!mgr->InitAnalysis()) {
 832     Error("MakeAOD", "Failed to initialize analysis train!");
 833     return;
 834   }
 835   // Skip terminate if we're so requested and not in Proof or full mode
 836   mgr->SetSkipTerminate(false);
 837   // Some informative output
 838   mgr->PrintStatus();
 839   if (proof) mgr->SetDebugLevel(3);
 840   if (mgr->GetDebugLevel() < 1 && !proof)
 841     mgr->SetUseProgressBar(kTRUE,100);
 842
 843   // Run the train
 844   t.Start();
 845   Printf("=== RUNNING ANALYSIS on %9d events ==========================",
 846          nEvents);
 847   mgr->StartAnalysis(proof ? "proof" : "local", chain, nEvents);
 848   t.Stop();
 849   t.Print();
 850 }
 851 #endif