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