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