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    * Default constructor
 162    *
 163    * Used by ROOT I/O sub-system - do not use
 164    */
 165   AliAODForwardMult();
 166   /**
 167    * Constructor
 168    *
 169    * @param isMC Whether this was from MC or not
 170    */
 171   AliAODForwardMult(Bool_t isMC);
 172   /**
 173    * Destructor
 174    */
 175   virtual ~AliAODForwardMult() {} // Destructor
 176   /**
 177    * Initialize
 178    *
 179    * @param etaAxis  Pseudo-rapidity axis
 180    */
 181   void Init(const TAxis& etaAxis);
 182   /**
 183    * Get the @f$ d^2N_{ch}/d\eta d\phi@f$ histogram,
 184    *
 185    * @return @f$ d^2N_{ch}/d\eta d\phi@f$ histogram,
 186    */
 187   const TH2D& GetHistogram() const { return fHist; } // Get histogram
 188   /**
 189    * Get the @f$ d^2N_{ch}/d\eta d\phi@f$ histogram,
 190    *
 191    * @return @f$ d^2N_{ch}/d\eta d\phi@f$ histogram,
 192    */
 193   TH2D& GetHistogram() { return fHist; } // Get histogram
 194   /**
 195    * Get the trigger bits
 196    *
 197    * @return Trigger bits
 198    */
 199   UInt_t GetTriggerBits() const { return fTriggers; } // Get triggers
 200   /**
 201    * Set the trigger mask
 202    *
 203    * @param trg Trigger mask
 204    */
 205   void SetTriggerMask(UInt_t trg) { fTriggers = trg; } // Set triggers
 206   /**
 207    * Set bit(s) in trigger mask
 208    *
 209    * @param bits bit(s) to set
 210    */
 211   void SetTriggerBits(UInt_t bits) { fTriggers |= bits; } // Set trigger bits
 212   /**
 213    * Check if all bit(s) are set in the trigger mask.  Note, this is
 214    * an @e and between the bits.  If you need an @e or you should use
 215    * the member function IsTriggerOrBits
 216    *
 217    * @param bits Bits to test for
 218    *
 219    * @return true if all enabled bits in the argument is also set in
 220    * the trigger word
 221    */
 222   Bool_t IsTriggerBits(UInt_t bits) const;
 223   /**
 224    * Check if any of bit(s) are enabled in the trigger word.  This is
 225    * an @e or between the selected bits.  If you need and @a and you
 226    * should use the member function IsTriggerBits;
 227    *
 228    * @param bits Bits to check for
 229    *
 230    * @return true if any of the enabled bits in the arguments are also
 231    * enabled in the trigger mask
 232    */
 233   Bool_t IsTriggerOrBits(UInt_t bits) const;
 234   /**
 235    * Whether we have any trigger bits
 236    *
 237    * @return true if we have some trigger
 238    */
 239   Bool_t HasTrigger() const { return fTriggers != 0; } // Check for triggers
 240   /**
 241    * Clear all data
 242    *
 243    * @param option  Passed on to TH2::Reset verbatim
 244    */
 245   void Clear(Option_t* option="");
 246   /**
 247    * browse this object
 248    *
 249    * @param b Browser
 250    */
 251   void Browse(TBrowser* b);
 252   /**
 253    * This is a folder
 254    *
 255    * @return Always true
 256    */
 257   Bool_t IsFolder() const { return kTRUE; } // Always true
 258   /**
 259    * Print content
 260    *
 261    * @param option Passed verbatim to TH2::Print
 262    */
 263   void Print(Option_t* option="") const;
 264   /**
 265    * Set the z coordinate of the interaction point
 266    *
 267    * @param ipZ Interaction point z coordinate
 268    */
 269   void SetIpZ(Float_t ipZ) { fIpZ = ipZ; } // Set Ip's Z coordinate
 270   /**
 271    * Set the center of mass energy per nucleon-pair.  This is stored
 272    * in the (0,0) of the histogram
 273    *
 274    * @param sNN Center of mass energy per nucleon pair (GeV)
 275    */
 276   void SetSNN(UShort_t sNN);
 277   /**
 278    * Get the collision system number
 279    * - 0: Unknown
 280    * - 1: pp
 281    * - 2: PbPb
 282    *
 283    * @param sys Collision system number
 284    */
 285   void SetSystem(UShort_t sys);
 286   /**
 287    * Set the event centrality
 288    *
 289    * @param c Centrality
 290    */
 291   void SetCentrality(Float_t c) { fCentrality = c; }
 292   /**
 293    * Set the z coordinate of the interaction point
 294    *
 295    * @return Interaction point z coordinate
 296    */
 297   Float_t GetIpZ() const { return fIpZ; } // Get Ip's Z coordinate
 298   /**
 299    * Check if we have a valid z coordinate of the interaction point
 300    *
 301    * @return True if we have a valid interaction point z coordinate
 302    */
 303   Bool_t HasIpZ() const;
 304   /**
 305    * Get the center of mass energy per nucleon pair (GeV)
 306    *
 307    * @return Center of mass energy per nucleon pair (GeV)
 308    */
 309   UShort_t GetSNN() const;
 310   /**
 311    * Get the collision system number
 312    * - 0: Unknown
 313    * - 1: pp
 314    * - 2: PbPb
 315    *
 316    * @return Collision system number
 317    */
 318   UShort_t GetSystem() const;
 319   /**
 320    * Check if the z coordinate of the interaction point is within the
 321    * given limits.  Note that the convention used corresponds to the
 322    * convention used in ROOTs TAxis.
 323    *
 324    * @param low  Lower cut (inclusive)
 325    * @param high Upper cut (exclusive)
 326    *
 327    * @return true if @f$ low \ge ipz < high@f$
 328    */
 329   Bool_t InRange(Float_t low, Float_t high) const;
 330   /**
 331    * Get the event centrality
 332    *
 333    *
 334    * @return
 335    */
 336   Float_t GetCentrality() const { return fCentrality; }
 337   /**
 338    * Check if we have a valid centrality
 339    *
 340    *
 341    * @return
 342    */
 343   Bool_t  HasCentrality() const { return !(fCentrality  < 0); }
 344   /**
 345    * Get the number of SPD clusters seen in @f$ |\eta|<1@f$
 346    *
 347    * @return Number of SPD clusters seen
 348    */
 349   UShort_t GetNClusters() const { return fNClusters; }
 350   /**
 351    * Set the number of SPD clusters seen in @f$ |\eta|<1@f$
 352    *
 353    * @param n Number of SPD clusters
 354    */
 355   void SetNClusters(UShort_t n) { fNClusters = n; }
 356   /**
 357    * Get the name of the object
 358    *
 359    * @return Name of object
 360    */
 361   const Char_t* GetName() const { return (fIsMC ? "ForwardMC" : "Forward"); }
 362   /**
 363    * Check if event meets the passses requirements.
 364    *
 365    * It returns true if @e all of the following is true
 366    *
 367    * - The trigger is within the bit mask passed.
 368    * - The vertex is within the specified limits.
 369    * - The centrality is within the specified limits, or if lower
 370    *   limit is equal to or larger than the upper limit.
 371    *
 372    * Note, for data with out a centrality estimate (e.g., pp), one
 373    * must pass equal centrality cuts, or no data will be accepted.  In
 374    * other words, for pp data, always pass cMin=0, cMax=0
 375    *
 376    * If a histogram is passed in the last parameter, then that
 377    * histogram is filled with the trigger bits.
 378    *
 379    * @param triggerMask  Trigger mask
 380    * @param vzMin        Minimum @f$ v_z@f$ (in centimeters)
 381    * @param vzMax        Maximum @f$ v_z@f$ (in centimeters)
 382    * @param cMin         Minimum centrality (in percent)
 383    * @param cMax         Maximum centrality (in percent)
 384    * @param hist         Histogram to fill
 385    *
 386    * @return @c true if the event meets the requirements
 387    */
 388   Bool_t CheckEvent(Int_t    triggerMask=kInel,
 389                     Double_t vzMin=-10, Double_t vzMax=10,
 390                     UShort_t cMin=0,    UShort_t cMax=100,
 391                     TH1*     hist=0) const;
 392   /**
 393    * Get a string correspondig to the trigger mask
 394    *
 395    * @param mask Trigger mask
 396    *
 397    * @return Static string (copy before use)
 398    */
 399   static const Char_t* GetTriggerString(UInt_t mask);
 400   /**
 401    * Make a histogram to record triggers in.
 402    *
 403    * The bins defined by the trigger enumeration in this class.  One
 404    * can use this enumeration to retrieve the number of triggers for
 405    * each class.
 406    *
 407    * @param name Name of the histogram
 408    * @param mask Trigger mask
 409    *
 410    * @return Newly allocated histogram
 411    */
 412   static TH1I* MakeTriggerHistogram(const char* name="triggers",
 413                                     Int_t mask=0);
 414   /**
 415    * Utility function to make a trigger mask from the passed string.
 416    *
 417    * The string is a comma or space seperated list of case-insensitive
 418    * strings
 419    *
 420    * - INEL
 421    * - INEL>0
 422    * - NSD
 423    *
 424    * @param what Which triggers to put in the mask.
 425    *
 426    * @return The generated trigger mask.
 427    */
 428   static UInt_t MakeTriggerMask(const char* what);
 429 protected:
 430   /** From MC or not */
 431   Bool_t   fIsMC;       // Whether this is from MC
 432   /** Histogram of @f$d^2N_{ch}/(d\eta d\phi)@f$ for this event */
 433   TH2D     fHist;       // Histogram of d^2N_{ch}/(deta dphi) for this event
 434   /** Trigger bits */
 435   UInt_t   fTriggers;   // Trigger bit mask
 436   /** Interaction point @f$z@f$ coordinate */
 437   Float_t  fIpZ;        // Z coordinate of the interaction point
 438   /** Centrality */
 439   Float_t  fCentrality; // Event centrality
 440   /** Number of clusters in @f$|\eta|<1@f$ */
 441   UShort_t fNClusters;  // Number of SPD clusters in |eta|<1
 442   /** Invalid value for interaction point @f$z@f$ coordiante */
 443   static const Float_t fgkInvalidIpZ; // Invalid IpZ value
 444   ClassDef(AliAODForwardMult,3); // AOD forward multiplicity
 445 };
 446
 447 //____________________________________________________________________
 448 inline Bool_t
 449 AliAODForwardMult::InRange(Float_t low, Float_t high) const
 450 {
 451   return HasIpZ() && fIpZ >= low && fIpZ < high;
 452 }
 453
 454 //____________________________________________________________________
 455 inline Bool_t
 456 AliAODForwardMult::IsTriggerBits(UInt_t bits) const
 457 {
 458   return HasTrigger() && ((fTriggers & bits) == bits);
 459 }
 460 //____________________________________________________________________
 461 inline Bool_t
 462 AliAODForwardMult::IsTriggerOrBits(UInt_t bits) const
 463 {
 464   return HasTrigger() && ((fTriggers & bits) != 0);
 465 }
 466
 467 #endif
 468 // Local Variables:
 469 //  mode: C++
 470 // End:
 471