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
  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 */
 110     kInel        = 0x0001,
 111     /** In-elastic collision with at least one SPD tracklet */
 112     kInelGt0     = 0x0002,
 113     /** Non-single diffractive collision */
 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   };
 138   /**
 139    * Bin numbers in trigger histograms
 140    */
 141   enum {
 142     kBinAll=1,
 143     kBinInel,
 144     kBinInelGt0,
 145     kBinNSD,
 146     kBinV0AND,
 147     kBinA,
 148     kBinB,
 149     kBinC,
 150     kBinE,
 151     kBinSatellite,
 152     kBinPileUp,
 153     kBinMCNSD,
 154     kBinOffline,
 155     kBinNClusterGt0,
 156     kWithTrigger,
 157     kWithVertex,
 158     kAccepted
 159   };
 160   /**
 161    * User bits of these objects (bits 14-23 can be used)
 162    */
 163   enum {
 164     /** Secondary correction maps where applied */
 165     kSecondary           = (1 << 14),
 166     /** Vertex bias correction was applied */
 167     kVertexBias          = (1 << 15),
 168     /** Acceptance correction was applied */
 169     kAcceptance          = (1 << 16),
 170     /** Merging efficiency correction was applied */
 171     kMergingEfficiency   = (1 << 17),
 172     /** Signal in overlaps is the sum */
 173     kSum                 = (1 << 18),
 174     /** Used eta dependent empirical correction - to be implemented */
 175     kEmpirical           = (1 << 19)
 176   };
 177   /**
 178    * Return codes of CheckEvent
 179    */
 180   enum ECheckStatus {
 181     /** Event accepted by cuts */
 182     kGoodEvent = 0,
 183     /** Event centrality not in range */
 184     kWrongCentrality,
 185     /** Event trigger isn't in the supplied mask */
 186     kWrongTrigger,
 187     /** Event is a pile-up event */
 188     kIsPileup,
 189     /** Event has no interaction point information */
 190     kNoVertex,
 191     /** Event interaction point is out of range */
 192     kWrongVertex
 193   };
 194
 195   /**
 196    * Default constructor
 197    *
 198    * Used by ROOT I/O sub-system - do not use
 199    */
 200   AliAODForwardMult();
 201   /**
 202    * Constructor
 203    *
 204    * @param isMC Whether this was from MC or not
 205    */
 206   AliAODForwardMult(Bool_t isMC);
 207   /**
 208    * Destructor
 209    */
 210   virtual ~AliAODForwardMult() {} // Destructor
 211   /**
 212    * Initialize
 213    *
 214    * @param etaAxis  Pseudo-rapidity axis
 215    */
 216   void Init(const TAxis& etaAxis);
 217   /**
 218    * Get the @f$ d^2N_{ch}/d\eta d\phi@f$ histogram,
 219    *
 220    * @return @f$ d^2N_{ch}/d\eta d\phi@f$ histogram,
 221    */
 222   const TH2D& GetHistogram() const { return fHist; } // Get histogram
 223   /**
 224    * Get the @f$ d^2N_{ch}/d\eta d\phi@f$ histogram,
 225    *
 226    * @return @f$ d^2N_{ch}/d\eta d\phi@f$ histogram,
 227    */
 228   TH2D& GetHistogram() { return fHist; } // Get histogram
 229   /**
 230    * Get the trigger bits
 231    *
 232    * @return Trigger bits
 233    */
 234   UInt_t GetTriggerBits() const { return fTriggers; } // Get triggers
 235   /**
 236    * Set the trigger mask
 237    *
 238    * @param trg Trigger mask
 239    */
 240   void SetTriggerMask(UInt_t trg) { fTriggers = trg; } // Set triggers
 241   /**
 242    * Set bit(s) in trigger mask
 243    *
 244    * @param bits bit(s) to set
 245    */
 246   void SetTriggerBits(UInt_t bits) { fTriggers |= bits; } // Set trigger bits
 247   /**
 248    * Check if all bit(s) are set in the trigger mask.  Note, this is
 249    * an @e and between the bits.  If you need an @e or you should use
 250    * the member function IsTriggerOrBits
 251    *
 252    * @param bits Bits to test for
 253    *
 254    * @return true if all enabled bits in the argument is also set in
 255    * the trigger word
 256    */
 257   Bool_t IsTriggerBits(UInt_t bits) const;
 258   /**
 259    * Check if any of bit(s) are enabled in the trigger word.  This is
 260    * an @e or between the selected bits.  If you need and @a and you
 261    * should use the member function IsTriggerBits;
 262    *
 263    * @param bits Bits to check for
 264    *
 265    * @return true if any of the enabled bits in the arguments are also
 266    * enabled in the trigger mask
 267    */
 268   Bool_t IsTriggerOrBits(UInt_t bits) const;
 269   /**
 270    * Whether we have any trigger bits
 271    *
 272    * @return true if we have some trigger
 273    */
 274   Bool_t HasTrigger() const { return fTriggers != 0; } // Check for triggers
 275   /**
 276    * Clear all data
 277    *
 278    * @param option  Passed on to TH2::Reset verbatim
 279    */
 280   void Clear(Option_t* option="");
 281   /**
 282    * browse this object
 283    *
 284    * @param b Browser
 285    */
 286   void Browse(TBrowser* b);
 287   /**
 288    * This is a folder
 289    *
 290    * @return Always true
 291    */
 292   Bool_t IsFolder() const { return kTRUE; } // Always true
 293
 294   Bool_t IsSecondaryCorrected() const { return TestBit(kSecondary); }
 295   Bool_t IsVertexBiasCorrected() const { return TestBit(kVertexBias); }
 296   Bool_t IsAcceptanceCorrected() const { return TestBit(kAcceptance); }
 297   Bool_t IsMergingEfficiencyCorrected() const {
 298     return TestBit(kMergingEfficiency); }
 299   Bool_t IsEmpiricalCorrected() const { return TestBit(kEmpirical); }
 300   Bool_t IsSumSignal() const { return TestBit(kSum); }
 301   /**
 302    * Print content
 303    *
 304    * @param option Passed verbatim to TH2::Print
 305    */
 306   void Print(Option_t* option="") const;
 307   /**
 308    * Set the z coordinate of the interaction point
 309    *
 310    * @param ipZ Interaction point z coordinate
 311    */
 312   void SetIpZ(Float_t ipZ) { fIpZ = ipZ; } // Set Ip's Z coordinate
 313   /**
 314    * Set the center of mass energy per nucleon-pair.  This is stored
 315    * in the (0,0) of the histogram
 316    *
 317    * @param sNN Center of mass energy per nucleon pair (GeV)
 318    */
 319   void SetSNN(UShort_t sNN);
 320   /**
 321    * Get the collision system number
 322    * - 0: Unknown
 323    * - 1: pp
 324    * - 2: PbPb
 325    *
 326    * @param sys Collision system number
 327    */
 328   void SetSystem(UShort_t sys);
 329   /**
 330    * Set the event centrality
 331    *
 332    * @param c Centrality
 333    */
 334   void SetCentrality(Float_t c) { fCentrality = c; }
 335   /**
 336    * Set the z coordinate of the interaction point
 337    *
 338    * @return Interaction point z coordinate
 339    */
 340   Float_t GetIpZ() const { return fIpZ; } // Get Ip's Z coordinate
 341   /**
 342    * Check if we have a valid z coordinate of the interaction point
 343    *
 344    * @return True if we have a valid interaction point z coordinate
 345    */
 346   Bool_t HasIpZ() const;
 347   /**
 348    * Get the center of mass energy per nucleon pair (GeV)
 349    *
 350    * @return Center of mass energy per nucleon pair (GeV)
 351    */
 352   UShort_t GetSNN() const;
 353   /**
 354    * Get the collision system number
 355    * - 0: Unknown
 356    * - 1: pp
 357    * - 2: PbPb
 358    *
 359    * @return Collision system number
 360    */
 361   UShort_t GetSystem() const;
 362   /**
 363    * Check if the z coordinate of the interaction point is within the
 364    * given limits.  Note that the convention used corresponds to the
 365    * convention used in ROOTs TAxis.
 366    *
 367    * @param low  Lower cut (inclusive)
 368    * @param high Upper cut (exclusive)
 369    *
 370    * @return true if @f$ low \ge ipz < high@f$
 371    */
 372   Bool_t InRange(Float_t low, Float_t high) const;
 373   /**
 374    * Get the event centrality
 375    *
 376    *
 377    * @return
 378    */
 379   Float_t GetCentrality() const { return fCentrality; }
 380   /**
 381    * Check if we have a valid centrality
 382    *
 383    *
 384    * @return
 385    */
 386   Bool_t  HasCentrality() const { return !(fCentrality  < 0); }
 387   /**
 388    * Get the number of SPD clusters seen in @f$ |\eta|<1@f$
 389    *
 390    * @return Number of SPD clusters seen
 391    */
 392   UShort_t GetNClusters() const { return fNClusters; }
 393   /**
 394    * Set the number of SPD clusters seen in @f$ |\eta|<1@f$
 395    *
 396    * @param n Number of SPD clusters
 397    */
 398   void SetNClusters(UShort_t n) { fNClusters = n; }
 399   /**
 400    * Get the name of the object
 401    *
 402    * @return Name of object
 403    */
 404   const Char_t* GetName() const { return (fIsMC ? "ForwardMC" : "Forward"); }
 405   /**
 406    * Check if event meets the passses requirements.
 407    *
 408    * It returns true if @e all of the following is true
 409    *
 410    * - The trigger is within the bit mask passed.
 411    * - The vertex is within the specified limits.
 412    * - The centrality is within the specified limits, or if lower
 413    *   limit is equal to or larger than the upper limit.
 414    *
 415    * Note, for data with out a centrality estimate (e.g., pp), one
 416    * must pass equal centrality cuts, or no data will be accepted.  In
 417    * other words, for pp data, always pass cMin=0, cMax=0
 418    *
 419    * If a histogram is passed in the last parameter, then that
 420    * histogram is filled with the trigger bits.
 421    *
 422    * @param triggerMask  Trigger mask
 423    * @param vzMin        Minimum @f$ v_z@f$ (in centimeters)
 424    * @param vzMax        Maximum @f$ v_z@f$ (in centimeters)
 425    * @param cMin         Minimum centrality (in percent)
 426    * @param cMax         Maximum centrality (in percent)
 427    * @param hist         Histogram to fill
 428    * @param status       Histogram to fill
 429    *
 430    *
 431    * @return @c true if the event meets the requirements
 432    */
 433   Bool_t CheckEvent(Int_t    triggerMask=kInel,
 434                     Double_t vzMin=-10, Double_t vzMax=10,
 435                     UShort_t cMin=0,    UShort_t cMax=100,
 436                     TH1*     hist=0,
 437                     TH1*     status=0) const;
 438   /**
 439    * Get a string correspondig to the trigger mask
 440    *
 441    * @param mask Trigger mask
 442    *
 443    * @return Static string (copy before use)
 444    */
 445   static const Char_t* GetTriggerString(UInt_t mask);
 446   /**
 447    * Make a histogram to record triggers in.
 448    *
 449    * The bins defined by the trigger enumeration in this class.  One
 450    * can use this enumeration to retrieve the number of triggers for
 451    * each class.
 452    *
 453    * @param name Name of the histogram
 454    * @param mask Trigger mask
 455    *
 456    * @return Newly allocated histogram
 457    */
 458   static TH1I* MakeTriggerHistogram(const char* name="triggers",
 459                                     Int_t mask=0);
 460   /**
 461    * Make a histogram to record status in.
 462    *
 463    * The bins defined by the status enumeration in this class.
 464    *
 465    * @param name Name of the histogram
 466    *
 467    * @return Newly allocated histogram
 468    */
 469   static TH1I* MakeStatusHistogram(const char* name="status");
 470   /**
 471    * Utility function to make a trigger mask from the passed string.
 472    *
 473    * The string is a comma or space seperated list of case-insensitive
 474    * strings
 475    *
 476    * - INEL
 477    * - INEL>0
 478    * - NSD
 479    *
 480    * @param what Which triggers to put in the mask.
 481    *
 482    * @return The generated trigger mask.
 483    */
 484   static UInt_t MakeTriggerMask(const char* what);
 485 protected:
 486   /** From MC or not */
 487   Bool_t   fIsMC;       // Whether this is from MC
 488   /** Histogram of @f$d^2N_{ch}/(d\eta d\phi)@f$ for this event */
 489   TH2D     fHist;       // Histogram of d^2N_{ch}/(deta dphi) for this event
 490   /** Trigger bits */
 491   UInt_t   fTriggers;   // Trigger bit mask
 492   /** Interaction point @f$z@f$ coordinate */
 493   Float_t  fIpZ;        // Z coordinate of the interaction point
 494   /** Centrality */
 495   Float_t  fCentrality; // Event centrality
 496   /** Number of clusters in @f$|\eta|<1@f$ */
 497   UShort_t fNClusters;  // Number of SPD clusters in |eta|<1
 498   /** Invalid value for interaction point @f$z@f$ coordiante */
 499   static const Float_t fgkInvalidIpZ; // Invalid IpZ value
 500   ClassDef(AliAODForwardMult,5); // AOD forward multiplicity
 501 };
 502
 503 //____________________________________________________________________
 504 inline Bool_t
 505 AliAODForwardMult::InRange(Float_t low, Float_t high) const
 506 {
 507   return HasIpZ() && fIpZ >= low && fIpZ < high;
 508 }
 509
 510 //____________________________________________________________________
 511 inline Bool_t
 512 AliAODForwardMult::IsTriggerBits(UInt_t bits) const
 513 {
 514   return HasTrigger() && ((fTriggers & bits) == bits);
 515 }
 516 //____________________________________________________________________
 517 inline Bool_t
 518 AliAODForwardMult::IsTriggerOrBits(UInt_t bits) const
 519 {
 520   return HasTrigger() && ((fTriggers & bits) != 0);
 521 }
 522
 523 #endif
 524 // Local Variables:
 525 //  mode: C++
 526 // End:
 527