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