PWG2/FLOW/AliFlowCommon/AliFlowEventSimple.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 #include "Riostream.h"
  17 #include "TObjArray.h"
  18 #include "TFile.h"
  19 #include "TList.h"
  20 #include "TMath.h"
  21 #include "TH1F.h"
  22 #include "TH1D.h"
  23 #include "TProfile.h"
  24 #include "TBrowser.h"
  25 #include "AliFlowVector.h"
  26 #include "AliFlowTrackSimple.h"
  27 #include "AliFlowEventSimple.h"
  28
  29 /**************************************
  30  * AliFlowEventSimple: A simple event *
  31  *  for flow analysis                 *
  32  *                                    *
  33  * authors: Naomi van der Kolk        *
  34  *           (kolk@nikhef.nl)         *
  35  *          Ante Bilandzic            *
  36  *           (anteb@nikhef.nl)        *
  37  * ***********************************/
  38
  39 ClassImp(AliFlowEventSimple)
  40
  41 //-----------------------------------------------------------------------
  42
  43 AliFlowEventSimple::AliFlowEventSimple():
  44   fTrackCollection(NULL),
  45   fNumberOfTracks(0),
  46   fEventNSelTracksIntFlow(0),
  47   fNumberOfTracksWrap(NULL),
  48   fEventNSelTracksIntFlowWrap(NULL),
  49   fMCReactionPlaneAngleWrap(NULL)
  50 {
  51   cout << "AliFlowEventSimple: Default constructor to be used only by root for io" << endl;
  52 }
  53
  54 //-----------------------------------------------------------------------
  55
  56 AliFlowEventSimple::AliFlowEventSimple(Int_t aLenght):
  57   fTrackCollection(NULL),
  58   fNumberOfTracks(0),
  59   fEventNSelTracksIntFlow(0),
  60   fNumberOfTracksWrap(NULL),
  61   fEventNSelTracksIntFlowWrap(NULL),
  62   fMCReactionPlaneAngleWrap(NULL)
  63 {
  64   //constructor
  65   fTrackCollection =  new TObjArray(aLenght) ;
  66 }
  67
  68 //-----------------------------------------------------------------------
  69
  70 AliFlowEventSimple::AliFlowEventSimple(const AliFlowEventSimple& anEvent):
  71   TObject(),
  72   fTrackCollection(anEvent.fTrackCollection),
  73   fNumberOfTracks(anEvent.fNumberOfTracks),
  74   fEventNSelTracksIntFlow(anEvent.fEventNSelTracksIntFlow),
  75   fMCReactionPlaneAngle(anEvent.fMCReactionPlaneAngle),
  76   fNumberOfTracksWrap(anEvent.fNumberOfTracksWrap),
  77   fEventNSelTracksIntFlowWrap(anEvent.fEventNSelTracksIntFlowWrap),
  78   fMCReactionPlaneAngleWrap(anEvent.fMCReactionPlaneAngleWrap)
  79 {
  80   //copy constructor
  81 }
  82
  83 //-----------------------------------------------------------------------
  84
  85 AliFlowEventSimple& AliFlowEventSimple::operator=(const AliFlowEventSimple& anEvent)
  86 {
  87   *fTrackCollection = *anEvent.fTrackCollection ;
  88   fNumberOfTracks = anEvent.fNumberOfTracks;
  89   fEventNSelTracksIntFlow = anEvent.fEventNSelTracksIntFlow;
  90   fMCReactionPlaneAngle = anEvent.fMCReactionPlaneAngle;
  91   fNumberOfTracksWrap = anEvent.fNumberOfTracksWrap;
  92   fEventNSelTracksIntFlowWrap = anEvent.fEventNSelTracksIntFlowWrap;
  93   fMCReactionPlaneAngleWrap=anEvent.fMCReactionPlaneAngleWrap;
  94
  95   return *this;
  96 }
  97
  98 //-----------------------------------------------------------------------
  99
 100 AliFlowEventSimple::~AliFlowEventSimple()
 101 {
 102   //destructor
 103   if (fTrackCollection) fTrackCollection->Delete(); delete fTrackCollection;
 104   if (fNumberOfTracksWrap) delete fNumberOfTracksWrap;
 105   if (fEventNSelTracksIntFlowWrap) delete fEventNSelTracksIntFlowWrap;
 106   if (fMCReactionPlaneAngleWrap) delete fMCReactionPlaneAngleWrap;
 107 }
 108
 109 //-----------------------------------------------------------------------
 110
 111 AliFlowTrackSimple* AliFlowEventSimple::GetTrack(Int_t i)
 112 {
 113   //get track i from collection
 114   AliFlowTrackSimple* pTrack = (AliFlowTrackSimple*)TrackCollection()->At(i) ;
 115   return pTrack;
 116 }
 117
 118 //-----------------------------------------------------------------------
 119 AliFlowVector AliFlowEventSimple::GetQ(Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights)
 120 {
 121   // calculate Q-vector in harmonic n without weights (default harmonic n=2)
 122   Double_t dQX = 0.;
 123   Double_t dQY = 0.;
 124   AliFlowVector vQ;
 125   vQ.Set(0.,0.);
 126
 127   Int_t iOrder = n;
 128   Double_t iUsedTracks = 0;
 129   Double_t dPhi=0.;
 130   Double_t dPt=0.;
 131   Double_t dEta=0.;
 132
 133   AliFlowTrackSimple* pTrack = NULL;
 134
 135   Int_t nBinsPhi=0;
 136   Double_t dBinWidthPt=0.;
 137   Double_t dPtMin=0.;
 138   Double_t dBinWidthEta=0.;
 139   Double_t dEtaMin=0.;
 140
 141   Double_t wPhi=1.; // weight Phi
 142   Double_t wPt=1.;  // weight Pt
 143   Double_t wEta=1.; // weight Eta
 144
 145   TH1F *phiWeights = NULL;
 146   TH1D *ptWeights  = NULL;
 147   TH1D *etaWeights = NULL;
 148
 149   Double_t dSumOfWeightsToPower2 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 2)
 150   Double_t dSumOfWeightsToPower3 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 3)
 151   Double_t dSumOfWeightsToPower4 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 4)
 152   Double_t dSumOfWeightsToPower5 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 5)
 153   Double_t dSumOfWeightsToPower6 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 6)
 154   Double_t dSumOfWeightsToPower7 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 7)
 155   Double_t dSumOfWeightsToPower8 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 8)
 156
 157   if(weightsList)
 158   {
 159    if(usePhiWeights)
 160    {
 161     phiWeights = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights"));
 162     if(phiWeights) nBinsPhi = phiWeights->GetNbinsX();
 163    }
 164    if(usePtWeights)
 165    {
 166     ptWeights = dynamic_cast<TH1D *>(weightsList->FindObject("pt_weights"));
 167     if(ptWeights)
 168     {
 169      dBinWidthPt = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
 170      dPtMin = (ptWeights->GetXaxis())->GetXmin();
 171     }
 172    }
 173    if(useEtaWeights)
 174    {
 175     etaWeights = dynamic_cast<TH1D *>(weightsList->FindObject("eta_weights"));
 176     if(etaWeights)
 177     {
 178      dBinWidthEta = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
 179      dEtaMin = (etaWeights->GetXaxis())->GetXmin();
 180     }
 181    }
 182   } // end of if(weightsList)
 183
 184   // loop over tracks
 185   for(Int_t i=0;i<fNumberOfTracks;i++)
 186   {
 187    pTrack = (AliFlowTrackSimple*)TrackCollection()->At(i);
 188    if(pTrack)
 189    {
 190     if(pTrack->UseForIntegratedFlow())
 191     {
 192      dPhi = pTrack->Phi();
 193      dPt  = pTrack->Pt();
 194      dEta = pTrack->Eta();
 195
 196      // determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
 197      if(phiWeights && nBinsPhi)
 198      {
 199       wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi())));
 200      }
 201      // determine v'(pt) weight:
 202      if(ptWeights && dBinWidthPt)
 203      {
 204       wPt=ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-dPtMin)/dBinWidthPt)));
 205      }
 206      // determine v'(eta) weight:
 207      if(etaWeights && dBinWidthEta)
 208      {
 209       wEta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-dEtaMin)/dBinWidthEta)));
 210      }
 211
 212      // building up the weighted Q-vector:
 213      dQX += wPhi*wPt*wEta*TMath::Cos(iOrder*dPhi);
 214      dQY += wPhi*wPt*wEta*TMath::Sin(iOrder*dPhi);
 215
 216      // weighted multiplicity:
 217      iUsedTracks+=wPhi*wPt*wEta;
 218
 219      // weights raised to various powers are summed up:
 220      dSumOfWeightsToPower2+=pow(wPhi*wPt*wEta, 2);
 221      dSumOfWeightsToPower3+=pow(wPhi*wPt*wEta, 3);
 222      dSumOfWeightsToPower4+=pow(wPhi*wPt*wEta, 4);
 223      dSumOfWeightsToPower5+=pow(wPhi*wPt*wEta, 5);
 224      dSumOfWeightsToPower6+=pow(wPhi*wPt*wEta, 6);
 225      dSumOfWeightsToPower7+=pow(wPhi*wPt*wEta, 7);
 226      dSumOfWeightsToPower8+=pow(wPhi*wPt*wEta, 8);
 227
 228     } // end of if (pTrack->UseForIntegratedFlow())
 229    } // end of if (pTrack)
 230    else {cerr << "no particle!!!"<<endl;}
 231   } // loop over particles
 232
 233   vQ.Set(dQX,dQY);
 234   vQ.SetMult(iUsedTracks);
 235   vQ.SetSumOfWeightsToPower2(dSumOfWeightsToPower2);
 236   vQ.SetSumOfWeightsToPower3(dSumOfWeightsToPower3);
 237   vQ.SetSumOfWeightsToPower4(dSumOfWeightsToPower4);
 238   vQ.SetSumOfWeightsToPower5(dSumOfWeightsToPower5);
 239   vQ.SetSumOfWeightsToPower6(dSumOfWeightsToPower6);
 240   vQ.SetSumOfWeightsToPower7(dSumOfWeightsToPower7);
 241   vQ.SetSumOfWeightsToPower8(dSumOfWeightsToPower8);
 242
 243   return vQ;
 244
 245 }
 246
 247 //-----------------------------------------------------------------------
 248 void AliFlowEventSimple::Print(Option_t *option) const
 249 {
 250   //   -*-*-*-*-*Print some global quantities for this histogram collection class *-*-*-*-*-*-*-*
 251   //             ===============================================
 252   //   printf( "TH1.Print Name  = %s, Entries= %d, Total sum= %g\n",GetName(),Int_t(fEntries),GetSumOfWeights());
 253   printf( "Class.Print Name = %s, Total number of tracks= %d, Number of selected tracks= %d, MC EventPlaneAngle= %d",
 254           GetName(),fNumberOfTracks, fEventNSelTracksIntFlow, fMCReactionPlaneAngle );
 255
 256   if (fTrackCollection) {
 257     fTrackCollection->Print(option);
 258   }
 259   else {
 260     printf( "Empty track collection \n");
 261   }
 262 }
 263
 264 //-----------------------------------------------------------------------
 265  void AliFlowEventSimple::Browse(TBrowser *b)
 266 {
 267   if (!b) return;
 268   if (!fNumberOfTracksWrap) {
 269     fNumberOfTracksWrap = new TParameter<int>("fNumberOfTracks", fNumberOfTracks);
 270     b->Add(fNumberOfTracksWrap);
 271   }
 272   if (!fEventNSelTracksIntFlowWrap) {
 273     fEventNSelTracksIntFlowWrap = new TParameter<int>("fEventNSelTracksIntFlow", fEventNSelTracksIntFlow);
 274     b->Add(fEventNSelTracksIntFlowWrap);
 275   }
 276   if (!fMCReactionPlaneAngleWrap) {
 277      fMCReactionPlaneAngleWrap = new TParameter<double>(" fMCReactionPlaneAngle",  fMCReactionPlaneAngle);
 278     b->Add( fMCReactionPlaneAngleWrap);
 279   }
 280   if (fTrackCollection) b->Add(fTrackCollection,"AliFlowTracksSimple");
 281 }
 282
 283
 284
 285
 286
 287
 288
 289
 290
 291
 292