2c03bdf9233cfff1d8b950965bb910a6b55a9c7f
[u/mrichter/AliRoot.git] / 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 "AliFlowVector.h"
25 #include "AliFlowTrackSimple.h"
26 #include "AliFlowEventSimple.h"
27
28 /**************************************
29  * AliFlowEventSimple: A simple event *
30  *  for flow analysis                 * 
31  *                                    * 
32  * authors: Naomi van der Kolk        *
33  *           (kolk@nikhef.nl)         *  
34  *          Ante Bilandzic            *
35  *           (anteb@nikhef.nl)        * 
36  * ***********************************/
37  
38 ClassImp(AliFlowEventSimple)
39
40 //-----------------------------------------------------------------------
41
42   AliFlowEventSimple::AliFlowEventSimple(Int_t aLenght):
43     fTrackCollection(NULL),
44     fNumberOfTracks(0),
45     fEventNSelTracksIntFlow(0)
46 {
47   //constructor 
48   fTrackCollection =  new TObjArray(aLenght) ;
49 }
50
51 //-----------------------------------------------------------------------
52
53 AliFlowEventSimple::AliFlowEventSimple(const AliFlowEventSimple& anEvent):
54   TObject(),
55   fTrackCollection(anEvent.fTrackCollection),
56   fNumberOfTracks(anEvent.fNumberOfTracks),
57   fEventNSelTracksIntFlow(anEvent.fEventNSelTracksIntFlow)
58 {
59   //copy constructor 
60 }
61
62 //-----------------------------------------------------------------------
63
64 AliFlowEventSimple& AliFlowEventSimple::operator=(const AliFlowEventSimple& anEvent)
65 {
66   *fTrackCollection = *anEvent.fTrackCollection ;
67   fNumberOfTracks = anEvent.fNumberOfTracks;
68   fEventNSelTracksIntFlow = anEvent.fEventNSelTracksIntFlow;
69   
70   return *this;
71 }
72
73 //----------------------------------------------------------------------- 
74
75 AliFlowEventSimple::~AliFlowEventSimple()
76 {
77   //destructor
78   fTrackCollection->Delete() ; delete fTrackCollection ;
79 }
80
81 //----------------------------------------------------------------------- 
82
83 AliFlowTrackSimple* AliFlowEventSimple::GetTrack(Int_t i)
84 {
85   //get track i from collection
86   AliFlowTrackSimple* pTrack = (AliFlowTrackSimple*)TrackCollection()->At(i) ;
87   return pTrack;
88 }
89
90 //-----------------------------------------------------------------------   
91 AliFlowVector AliFlowEventSimple::GetQ(Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights) 
92 {
93   // calculate Q-vector in harmonic n without weights (default harmonic n=2)  
94   Double_t dQX = 0.;
95   Double_t dQY = 0.;
96   AliFlowVector vQ;
97   vQ.Set(0.,0.);
98   
99   Int_t iOrder = n;
100   Double_t iUsedTracks = 0;
101   Double_t dPhi=0.;
102   Double_t dPt=0.;
103   Double_t dEta=0.;
104   
105   AliFlowTrackSimple* pTrack = NULL;
106  
107   Int_t nBinsPhi=0; 
108   Double_t dBinWidthPt=0.;
109   Double_t dPtMin=0.;
110   Double_t dBinWidthEta=0.;
111   Double_t dEtaMin=0.;
112  
113   Double_t wPhi=1.; // weight Phi  
114   Double_t wPt=1.;  // weight Pt 
115   Double_t wEta=1.; // weight Eta 
116   
117   TH1F *phiWeights = NULL;
118   TH1D *ptWeights  = NULL;
119   TH1D *etaWeights = NULL;
120
121   Double_t dSumOfWeightsToPower2 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 2)
122   Double_t dSumOfWeightsToPower3 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 3)
123   Double_t dSumOfWeightsToPower4 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 4)
124   Double_t dSumOfWeightsToPower5 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 5)
125   Double_t dSumOfWeightsToPower6 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 6)
126   Double_t dSumOfWeightsToPower7 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 7)
127   Double_t dSumOfWeightsToPower8 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 8) 
128
129   if(weightsList)
130   {
131    if(usePhiWeights)
132    {
133     phiWeights = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights"));
134     if(phiWeights) nBinsPhi = phiWeights->GetNbinsX();
135    }          
136    if(usePtWeights)
137    {
138     ptWeights = dynamic_cast<TH1D *>(weightsList->FindObject("pt_weights"));
139     if(ptWeights)
140     {
141      dBinWidthPt = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
142      dPtMin = (ptWeights->GetXaxis())->GetXmin();
143     } 
144    }       
145    if(useEtaWeights)
146    {
147     etaWeights = dynamic_cast<TH1D *>(weightsList->FindObject("eta_weights"));
148     if(etaWeights)
149     {
150      dBinWidthEta = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
151      dEtaMin = (etaWeights->GetXaxis())->GetXmin();
152     } 
153    }          
154   } // end of if(weightsList)
155   
156   // loop over tracks    
157   for(Int_t i=0;i<fNumberOfTracks;i++)                               
158   {
159    pTrack = (AliFlowTrackSimple*)TrackCollection()->At(i); 
160    if(pTrack)
161    {
162     if(pTrack->UseForIntegratedFlow()) 
163     {
164      dPhi = pTrack->Phi();
165      dPt  = pTrack->Pt();
166      dEta = pTrack->Eta();
167      
168      // determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
169      if(phiWeights && nBinsPhi)
170      {
171       wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi())));
172      }
173      // determine v'(pt) weight:    
174      if(ptWeights && dBinWidthPt)
175      {
176       wPt=ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-dPtMin)/dBinWidthPt))); 
177      }            
178      // determine v'(eta) weight:    
179      if(etaWeights && dBinWidthEta)
180      {
181       wEta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-dEtaMin)/dBinWidthEta))); 
182      } 
183
184      // building up the weighted Q-vector:       
185      dQX += wPhi*wPt*wEta*TMath::Cos(iOrder*dPhi);
186      dQY += wPhi*wPt*wEta*TMath::Sin(iOrder*dPhi); 
187     
188      // weighted multiplicity:
189      iUsedTracks+=wPhi*wPt*wEta;
190     
191      // weights raised to various powers are summed up:
192      dSumOfWeightsToPower2+=pow(wPhi*wPt*wEta, 2); 
193      dSumOfWeightsToPower3+=pow(wPhi*wPt*wEta, 3); 
194      dSumOfWeightsToPower4+=pow(wPhi*wPt*wEta, 4); 
195      dSumOfWeightsToPower5+=pow(wPhi*wPt*wEta, 5); 
196      dSumOfWeightsToPower6+=pow(wPhi*wPt*wEta, 6); 
197      dSumOfWeightsToPower7+=pow(wPhi*wPt*wEta, 7); 
198      dSumOfWeightsToPower8+=pow(wPhi*wPt*wEta, 8); 
199      
200     } // end of if (pTrack->UseForIntegratedFlow())
201    } // end of if (pTrack)
202    else {cerr << "no particle!!!"<<endl;}
203   } // loop over particles
204     
205   vQ.Set(dQX,dQY);
206   vQ.SetMult(iUsedTracks);
207   vQ.SetSumOfWeightsToPower2(dSumOfWeightsToPower2);
208   vQ.SetSumOfWeightsToPower3(dSumOfWeightsToPower3);
209   vQ.SetSumOfWeightsToPower4(dSumOfWeightsToPower4);
210   vQ.SetSumOfWeightsToPower5(dSumOfWeightsToPower5);
211   vQ.SetSumOfWeightsToPower6(dSumOfWeightsToPower6);
212   vQ.SetSumOfWeightsToPower7(dSumOfWeightsToPower7);
213   vQ.SetSumOfWeightsToPower8(dSumOfWeightsToPower8);
214
215   return vQ;
216   
217 }
218
219
220
221
222
223
224
225
226
227