PWGLF/FORWARD/analysis2/AliAODForwardMult.h

   1 //
   2 // See implementation or Doxygen comments for more information
   3 //
   4 #ifndef ALIAODFORWARDMULT_H
   5 #define ALIAODFORWARDMULT_H
   6 /**
   7  * @file   AliAODForwardMult.h
   8  * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
   9  * @date   Wed Mar 23 13:58:00 2011
  10  *
  11  * @brief  Per-event @f$ N_{ch}@f$ per @f$(\eta,\varphi)@f$ bin
  12  *
  13  * @ingroup pwglf_forward_aod
  14  *
  15  */
  16 #include <TObject.h>
  17 #include <TH2D.h>
  18 class TBrowser;
  19 class TH1I;
  20 /**
  21  * Class that contains the forward multiplicity data per event
  22  *
  23  * This class contains a histogram of
  24  * @f[
  25  *   \frac{d^2N_{ch}}{d\eta d\phi}\quad,
  26  * @f]
  27  * as well as a trigger mask for each analysed event.
  28  *
  29  * The eta acceptance of the event is stored in the underflow bins of
  30  * the histogram.  So to build the final histogram, one needs to
  31  * correct for this acceptance (properly weighted by the events), and
  32  * the vertex efficiency.  This simply boils down to defining a 2D
  33  * histogram and summing the event histograms in that histogram.  One
  34  * should of course also do proper book-keeping of the accepted event.
  35  *
  36  * @code
  37  * TTree* GetAODTree()
  38  * {
  39  *    TFile* file = TFile::Open("AliAODs.root","READ");
  40  *    TTree* tree = static_cast<TTree*>(file->Get("aodTree"));
  41  *    return tree;
  42  * }
  43  *
  44  * void Analyse()
  45  * {
  46  *   TH2D*              sum        = 0;                  // Summed hist
  47  *   TTree*             tree       = GetAODTree();       // AOD tree
  48  *   AliAODForwardMult* mult       = 0;                  // AOD object
  49  *   Int_t              nTriggered = 0;                  // # of triggered ev.
  50  *   Int_t              nWithVertex= 0;                  // # of ev. w/vertex
  51  *   Int_t              nAccepted  = 0;                  // # of ev. used
  52  *   Int_t              nAvailable = tree->GetEntries(); // How many entries
  53  *   Float_t            vzLow      = -10;                // Lower ip cut
  54  *   Float_t            vzHigh     =  10;                // Upper ip cut
  55  *   Int_t              mask       = AliAODForwardMult::kInel;// Trigger mask
  56  *   tree->SetBranchAddress("forward", &forward);        // Set the address
  57  *
  58  *   for (int i = 0; i < nAvailable; i++) {
  59  *     // Create sum histogram on first event - to match binning to input
  60  *     if (!sum) sum = static_cast<TH2D*>(mult->Clone("d2ndetadphi"));
  61  *
  62  *     tree->GetEntry(i);
  63  *
  64  *     // Other trigger/event requirements could be defined
  65  *     if (!mult->IsTriggerBits(mask)) continue;
  66  *     nTriggered++;
  67  *
  68  *     // Check if we have vertex
  69  *     if (!mult->HasIpZ()) continue;
  70  *     nWithVertex++;
  71  *
  72  *     // Select vertex range (in centimeters)
  73  *     if (!mult->InRange(vzLow, vzHigh) continue;
  74  *     nAccepted++;
  75  *
  76  *     // Add contribution from this event
  77  *     sum->Add(&(mult->GetHistogram()));
  78  *   }
  79  *
  80  *   // Get acceptance normalisation from underflow bins
  81  *   TH1D* norm   = sum->ProjectionX("norm", 0, 1, "");
  82  *   // Project onto eta axis - _ignoring_underflow_bins_!
  83  *   TH1D* dndeta = sum->Projection("dndeta", 1, -1, "e");
  84  *   // Normalize to the acceptance
  85  *   dndeta->Divide(norm);
  86  *   // Scale by the vertex efficiency
  87  *   dndeta->Scale(Double_t(nWithVertex)/nTriggered, "width");
  88  *   // And draw the result
  89  *   dndeta->Draw();
  90  * }
  91  * @endcode
  92  *
  93  * The above code will draw the final @f$ dN_{ch}/d\eta@f$ for the
  94  * selected event class and vertex range
  95  *
  96  * The histogram can be used as input for other kinds of analysis too,
  97  * like flow, event-plane, centrality, and so on.
  98  *
  99  * @ingroup pwglf_forward
 100  * @ingroup pwglf_forward_aod
 101  */
 102 class AliAODForwardMult : public TObject
 103 {
 104 public:
 105   /**
 106    * Bits of the trigger pattern
 107    */
 108   enum {
 109     /** In-elastic collision - really MBOR*/
 110     kInel        = 0x0001,
 111     /** In-elastic collision with at least one SPD tracklet */
 112     kInelGt0     = 0x0002,
 113     /** Non-single diffractive collision - (V0AND||FASTOR>5) */
 114     kNSD         = 0x0004,
 115     /** Empty bunch crossing */
 116     kEmpty       = 0x0008,
 117     /** A-side trigger */
 118     kA           = 0x0010,
 119     /** B(arrel) trigger */
 120     kB           = 0x0020,
 121     /** C-side trigger */
 122     kC           = 0x0080,
 123     /** Empty trigger */
 124     kE           = 0x0100,
 125     /** pileup from SPD */
 126     kPileUp      = 0x0200,
 127     /** true NSD from MC */
 128     kMCNSD       = 0x0400,
 129     /** Offline MB triggered */
 130     kOffline     = 0x0800,
 131     /** At least one SPD cluster */
 132     kNClusterGt0 = 0x1000,
 133     /** V0-AND trigger */
 134     kV0AND       = 0x2000,
 135     /** Satellite event */
 136     kSatellite   = 0x4000,
 137     /** SPD outlier event */
 138     kSPDOutlier  = 0x8000,
 139     /** SPD pile-up */
 140     kPileupSPD   = 0x10000,
 141     /** Track pile-up */
 142     kPileupTrack = 0x20000,
 143     /** Out of bunch pileup */
 144     kPileupBC    = 0x40000
 145
 146   };
 147   /**
 148    * Bin numbers in trigger histograms
 149    */
 150   enum {
 151     kBinAll=1,
 152     kBinInel,
 153     kBinInelGt0,
 154     kBinNSD,
 155     kBinV0AND,
 156     kBinA,
 157     kBinB,
 158     kBinC,
 159     kBinE,
 160     kBinSatellite,
 161     kBinPileUp,
 162     kBinMCNSD,
 163     kBinOffline,
 164     kBinNClusterGt0,
 165     kWithTrigger,
 166     kWithVertex,
 167     kAccepted
 168   };
 169   /**
 170    * User bits of these objects (bits 14-23 can be used)
 171    */
 172   enum {
 173     /** Secondary correction maps where applied */
 174     kSecondary           = (1 << 14),
 175     /** Vertex bias correction was applied */
 176     kVertexBias          = (1 << 15),
 177     /** Acceptance correction was applied */
 178     kAcceptance          = (1 << 16),
 179     /** Merging efficiency correction was applied */
 180     kMergingEfficiency   = (1 << 17),
 181     /** Signal in overlaps is the sum */
 182     kSum                 = (1 << 18),
 183     /** Used eta dependent empirical correction - to be implemented */
 184     kEmpirical           = (1 << 19)
 185   };
 186   /**
 187    * Return codes of CheckEvent
 188    */
 189   enum ECheckStatus {
 190     /** Event accepted by cuts */
 191     kGoodEvent = 0,
 192     /** Event centrality not in range */
 193     kWrongCentrality,
 194     /** Event trigger isn't in the supplied mask */
 195     kWrongTrigger,
 196     /** Event is a pile-up event */
 197     kIsPileup,
 198     /** Event has no interaction point information */
 199     kNoVertex,
 200     /** Event interaction point is out of range */
 201     kWrongVertex,
 202     /** Outlier */
 203     kOutlierEvent
 204   };
 205
 206   /**
 207    * Default constructor
 208    *
 209    * Used by ROOT I/O sub-system - do not use
 210    */
 211   AliAODForwardMult();
 212   /**
 213    * Constructor
 214    *
 215    * @param isMC Whether this was from MC or not
 216    */
 217   AliAODForwardMult(Bool_t isMC);
 218   /**
 219    * Destructor
 220    */
 221   virtual ~AliAODForwardMult() {} // Destructor
 222   /**
 223    * Initialize
 224    *
 225    * @param etaAxis  Pseudo-rapidity axis
 226    */
 227   void Init(const TAxis& etaAxis);
 228   /**
 229    * Get the @f$ d^2N_{ch}/d\eta d\phi@f$ histogram,
 230    *
 231    * @return @f$ d^2N_{ch}/d\eta d\phi@f$ histogram,
 232    */
 233   const TH2D& GetHistogram() const { return fHist; } // Get histogram
 234   /**
 235    * Get the @f$ d^2N_{ch}/d\eta d\phi@f$ histogram,
 236    *
 237    * @return @f$ d^2N_{ch}/d\eta d\phi@f$ histogram,
 238    */
 239   TH2D& GetHistogram() { return fHist; } // Get histogram
 240   /**
 241    * Get the trigger bits
 242    *
 243    * @return Trigger bits
 244    */
 245   UInt_t GetTriggerBits() const { return fTriggers; } // Get triggers
 246   /**
 247    * Set the trigger mask
 248    *
 249    * @param trg Trigger mask
 250    */
 251   void SetTriggerMask(UInt_t trg) { fTriggers = trg; } // Set triggers
 252   /**
 253    * Set bit(s) in trigger mask
 254    *
 255    * @param bits bit(s) to set
 256    */
 257   void SetTriggerBits(UInt_t bits) { fTriggers |= bits; } // Set trigger bits
 258   /**
 259    * Check if all bit(s) are set in the trigger mask.  Note, this is
 260    * an @e and between the bits.  If you need an @e or you should use
 261    * the member function IsTriggerOrBits
 262    *
 263    * @param bits Bits to test for
 264    *
 265    * @return true if all enabled bits in the argument is also set in
 266    * the trigger word
 267    */
 268   Bool_t IsTriggerBits(UInt_t bits) const;
 269   /**
 270    * Check if any of bit(s) are enabled in the trigger word.  This is
 271    * an @e or between the selected bits.  If you need and @a and you
 272    * should use the member function IsTriggerBits;
 273    *
 274    * @param bits Bits to check for
 275    *
 276    * @return true if any of the enabled bits in the arguments are also
 277    * enabled in the trigger mask
 278    */
 279   Bool_t IsTriggerOrBits(UInt_t bits) const;
 280   /**
 281    * Whether we have any trigger bits
 282    *
 283    * @return true if we have some trigger
 284    */
 285   Bool_t HasTrigger() const { return fTriggers != 0; } // Check for triggers
 286   /**
 287    * Clear all data
 288    *
 289    * @param option  Passed on to TH2::Reset verbatim
 290    */
 291   void Clear(Option_t* option="");
 292   /**
 293    * browse this object
 294    *
 295    * @param b Browser
 296    */
 297   void Browse(TBrowser* b);
 298   /**
 299    * This is a folder
 300    *
 301    * @return Always true
 302    */
 303   Bool_t IsFolder() const { return kTRUE; } // Always true
 304   /**
 305    * Check if the data has been secondary corrected by MC maps
 306    *
 307    * @return true if secondary corrected via MC maps
 308    */
 309   Bool_t IsSecondaryCorrected() const { return TestBit(kSecondary); }
 310   /**
 311    * Check if vertex bias correction was applied
 312    *
 313    * @return true if MC based vertex bias was applied
 314    */
 315   Bool_t IsVertexBiasCorrected() const { return TestBit(kVertexBias); }
 316   /**
 317    * Check if acceptance correction (from dead strips) was applied
 318    *
 319    * @return true if acceptance correction was applied
 320    */
 321   Bool_t IsAcceptanceCorrected() const { return TestBit(kAcceptance); }
 322   /**
 323    * Check if merging efficiency (from MC) was applied
 324    *
 325    * @return true if merging efficiency was applied
 326    */
 327   Bool_t IsMergingEfficiencyCorrected() const {
 328     return TestBit(kMergingEfficiency); }
 329   /**
 330    * Check if an empirical correction was applied in the event level.
 331    *
 332    * @return True if the empirical correction was applied per event.
 333    */
 334   Bool_t IsEmpiricalCorrected() const { return TestBit(kEmpirical); }
 335   /**
 336    * Check if the output is the sum (not average) in regions of
 337    * overlap between detectors.
 338    *
 339    *
 340    * @return true if the sum (not average) is stored in overlap
 341    * regions.
 342    */
 343   Bool_t IsSumSignal() const { return TestBit(kSum); }
 344   /**
 345    * Print content
 346    *
 347    * @param option Passed verbatim to TH2::Print
 348    */
 349   void Print(Option_t* option="") const;
 350   /**
 351    * Set the z coordinate of the interaction point
 352    *
 353    * @param ipZ Interaction point z coordinate
 354    */
 355   void SetIpZ(Float_t ipZ) { fIpZ = ipZ; } // Set Ip's Z coordinate
 356   /**
 357    * Set the center of mass energy per nucleon-pair.  This is stored
 358    * in the (0,0) of the histogram
 359    *
 360    * @param sNN Center of mass energy per nucleon pair (GeV)
 361    */
 362   void SetSNN(UShort_t sNN);
 363   /**
 364    * Get the collision system number
 365    * - 0: Unknown
 366    * - 1: pp
 367    * - 2: PbPb
 368    *
 369    * @param sys Collision system number
 370    */
 371   void SetSystem(UShort_t sys);
 372   /**
 373    * Set the event centrality
 374    *
 375    * @param c Centrality
 376    */
 377   void SetCentrality(Float_t c) { fCentrality = c; }
 378   /**
 379    * Set the z coordinate of the interaction point
 380    *
 381    * @return Interaction point z coordinate
 382    */
 383   Float_t GetIpZ() const { return fIpZ; } // Get Ip's Z coordinate
 384   /**
 385    * Check if we have a valid z coordinate of the interaction point
 386    *
 387    * @return True if we have a valid interaction point z coordinate
 388    */
 389   Bool_t HasIpZ() const;
 390   /**
 391    * Get the center of mass energy per nucleon pair (GeV)
 392    *
 393    * @return Center of mass energy per nucleon pair (GeV)
 394    */
 395   UShort_t GetSNN() const;
 396   /**
 397    * Get the collision system number
 398    * - 0: Unknown
 399    * - 1: pp
 400    * - 2: PbPb
 401    *
 402    * @return Collision system number
 403    */
 404   UShort_t GetSystem() const;
 405   /**
 406    * Check if the z coordinate of the interaction point is within the
 407    * given limits.  Note that the convention used corresponds to the
 408    * convention used in ROOTs TAxis.
 409    *
 410    * @param low  Lower cut (inclusive)
 411    * @param high Upper cut (exclusive)
 412    *
 413    * @return true if @f$ low \ge ipz < high@f$
 414    */
 415   Bool_t InRange(Float_t low, Float_t high) const;
 416   /**
 417    * Get the event centrality
 418    *
 419    * @return Centrality
 420    */
 421   Float_t GetCentrality() const { return fCentrality; }
 422   /**
 423    * Check if we have a valid centrality
 424    *
 425    * @return true if the centrality information is valid
 426    */
 427   Bool_t  HasCentrality() const { return !(fCentrality  < 0); }
 428   /**
 429    * Get the number of SPD clusters seen in @f$ |\eta|<1@f$
 430    *
 431    * @return Number of SPD clusters seen
 432    */
 433   UShort_t GetNClusters() const { return fNClusters; }
 434   /**
 435    * Set the number of SPD clusters seen in @f$ |\eta|<1@f$
 436    *
 437    * @param n Number of SPD clusters
 438    */
 439   void SetNClusters(UShort_t n) { fNClusters = n; }
 440   /**
 441    * Get the name of the object
 442    *
 443    * @return Name of object
 444    */
 445   const Char_t* GetName() const { return (fIsMC ? "ForwardMC" : "Forward"); }
 446   /**
 447    * Check if event meets the passses requirements.
 448    *
 449    * It returns true if @e all of the following is true
 450    *
 451    * - The trigger is within the bit mask passed.
 452    * - The vertex is within the specified limits.
 453    * - The centrality is within the specified limits, or if lower
 454    *   limit is equal to or larger than the upper limit.
 455    *
 456    * Note, for data with out a centrality estimate (e.g., pp), one
 457    * must pass equal centrality cuts, or no data will be accepted.  In
 458    * other words, for pp data, always pass cMin=0, cMax=0
 459    *
 460    * If a histogram is passed in the last parameter, then that
 461    * histogram is filled with the trigger bits.
 462    *
 463    * @param triggerMask  Trigger mask
 464    * @param vzMin        Minimum @f$ v_z@f$ (in centimeters)
 465    * @param vzMax        Maximum @f$ v_z@f$ (in centimeters)
 466    * @param cMin         Minimum centrality (in percent)
 467    * @param cMax         Maximum centrality (in percent)
 468    * @param hist         Histogram to fill
 469    * @param status       Histogram to fill
 470    * @param removePileup If true, do not accept pile-up events (default)
 471    *
 472    * @return @c true if the event meets the requirements
 473    */
 474   Bool_t CheckEvent(Int_t    triggerMask=kInel,
 475                     Double_t vzMin=-10, Double_t vzMax=10,
 476                     UShort_t cMin=0,    UShort_t cMax=100,
 477                     TH1*     hist=0,
 478                     TH1*     status=0,
 479                     Bool_t   removePileup=true) const;
 480   /**
 481    * Get a string correspondig to the trigger mask
 482    *
 483    * @param mask Trigger mask
 484    *
 485    * @return Static string (copy before use)
 486    */
 487   static const Char_t* GetTriggerString(UInt_t mask);
 488   /**
 489    * Make a histogram to record triggers in.
 490    *
 491    * The bins defined by the trigger enumeration in this class.  One
 492    * can use this enumeration to retrieve the number of triggers for
 493    * each class.
 494    *
 495    * @param name Name of the histogram
 496    * @param mask Trigger mask
 497    *
 498    * @return Newly allocated histogram
 499    */
 500   static TH1I* MakeTriggerHistogram(const char* name="triggers",
 501                                     Int_t mask=0);
 502   /**
 503    * Make a histogram to record status in.
 504    *
 505    * The bins defined by the status enumeration in this class.
 506    *
 507    * @param name Name of the histogram
 508    *
 509    * @return Newly allocated histogram
 510    */
 511   static TH1I* MakeStatusHistogram(const char* name="status");
 512   /**
 513    * Utility function to make a trigger mask from the passed string.
 514    *
 515    * The string is a comma or space seperated list of case-insensitive
 516    * strings
 517    *
 518    * - INEL
 519    * - INEL>0
 520    * - NSD
 521    *
 522    * @param what Which triggers to put in the mask.
 523    *
 524    * @return The generated trigger mask.
 525    */
 526   static UInt_t MakeTriggerMask(const char* what);
 527 protected:
 528   /** From MC or not */
 529   Bool_t   fIsMC;       // Whether this is from MC
 530   /** Histogram of @f$d^2N_{ch}/(d\eta d\phi)@f$ for this event */
 531   TH2D     fHist;       // Histogram of d^2N_{ch}/(deta dphi) for this event
 532   /** Trigger bits */
 533   UInt_t   fTriggers;   // Trigger bit mask
 534   /** Interaction point @f$z@f$ coordinate */
 535   Float_t  fIpZ;        // Z coordinate of the interaction point
 536   /** Centrality */
 537   Float_t  fCentrality; // Event centrality
 538   /** Number of clusters in @f$|\eta|<1@f$ */
 539   UShort_t fNClusters;  // Number of SPD clusters in |eta|<1
 540   /** Invalid value for interaction point @f$z@f$ coordiante */
 541   static const Float_t fgkInvalidIpZ; // Invalid IpZ value
 542   ClassDef(AliAODForwardMult,5); // AOD forward multiplicity
 543 };
 544
 545 //____________________________________________________________________
 546 inline Bool_t
 547 AliAODForwardMult::InRange(Float_t low, Float_t high) const
 548 {
 549   return HasIpZ() && fIpZ >= low && fIpZ < high;
 550 }
 551
 552 //____________________________________________________________________
 553 inline Bool_t
 554 AliAODForwardMult::IsTriggerBits(UInt_t bits) const
 555 {
 556   return HasTrigger() && ((fTriggers & bits) == bits);
 557 }
 558 //____________________________________________________________________
 559 inline Bool_t
 560 AliAODForwardMult::IsTriggerOrBits(UInt_t bits) const
 561 {
 562   return HasTrigger() && ((fTriggers & bits) != 0);
 563 }
 564
 565 #endif
 566 // Local Variables:
 567 //  mode: C++
 568 // End:
 569