]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/FLOW/Base/AliFlowEventSimple.cxx
generatio of various plaes in on the fly
[u/mrichter/AliRoot.git] / PWG / FLOW / Base / AliFlowEventSimple.cxx
CommitLineData
f1d945a1 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
929098e4 16/*****************************************************************
17 AliFlowEventSimple: A simple event
18 for flow analysis
19
20 origin: Naomi van der Kolk (kolk@nikhef.nl)
21 Ante Bilandzic (anteb@nikhef.nl)
22 Raimond Snellings (Raimond.Snellings@nikhef.nl)
23 mods: Mikolaj Krzewicki (mikolaj.krzewicki@cern.ch)
24*****************************************************************/
25
f1d945a1 26#include "Riostream.h"
27#include "TObjArray.h"
26c4cbb9 28#include "TFile.h"
03a02aca 29#include "TList.h"
929098e4 30#include "TTree.h"
31#include "TParticle.h"
f1d945a1 32#include "TMath.h"
26c4cbb9 33#include "TH1F.h"
34#include "TH1D.h"
bc231a12 35#include "TF1.h"
26c4cbb9 36#include "TProfile.h"
929098e4 37#include "TParameter.h"
c076fda8 38#include "TBrowser.h"
f1d945a1 39#include "AliFlowVector.h"
40#include "AliFlowTrackSimple.h"
929098e4 41#include "AliFlowTrackSimpleCuts.h"
701f71c1 42#include "AliFlowEventSimple.h"
7382279b 43#include "TRandom.h"
f1d945a1 44
3a7af7bd 45using std::cout;
46using std::endl;
47using std::cerr;
f1d945a1 48ClassImp(AliFlowEventSimple)
49
50//-----------------------------------------------------------------------
46bec39c 51AliFlowEventSimple::AliFlowEventSimple():
52 fTrackCollection(NULL),
85d2ee8d 53 fReferenceMultiplicity(0),
46bec39c 54 fNumberOfTracks(0),
bc231a12 55 fNumberOfRPs(0),
df23c5ae 56 fNumberOfPOIs(0),
8d1c89ad 57 fUseGlauberMCSymmetryPlanes(kFALSE),
58 fUseExternalSymmetryPlanes(kFALSE),
59 fPsi1(0.),
60 fPsi2(0.),
61 fPsi3(0.),
62 fPsi4(0.),
63 fPsi5(0.),
64 fPsi1Psi3(0x0),
65 fPsi2Psi4(0x0),
66 fPsi3Psi5(0x0),
7183fe85 67 fMCReactionPlaneAngle(0.),
929098e4 68 fMCReactionPlaneAngleIsSet(kFALSE),
bc231a12 69 fAfterBurnerPrecision(0.001),
9fa64edc 70 fUserModified(kFALSE),
c076fda8 71 fNumberOfTracksWrap(NULL),
bc231a12 72 fNumberOfRPsWrap(NULL),
df23c5ae 73 fNumberOfPOIsWrap(NULL),
ad4d7e4e 74 fMCReactionPlaneAngleWrap(NULL),
75 fShuffledIndexes(NULL),
76 fShuffleTracks(kFALSE)
46bec39c 77{
78 cout << "AliFlowEventSimple: Default constructor to be used only by root for io" << endl;
79}
80
81//-----------------------------------------------------------------------
9fa64edc 82AliFlowEventSimple::AliFlowEventSimple( Int_t n,
83 ConstructionMethod method,
bc231a12 84 TF1* ptDist,
85 Double_t phiMin,
86 Double_t phiMax,
87 Double_t etaMin,
88 Double_t etaMax):
9fa64edc 89 fTrackCollection(new TObjArray(n)),
a1c43d26 90 fReferenceMultiplicity(0),
bc231a12 91 fNumberOfTracks(0),
92 fNumberOfRPs(0),
df23c5ae 93 fNumberOfPOIs(0),
8d1c89ad 94 fUseGlauberMCSymmetryPlanes(kFALSE),
95 fUseExternalSymmetryPlanes(kFALSE),
96 fPsi1(0.),
97 fPsi2(0.),
98 fPsi3(0.),
99 fPsi4(0.),
100 fPsi5(0.),
101 fPsi1Psi3(0x0),
102 fPsi2Psi4(0x0),
103 fPsi3Psi5(0x0),
bc231a12 104 fMCReactionPlaneAngle(0.),
105 fMCReactionPlaneAngleIsSet(kFALSE),
106 fAfterBurnerPrecision(0.001),
9fa64edc 107 fUserModified(kFALSE),
bc231a12 108 fNumberOfTracksWrap(NULL),
109 fNumberOfRPsWrap(NULL),
df23c5ae 110 fNumberOfPOIsWrap(NULL),
ad4d7e4e 111 fMCReactionPlaneAngleWrap(NULL),
112 fShuffledIndexes(NULL),
113 fShuffleTracks(kFALSE)
bc231a12 114{
9fa64edc 115 //ctor
116 // if second argument is set to AliFlowEventSimple::kGenerate
117 // it generates n random tracks with given Pt distribution
118 // (a sane default is provided), phi and eta are uniform
119
120 if (method==kGenerate)
121 Generate(n,ptDist,phiMin,phiMax,etaMin,etaMax);
f1d945a1 122}
123
124//-----------------------------------------------------------------------
e35ddff0 125AliFlowEventSimple::AliFlowEventSimple(const AliFlowEventSimple& anEvent):
bc231a12 126 TObject(anEvent),
127 fTrackCollection((TObjArray*)(anEvent.fTrackCollection)->Clone()),
85d2ee8d 128 fReferenceMultiplicity(anEvent.fReferenceMultiplicity),
e35ddff0 129 fNumberOfTracks(anEvent.fNumberOfTracks),
bc231a12 130 fNumberOfRPs(anEvent.fNumberOfRPs),
df23c5ae 131 fNumberOfPOIs(anEvent.fNumberOfPOIs),
8d1c89ad 132 fUseGlauberMCSymmetryPlanes(anEvent.fUseGlauberMCSymmetryPlanes),
133 fUseExternalSymmetryPlanes(anEvent.fUseExternalSymmetryPlanes),
134 fPsi1(anEvent.fPsi1),
135 fPsi2(anEvent.fPsi2),
136 fPsi3(anEvent.fPsi3),
137 fPsi4(anEvent.fPsi4),
138 fPsi5(anEvent.fPsi5),
139 fPsi1Psi3(anEvent.fPsi1Psi3),
140 fPsi2Psi4(anEvent.fPsi2Psi4),
141 fPsi3Psi5(anEvent.fPsi3Psi5),
a12990bb 142 fMCReactionPlaneAngle(anEvent.fMCReactionPlaneAngle),
929098e4 143 fMCReactionPlaneAngleIsSet(anEvent.fMCReactionPlaneAngleIsSet),
bc231a12 144 fAfterBurnerPrecision(anEvent.fAfterBurnerPrecision),
9fa64edc 145 fUserModified(anEvent.fUserModified),
c076fda8 146 fNumberOfTracksWrap(anEvent.fNumberOfTracksWrap),
bc231a12 147 fNumberOfRPsWrap(anEvent.fNumberOfRPsWrap),
df23c5ae 148 fNumberOfPOIsWrap(anEvent.fNumberOfPOIsWrap),
ad4d7e4e 149 fMCReactionPlaneAngleWrap(anEvent.fMCReactionPlaneAngleWrap),
150 fShuffledIndexes(NULL),
151 fShuffleTracks(anEvent.fShuffleTracks)
f1d945a1 152{
929098e4 153 //copy constructor
f1d945a1 154}
155
156//-----------------------------------------------------------------------
e35ddff0 157AliFlowEventSimple& AliFlowEventSimple::operator=(const AliFlowEventSimple& anEvent)
f1d945a1 158{
124fb262 159 //assignment operator
2047680a 160 if (&anEvent==this) return *this; //check self-assignment
701f71c1 161 if (fTrackCollection) fTrackCollection->Delete();
124fb262 162 delete fTrackCollection;
bc231a12 163 fTrackCollection = (TObjArray*)(anEvent.fTrackCollection)->Clone(); //deep copy
85d2ee8d 164 fReferenceMultiplicity = anEvent.fReferenceMultiplicity;
e35ddff0 165 fNumberOfTracks = anEvent.fNumberOfTracks;
bc231a12 166 fNumberOfRPs = anEvent.fNumberOfRPs;
df23c5ae 167 fNumberOfPOIs = anEvent.fNumberOfPOIs;
8d1c89ad 168 fUseGlauberMCSymmetryPlanes = anEvent.fUseGlauberMCSymmetryPlanes;
169 fUseExternalSymmetryPlanes = anEvent.fUseExternalSymmetryPlanes;
170 fPsi1 = anEvent.fPsi1;
171 fPsi2 = anEvent.fPsi2;
172 fPsi3 = anEvent.fPsi3;
173 fPsi4 = anEvent.fPsi4;
174 fPsi5 = anEvent.fPsi5;
175 fPsi1Psi3 = anEvent.fPsi1Psi3;
176 fPsi2Psi4 = anEvent.fPsi2Psi4;
177 fPsi3Psi5 = anEvent.fPsi3Psi5;
a12990bb 178 fMCReactionPlaneAngle = anEvent.fMCReactionPlaneAngle;
929098e4 179 fMCReactionPlaneAngleIsSet = anEvent.fMCReactionPlaneAngleIsSet;
bc231a12 180 fAfterBurnerPrecision = anEvent.fAfterBurnerPrecision;
9fa64edc 181 fUserModified=anEvent.fUserModified;
929098e4 182 fNumberOfTracksWrap = anEvent.fNumberOfTracksWrap;
bc231a12 183 fNumberOfRPsWrap = anEvent.fNumberOfRPsWrap;
df23c5ae 184 fNumberOfPOIsWrap = anEvent.fNumberOfPOIsWrap;
a12990bb 185 fMCReactionPlaneAngleWrap=anEvent.fMCReactionPlaneAngleWrap;
ad4d7e4e 186 fShuffleTracks=anEvent.fShuffleTracks;
187 delete [] fShuffledIndexes;
f1d945a1 188 return *this;
f1d945a1 189}
190
929098e4 191//-----------------------------------------------------------------------
f1d945a1 192AliFlowEventSimple::~AliFlowEventSimple()
193{
194 //destructor
929098e4 195 if (fTrackCollection) fTrackCollection->Delete();
196 delete fTrackCollection;
7d27a354 197 delete fNumberOfTracksWrap;
198 delete fNumberOfRPsWrap;
df23c5ae 199 delete fNumberOfPOIsWrap;
7d27a354 200 delete fMCReactionPlaneAngleWrap;
ad4d7e4e 201 delete fShuffledIndexes;
f1d945a1 202}
203
8d1c89ad 204//-----------------------------------------------------------------------
205void AliFlowEventSimple::SetUseExternalSymmetryPlanes(TF1 *gPsi1Psi3,
206 TF1 *gPsi2Psi4,
207 TF1 *gPsi3Psi5) {
208 //Use symmetry planes, setup correlations between different Psi_n
209 fUseExternalSymmetryPlanes = kTRUE;
210
211 //Correlations between Psi_1 and Psi_3
212 if(gPsi1Psi3) fPsi1Psi3 = gPsi1Psi3;
213 else {
214 fPsi1Psi3 = new TF1("fPsi1Psi3","[0]*x+[1]",0.,2.*TMath::Pi());
215 fPsi1Psi3->SetParameter(0,1.);
216 fPsi1Psi3->SetParameter(1,0.);
217 }
218
219 //Correlations between Psi_2 and Psi_4
220 if(gPsi2Psi4) fPsi2Psi4 = gPsi2Psi4;
221 else {
222 fPsi2Psi4 = new TF1("fPsi2Psi4","[0]*x+[1]",0.,2.*TMath::Pi());
223 fPsi2Psi4->SetParameter(0,1.);
224 fPsi2Psi4->SetParameter(1,0.);
225 }
226
227 //Correlations between Psi_3 and Psi_5
228 if(gPsi3Psi5) fPsi3Psi5 = gPsi3Psi5;
229 else {
230 fPsi3Psi5 = new TF1("fPsi3Psi5","[0]*x+[1]",0.,2.*TMath::Pi());
231 fPsi3Psi5->SetParameter(0,1.);
232 fPsi3Psi5->SetParameter(1,0.);
233 }
234}
235
bc231a12 236//-----------------------------------------------------------------------
237void AliFlowEventSimple::Generate(Int_t nParticles,
6e027e29 238 TF1* ptDist,
239 Double_t phiMin,
240 Double_t phiMax,
241 Double_t etaMin,
242 Double_t etaMax)
bc231a12 243{
244 //generate nParticles random tracks uniform in phi and eta
245 //according to the specified pt distribution
6e027e29 246 if (!ptDist)
247 {
248 static TF1 ptdistribution("ptSpectra","x*TMath::Exp(-pow(0.13957*0.13957+x*x,0.5)/0.4)",0.1,10.);
249 ptDist=&ptdistribution;
250 }
9fa64edc 251
bc231a12 252 for (Int_t i=0; i<nParticles; i++)
253 {
701f71c1 254 AliFlowTrackSimple* track = new AliFlowTrackSimple();
255 track->SetPhi( gRandom->Uniform(phiMin,phiMax) );
256 track->SetEta( gRandom->Uniform(etaMin,etaMax) );
257 track->SetPt( ptDist->GetRandom() );
258 track->SetCharge( (gRandom->Uniform()-0.5<0)?-1:1 );
259 AddTrack(track);
bc231a12 260 }
8d1c89ad 261 if(fUseExternalSymmetryPlanes) {
262 Double_t betaParameter = gRandom->Gaus(0.,1.3);
263 fPsi1Psi3->SetParameter(1,betaParameter);
264
265 betaParameter = gRandom->Gaus(0.,0.9);
266 fPsi2Psi4->SetParameter(1,betaParameter);
267
268 betaParameter = gRandom->Gaus(0.,1.5);
269 fPsi3Psi5->SetParameter(1,betaParameter);
270
271 fPsi1 = gRandom->Uniform(2.*TMath::Pi());
272 fPsi2 = gRandom->Uniform(2.*TMath::Pi());
273 fPsi3 = fPsi1Psi3->Eval(fPsi1);
274 if(fPsi3 < 0) fPsi3 += 2.*TMath::Pi();
275 else if(fPsi3 > 2.*TMath::Pi()) fPsi3 -= 2.*TMath::Pi();
276 fPsi4 = fPsi2Psi4->Eval(fPsi2);
277 if(fPsi4 < 0) fPsi4 += 2.*TMath::Pi();
278 else if(fPsi4 > 2.*TMath::Pi()) fPsi4 -= 2.*TMath::Pi();
279 fPsi5 = fPsi3Psi5->Eval(fPsi3);
280 if(fPsi5 < 0) fPsi5 += 2.*TMath::Pi();
281 else if(fPsi5 > 2.*TMath::Pi()) fPsi5 -= 2.*TMath::Pi();
282
283 fMCReactionPlaneAngle=fPsi2;
284 fMCReactionPlaneAngleIsSet=kTRUE;
285 }
286 else {
287 fMCReactionPlaneAngle=gRandom->Uniform(0.0,TMath::TwoPi());
288 fMCReactionPlaneAngleIsSet=kTRUE;
289 }
9fa64edc 290 SetUserModified();
bc231a12 291}
292
929098e4 293//-----------------------------------------------------------------------
f1d945a1 294AliFlowTrackSimple* AliFlowEventSimple::GetTrack(Int_t i)
295{
296 //get track i from collection
7382279b 297 if (i>=fNumberOfTracks) return NULL;
ad4d7e4e 298 Int_t trackIndex=i;
299 //if asked use the shuffled index
300 if (fShuffleTracks)
301 {
ba4dd29f 302 if (!fShuffledIndexes) ShuffleTracks();
ad4d7e4e 303 trackIndex=fShuffledIndexes[i];
304 }
305 AliFlowTrackSimple* pTrack = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(trackIndex)) ;
e35ddff0 306 return pTrack;
f1d945a1 307}
308
ba4dd29f 309//-----------------------------------------------------------------------
310void AliFlowEventSimple::ShuffleTracks()
311{
312 //shuffle track indexes
313 if (!fShuffledIndexes)
314 {
315 //initialize the table with shuffeled indexes
316 fShuffledIndexes = new Int_t[fNumberOfTracks];
317 for (Int_t j=0; j<fNumberOfTracks; j++) { fShuffledIndexes[j]=j; }
318 }
319 //shuffle
320 std::random_shuffle(&fShuffledIndexes[0], &fShuffledIndexes[fNumberOfTracks]);
321}
322
7382279b 323//-----------------------------------------------------------------------
324void AliFlowEventSimple::AddTrack( AliFlowTrackSimple* track )
325{
7d27a354 326 //add a track, delete the old one if necessary
327 if (fNumberOfTracks < fTrackCollection->GetEntriesFast())
328 {
329 TObject* o = fTrackCollection->At(fNumberOfTracks);
330 if (o) delete o;
331 }
332 fTrackCollection->AddAtAndExpand(track,fNumberOfTracks);
7382279b 333 fNumberOfTracks++;
37b19b98 334 delete [] fShuffledIndexes; fShuffledIndexes=NULL;
7382279b 335}
336
7b5556ef 337//-----------------------------------------------------------------------
338AliFlowTrackSimple* AliFlowEventSimple::MakeNewTrack()
339{
340 AliFlowTrackSimple *t=dynamic_cast<AliFlowTrackSimple *>(fTrackCollection->RemoveAt(fNumberOfTracks));
341 if( !t ) { // If there was no track at the end of the list then create a new track
342 t=new AliFlowTrackSimple();
343 }
344
345 return t;
346}
347
929098e4 348//-----------------------------------------------------------------------
349AliFlowVector AliFlowEventSimple::GetQ(Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights)
f1d945a1 350{
929098e4 351 // calculate Q-vector in harmonic n without weights (default harmonic n=2)
e35ddff0 352 Double_t dQX = 0.;
353 Double_t dQY = 0.;
354 AliFlowVector vQ;
355 vQ.Set(0.,0.);
929098e4 356
9825d4a9 357 Int_t iOrder = n;
8eca5d19 358 Double_t sumOfWeights = 0.;
44e060e0 359 Double_t dPhi = 0.;
360 Double_t dPt = 0.;
361 Double_t dEta = 0.;
362 Double_t dWeight = 1.;
929098e4 363
26c4cbb9 364 AliFlowTrackSimple* pTrack = NULL;
929098e4 365
44e060e0 366 Int_t nBinsPhi = 0;
367 Double_t dBinWidthPt = 0.;
368 Double_t dPtMin = 0.;
369 Double_t dBinWidthEta = 0.;
370 Double_t dEtaMin = 0.;
929098e4 371
44e060e0 372 Double_t wPhi = 1.; // weight Phi
373 Double_t wPt = 1.; // weight Pt
374 Double_t wEta = 1.; // weight Eta
929098e4 375
03a02aca 376 TH1F *phiWeights = NULL;
ae733b3b 377 TH1D *ptWeights = NULL;
03a02aca 378 TH1D *etaWeights = NULL;
929098e4 379
03a02aca 380 if(weightsList)
26c4cbb9 381 {
929098e4 382 if(usePhiWeights)
383 {
384 phiWeights = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights"));
385 if(phiWeights) nBinsPhi = phiWeights->GetNbinsX();
386 }
387 if(usePtWeights)
03a02aca 388 {
929098e4 389 ptWeights = dynamic_cast<TH1D *>(weightsList->FindObject("pt_weights"));
390 if(ptWeights)
391 {
392 dBinWidthPt = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
393 dPtMin = (ptWeights->GetXaxis())->GetXmin();
394 }
395 }
396 if(useEtaWeights)
03a02aca 397 {
929098e4 398 etaWeights = dynamic_cast<TH1D *>(weightsList->FindObject("eta_weights"));
399 if(etaWeights)
400 {
401 dBinWidthEta = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
402 dEtaMin = (etaWeights->GetXaxis())->GetXmin();
403 }
404 }
03a02aca 405 } // end of if(weightsList)
929098e4 406
407 // loop over tracks
408 for(Int_t i=0; i<fNumberOfTracks; i++)
26c4cbb9 409 {
124fb262 410 pTrack = (AliFlowTrackSimple*)fTrackCollection->At(i);
929098e4 411 if(pTrack)
f1d945a1 412 {
929098e4 413 if(pTrack->InRPSelection())
414 {
415 dPhi = pTrack->Phi();
416 dPt = pTrack->Pt();
417 dEta = pTrack->Eta();
701f71c1 418 dWeight = pTrack->Weight();
929098e4 419
420 // determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
421 if(phiWeights && nBinsPhi)
422 {
423 wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi())));
424 }
425 // determine v'(pt) weight:
426 if(ptWeights && dBinWidthPt)
427 {
428 wPt=ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-dPtMin)/dBinWidthPt)));
429 }
430 // determine v'(eta) weight:
431 if(etaWeights && dBinWidthEta)
432 {
433 wEta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-dEtaMin)/dBinWidthEta)));
434 }
435
436 // building up the weighted Q-vector:
44e060e0 437 dQX += dWeight*wPhi*wPt*wEta*TMath::Cos(iOrder*dPhi);
438 dQY += dWeight*wPhi*wPt*wEta*TMath::Sin(iOrder*dPhi);
929098e4 439
440 // weighted multiplicity:
8eca5d19 441 sumOfWeights += dWeight*wPhi*wPt*wEta;
929098e4 442
443 } // end of if (pTrack->InRPSelection())
444 } // end of if (pTrack)
445 else
446 {
447 cerr << "no particle!!!"<<endl;
448 }
ae733b3b 449 } // loop over particles
929098e4 450
e35ddff0 451 vQ.Set(dQX,dQY);
8eca5d19 452 vQ.SetMult(sumOfWeights);
929098e4 453
e35ddff0 454 return vQ;
929098e4 455
5fef318d 456}
457
929098e4 458//-----------------------------------------------------------------------
459void AliFlowEventSimple::Get2Qsub(AliFlowVector* Qarray, Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights)
395fadba 460{
929098e4 461
462 // calculate Q-vector in harmonic n without weights (default harmonic n=2)
395fadba 463 Double_t dQX = 0.;
464 Double_t dQY = 0.;
929098e4 465
395fadba 466 Int_t iOrder = n;
8eca5d19 467 Double_t sumOfWeights = 0.;
395fadba 468 Double_t dPhi = 0.;
469 Double_t dPt = 0.;
470 Double_t dEta = 0.;
44e060e0 471 Double_t dWeight = 1.;
929098e4 472
395fadba 473 AliFlowTrackSimple* pTrack = NULL;
929098e4 474
44e060e0 475 Int_t iNbinsPhiSub0 = 0;
476 Int_t iNbinsPhiSub1 = 0;
395fadba 477 Double_t dBinWidthPt = 0.;
478 Double_t dPtMin = 0.;
479 Double_t dBinWidthEta= 0.;
480 Double_t dEtaMin = 0.;
929098e4 481
482 Double_t dWphi = 1.; // weight Phi
483 Double_t dWpt = 1.; // weight Pt
484 Double_t dWeta = 1.; // weight Eta
485
44e060e0 486 TH1F* phiWeightsSub0 = NULL;
487 TH1F* phiWeightsSub1 = NULL;
29195b69 488 TH1D* ptWeights = NULL;
489 TH1D* etaWeights = NULL;
929098e4 490
491 if(weightsList)
492 {
493 if(usePhiWeights)
494 {
44e060e0 495 phiWeightsSub0 = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights_sub0"));
44e060e0 496 if(phiWeightsSub0) {
13ff9ccd 497 iNbinsPhiSub0 = phiWeightsSub0->GetNbinsX();
44e060e0 498 }
13ff9ccd 499 phiWeightsSub1 = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights_sub1"));
44e060e0 500 if(phiWeightsSub1) {
13ff9ccd 501 iNbinsPhiSub1 = phiWeightsSub1->GetNbinsX();
29195b69 502 }
929098e4 503 }
504 if(usePtWeights)
505 {
29195b69 506 ptWeights = dynamic_cast<TH1D *>(weightsList->FindObject("pt_weights"));
929098e4 507 if(ptWeights)
508 {
509 dBinWidthPt = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
510 dPtMin = (ptWeights->GetXaxis())->GetXmin();
511 }
512 }
513 if(useEtaWeights)
514 {
515 etaWeights = dynamic_cast<TH1D *>(weightsList->FindObject("eta_weights"));
516 if(etaWeights)
517 {
518 dBinWidthEta = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
519 dEtaMin = (etaWeights->GetXaxis())->GetXmin();
520 }
521 }
395fadba 522 } // end of if(weightsList)
929098e4 523
b125a454 524 //loop over the two subevents
929098e4 525 for (Int_t s=0; s<2; s++)
526 {
527 // loop over tracks
528 for(Int_t i=0; i<fNumberOfTracks; i++)
529 {
124fb262 530 pTrack = (AliFlowTrackSimple*)fTrackCollection->At(i);
929098e4 531 if(pTrack)
532 {
533 if(pTrack->InRPSelection())
534 {
535 if (pTrack->InSubevent(s))
536 {
44e060e0 537 dPhi = pTrack->Phi();
538 dPt = pTrack->Pt();
539 dEta = pTrack->Eta();
540 dWeight = pTrack->Weight();
929098e4 541
542 // determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
13ff9ccd 543 //subevent 0
544 if(s == 0) {
545 if(phiWeightsSub0 && iNbinsPhiSub0) {
546 Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNbinsPhiSub0/TMath::TwoPi()));
547 //use the phi value at the center of the bin
548 dPhi = phiWeightsSub0->GetBinCenter(phiBin);
549 dWphi = phiWeightsSub0->GetBinContent(phiBin);
550 }
551 }
552 //subevent 1
553 else if (s == 1) {
554 if(phiWeightsSub1 && iNbinsPhiSub1) {
555 Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNbinsPhiSub1/TMath::TwoPi()));
556 //use the phi value at the center of the bin
557 dPhi = phiWeightsSub1->GetBinCenter(phiBin);
558 dWphi = phiWeightsSub1->GetBinContent(phiBin);
44e060e0 559 }
560 }
13ff9ccd 561
929098e4 562 // determine v'(pt) weight:
563 if(ptWeights && dBinWidthPt)
564 {
565 dWpt=ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-dPtMin)/dBinWidthPt)));
566 }
44e060e0 567
929098e4 568 // determine v'(eta) weight:
569 if(etaWeights && dBinWidthEta)
570 {
571 dWeta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-dEtaMin)/dBinWidthEta)));
572 }
573
574 // building up the weighted Q-vector:
44e060e0 575 dQX += dWeight*dWphi*dWpt*dWeta*TMath::Cos(iOrder*dPhi);
576 dQY += dWeight*dWphi*dWpt*dWeta*TMath::Sin(iOrder*dPhi);
929098e4 577
578 // weighted multiplicity:
8eca5d19 579 sumOfWeights+=dWeight*dWphi*dWpt*dWeta;
929098e4 580
581 } // end of subevent
582 } // end of if (pTrack->InRPSelection())
29195b69 583 } // end of if (pTrack)
929098e4 584 else
585 {
586 cerr << "no particle!!!"<<endl;
587 }
29195b69 588 } // loop over particles
589 Qarray[s].Set(dQX,dQY);
8eca5d19 590 Qarray[s].SetMult(sumOfWeights);
29195b69 591 //reset
8eca5d19 592 sumOfWeights = 0.;
29195b69 593 dQX = 0.;
594 dQY = 0.;
595 }
929098e4 596
395fadba 597}
598
599
929098e4 600//-----------------------------------------------------------------------
c076fda8 601void AliFlowEventSimple::Print(Option_t *option) const
602{
603 // -*-*-*-*-*Print some global quantities for this histogram collection class *-*-*-*-*-*-*-*
604 // ===============================================
605 // printf( "TH1.Print Name = %s, Entries= %d, Total sum= %g\n",GetName(),Int_t(fEntries),GetSumOfWeights());
df23c5ae 606 printf( "Class.Print Name = %s, #tracks= %d, Number of RPs= %d, Number of POIs= %d, MC EventPlaneAngle= %f\n",
607 GetName(),fNumberOfTracks, fNumberOfRPs, fNumberOfPOIs, fMCReactionPlaneAngle );
c076fda8 608
b4dba88d 609 TString optionstr(option);
610 if (!optionstr.Contains("all")) return;
929098e4 611 if (fTrackCollection)
612 {
c076fda8 613 fTrackCollection->Print(option);
614 }
929098e4 615 else
616 {
c076fda8 617 printf( "Empty track collection \n");
618 }
619}
620
929098e4 621//-----------------------------------------------------------------------
622void AliFlowEventSimple::Browse(TBrowser *b)
c076fda8 623{
cc0afcfc 624 //browse in TBrowser
c076fda8 625 if (!b) return;
929098e4 626 if (!fNumberOfTracksWrap)
627 {
c076fda8 628 fNumberOfTracksWrap = new TParameter<int>("fNumberOfTracks", fNumberOfTracks);
629 b->Add(fNumberOfTracksWrap);
630 }
bc231a12 631 if (!fNumberOfRPsWrap)
929098e4 632 {
bc231a12 633 fNumberOfRPsWrap = new TParameter<int>("fNumberOfRPs", fNumberOfRPs);
634 b->Add(fNumberOfRPsWrap);
c076fda8 635 }
df23c5ae 636 if (!fNumberOfPOIsWrap)
637 {
638 fNumberOfPOIsWrap = new TParameter<int>("fNumberOfPOIs", fNumberOfPOIs);
639 b->Add(fNumberOfPOIsWrap);
640 }
929098e4 641 if (!fMCReactionPlaneAngleWrap)
642 {
7183fe85 643 fMCReactionPlaneAngleWrap = new TParameter<double>(" fMCReactionPlaneAngle", fMCReactionPlaneAngle);
a12990bb 644 b->Add( fMCReactionPlaneAngleWrap);
645 }
c076fda8 646 if (fTrackCollection) b->Add(fTrackCollection,"AliFlowTracksSimple");
647}
648
929098e4 649//-----------------------------------------------------------------------
650AliFlowEventSimple::AliFlowEventSimple( TTree* inputTree,
651 const AliFlowTrackSimpleCuts* rpCuts,
652 const AliFlowTrackSimpleCuts* poiCuts):
653 fTrackCollection(NULL),
85d2ee8d 654 fReferenceMultiplicity(0),
929098e4 655 fNumberOfTracks(0),
bc231a12 656 fNumberOfRPs(0),
df23c5ae 657 fNumberOfPOIs(0),
8d1c89ad 658 fUseGlauberMCSymmetryPlanes(kFALSE),
659 fUseExternalSymmetryPlanes(kFALSE),
660 fPsi1(0.),
661 fPsi2(0.),
662 fPsi3(0.),
663 fPsi4(0.),
664 fPsi5(0.),
665 fPsi1Psi3(0x0),
666 fPsi2Psi4(0x0),
667 fPsi3Psi5(0x0),
929098e4 668 fMCReactionPlaneAngle(0.),
669 fMCReactionPlaneAngleIsSet(kFALSE),
bc231a12 670 fAfterBurnerPrecision(0.001),
9fa64edc 671 fUserModified(kFALSE),
929098e4 672 fNumberOfTracksWrap(NULL),
bc231a12 673 fNumberOfRPsWrap(NULL),
df23c5ae 674 fNumberOfPOIsWrap(NULL),
ad4d7e4e 675 fMCReactionPlaneAngleWrap(NULL),
676 fShuffledIndexes(NULL),
677 fShuffleTracks(kFALSE)
929098e4 678{
679 //constructor, fills the event from a TTree of kinematic.root files
680 //applies RP and POI cuts, tags the tracks
681
682 Int_t numberOfInputTracks = inputTree->GetEntries() ;
683 fTrackCollection = new TObjArray(numberOfInputTracks/2);
684
685 TParticle* pParticle = new TParticle();
686 inputTree->SetBranchAddress("Particles",&pParticle);
687
929098e4 688 for (Int_t i=0; i<numberOfInputTracks; i++)
689 {
690 inputTree->GetEntry(i); //get input particle
691
692 if (!pParticle) continue; //no particle
693 if (!pParticle->IsPrimary()) continue;
694
7382279b 695 Bool_t rpOK = rpCuts->PassesCuts(pParticle);
696 Bool_t poiOK = poiCuts->PassesCuts(pParticle);
929098e4 697
7382279b 698 if (rpOK || poiOK)
929098e4 699 {
700 AliFlowTrackSimple* pTrack = new AliFlowTrackSimple(pParticle);
701
702 //marking the particles used for int. flow:
7382279b 703 if(rpOK)
929098e4 704 {
705 pTrack->SetForRPSelection(kTRUE);
bc231a12 706 fNumberOfRPs++;
df23c5ae 707 cout<<"fNumberOfRPs = "<<fNumberOfRPs<<endl;
929098e4 708 }
709 //marking the particles used for diff. flow:
7382279b 710 if(poiOK)
929098e4 711 {
712 pTrack->SetForPOISelection(kTRUE);
df23c5ae 713 fNumberOfPOIs++;
714 cout<<"fNumberOfPOIs = "<<fNumberOfPOIs<<endl;
929098e4 715 }
716 //adding a particles which were used either for int. or diff. flow to the list
124fb262 717 AddTrack(pTrack);
929098e4 718 }
719 }//for i
720 delete pParticle;
721}
722
723//_____________________________________________________________________________
724void AliFlowEventSimple::CloneTracks(Int_t n)
725{
726 //clone every track n times to add non-flow
34b15925 727 if (n<=0) return; //no use to clone stuff zero or less times
bc231a12 728 Int_t ntracks = fNumberOfTracks;
729 fTrackCollection->Expand((n+1)*fNumberOfTracks);
730 for (Int_t i=0; i<n; i++)
929098e4 731 {
bc231a12 732 for (Int_t itrack=0; itrack<ntracks; itrack++)
929098e4 733 {
734 AliFlowTrackSimple* track = dynamic_cast<AliFlowTrackSimple*>(fTrackCollection->At(itrack));
735 if (!track) continue;
bc231a12 736 AddTrack(static_cast<AliFlowTrackSimple*>(track->Clone()));
929098e4 737 }
738 }
9fa64edc 739 SetUserModified();
929098e4 740}
7382279b 741
742//_____________________________________________________________________________
743void AliFlowEventSimple::ResolutionPt(Double_t res)
744{
745 //smear pt of all tracks by gaussian with sigma=res
746 for (Int_t i=0; i<fNumberOfTracks; i++)
747 {
748 AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
124fb262 749 if (track) track->ResolutionPt(res);
7382279b 750 }
9fa64edc 751 SetUserModified();
7382279b 752}
124fb262 753
754//_____________________________________________________________________________
34b15925 755void AliFlowEventSimple::TagSubeventsInEta( Double_t etaMinA,
756 Double_t etaMaxA,
757 Double_t etaMinB,
758 Double_t etaMaxB )
124fb262 759{
760 //Flag two subevents in given eta ranges
761 for (Int_t i=0; i<fNumberOfTracks; i++)
762 {
763 AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
3ca4688a 764 if (!track) continue;
765 track->ResetSubEventTags();
124fb262 766 Double_t eta=track->Eta();
767 if (eta >= etaMinA && eta <= etaMaxA) track->SetForSubevent(0);
768 if (eta >= etaMinB && eta <= etaMaxB) track->SetForSubevent(1);
769 }
770}
771
076df7bf 772//_____________________________________________________________________________
773void AliFlowEventSimple::TagSubeventsByCharge()
774{
775 //Flag two subevents in given eta ranges
776 for (Int_t i=0; i<fNumberOfTracks; i++)
777 {
778 AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
3ca4688a 779 if (!track) continue;
780 track->ResetSubEventTags();
076df7bf 781 Int_t charge=track->Charge();
782 if (charge<0) track->SetForSubevent(0);
783 if (charge>0) track->SetForSubevent(1);
784 }
785}
786
bc231a12 787//_____________________________________________________________________________
34b15925 788void AliFlowEventSimple::AddV1( Double_t v1 )
bc231a12 789{
790 //add v2 to all tracks wrt the reaction plane angle
791 for (Int_t i=0; i<fNumberOfTracks; i++)
792 {
793 AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
8d1c89ad 794 if (track) {
795 if((fUseExternalSymmetryPlanes)||(fUseGlauberMCSymmetryPlanes))
796 track->AddV1(v1, fPsi1, fAfterBurnerPrecision);
797 else
798 track->AddV1(v1, fMCReactionPlaneAngle, fAfterBurnerPrecision);
799 }
bc231a12 800 }
9fa64edc 801 SetUserModified();
bc231a12 802}
803
124fb262 804//_____________________________________________________________________________
34b15925 805void AliFlowEventSimple::AddV2( Double_t v2 )
124fb262 806{
244c607a 807 //add v2 to all tracks wrt the reaction plane angle
124fb262 808 for (Int_t i=0; i<fNumberOfTracks; i++)
809 {
810 AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
8d1c89ad 811 if (track) {
812 if((fUseExternalSymmetryPlanes)||(fUseGlauberMCSymmetryPlanes))
813 track->AddV2(v2, fPsi2, fAfterBurnerPrecision);
814 else
815 track->AddV2(v2, fMCReactionPlaneAngle, fAfterBurnerPrecision);
816 }
bc231a12 817 }
9fa64edc 818 SetUserModified();
bc231a12 819}
820
821//_____________________________________________________________________________
54089829 822void AliFlowEventSimple::AddV3( Double_t v3 )
823{
824 //add v3 to all tracks wrt the reaction plane angle
825 for (Int_t i=0; i<fNumberOfTracks; i++)
826 {
827 AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
8d1c89ad 828 if(track) {
829 if((fUseExternalSymmetryPlanes)||(fUseGlauberMCSymmetryPlanes))
830 track->AddV3(v3, fPsi3, fAfterBurnerPrecision);
831 else
832 track->AddV3(v3, fMCReactionPlaneAngle, fAfterBurnerPrecision);
833 }
54089829 834 }
835 SetUserModified();
836}
837
838//_____________________________________________________________________________
34b15925 839void AliFlowEventSimple::AddV4( Double_t v4 )
bc231a12 840{
841 //add v4 to all tracks wrt the reaction plane angle
842 for (Int_t i=0; i<fNumberOfTracks; i++)
843 {
844 AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
8d1c89ad 845 if(track) {
846 if((fUseExternalSymmetryPlanes)||(fUseGlauberMCSymmetryPlanes))
847 track->AddV4(v4, fPsi4, fAfterBurnerPrecision);
848 else
849 track->AddV4(v4, fMCReactionPlaneAngle, fAfterBurnerPrecision);
850 }
bc231a12 851 }
9fa64edc 852 SetUserModified();
bc231a12 853}
854
855//_____________________________________________________________________________
1a80f9f6 856void AliFlowEventSimple::AddV5( Double_t v5 )
857{
858 //add v4 to all tracks wrt the reaction plane angle
859 for (Int_t i=0; i<fNumberOfTracks; i++)
860 {
861 AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
8d1c89ad 862 if(track) {
863 if((fUseExternalSymmetryPlanes)||(fUseGlauberMCSymmetryPlanes))
864 track->AddV5(v5, fPsi5, fAfterBurnerPrecision);
865 else
866 track->AddV5(v5, fMCReactionPlaneAngle, fAfterBurnerPrecision);
867 }
1a80f9f6 868 }
869 SetUserModified();
870}
871
872//_____________________________________________________________________________
873void AliFlowEventSimple::AddFlow( Double_t v1, Double_t v2, Double_t v3, Double_t v4, Double_t v5,
874 Double_t rp1, Double_t rp2, Double_t rp3, Double_t rp4, Double_t rp5 )
875{
876 //add flow to all tracks wrt the reaction plane angle, for all harmonic separate angle
877 for (Int_t i=0; i<fNumberOfTracks; i++)
878 {
879 AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
880 if (track) track->AddFlow(v1,v2,v3,v4,v5,rp1,rp2,rp3,rp4,rp5,fAfterBurnerPrecision);
881 }
882 SetUserModified();
883}
884
885//_____________________________________________________________________________
886void AliFlowEventSimple::AddFlow( Double_t v1, Double_t v2, Double_t v3, Double_t v4, Double_t v5 )
bc231a12 887{
888 //add flow to all tracks wrt the reaction plane angle
889 for (Int_t i=0; i<fNumberOfTracks; i++)
890 {
891 AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
1a80f9f6 892 if (track) track->AddFlow(v1,v2,v3,v4,v5,fMCReactionPlaneAngle, fAfterBurnerPrecision);
124fb262 893 }
9fa64edc 894 SetUserModified();
124fb262 895}
896
99ff691b 897//_____________________________________________________________________________
898void AliFlowEventSimple::AddV2( TF1* ptDepV2 )
899{
900 //add v2 to all tracks wrt the reaction plane angle
901 for (Int_t i=0; i<fNumberOfTracks; i++)
902 {
903 AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
2a745a5f 904 if (!track) continue;
99ff691b 905 Double_t v2 = ptDepV2->Eval(track->Pt());
2a745a5f 906 track->AddV2(v2, fMCReactionPlaneAngle, fAfterBurnerPrecision);
99ff691b 907 }
908 SetUserModified();
909}
910
dd91b595 911//_____________________________________________________________________________
cc0afcfc 912void AliFlowEventSimple::TagRP( const AliFlowTrackSimpleCuts* cuts )
dd91b595 913{
914 //tag tracks as reference particles (RPs)
915 for (Int_t i=0; i<fNumberOfTracks; i++)
916 {
917 AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
918 if (!track) continue;
701f71c1 919 Bool_t pass=cuts->PassesCuts(track);
e0da66a2 920 Bool_t rpTrack=track->InRPSelection();
701f71c1 921 if (pass)
e0da66a2 922 {
923 if (!rpTrack) fNumberOfRPs++; //only increase if not already tagged
924 }
701f71c1 925 else
e0da66a2 926 {
927 if (rpTrack) fNumberOfRPs--; //only decrease if detagging
928 }
929 track->SetForRPSelection(pass);
dd91b595 930 }
931}
932
933//_____________________________________________________________________________
cc0afcfc 934void AliFlowEventSimple::TagPOI( const AliFlowTrackSimpleCuts* cuts )
dd91b595 935{
936 //tag tracks as particles of interest (POIs)
937 for (Int_t i=0; i<fNumberOfTracks; i++)
938 {
939 AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
940 if (!track) continue;
701f71c1 941 Bool_t pass=cuts->PassesCuts(track);
df23c5ae 942 Bool_t poiTrack=track->InPOISelection();
943 if (pass)
944 {
945 if (!poiTrack) fNumberOfPOIs++; //only increase if not already tagged
946 }
947 else
948 {
949 if (poiTrack) fNumberOfPOIs--; //only decrease if detagging
950 }
701f71c1 951 track->SetForPOISelection(pass);
dd91b595 952 }
953}
954
34b15925 955//_____________________________________________________________________________
956void AliFlowEventSimple::DefineDeadZone( Double_t etaMin,
957 Double_t etaMax,
958 Double_t phiMin,
959 Double_t phiMax )
960{
961 //mark tracks in given eta-phi region as dead
962 //by resetting the flow bits
963 for (Int_t i=0; i<fNumberOfTracks; i++)
964 {
965 AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
966 Double_t eta = track->Eta();
967 Double_t phi = track->Phi();
968 if (eta>etaMin && eta<etaMax && phi>phiMin && phi<phiMax)
1a43aca9 969 {
df23c5ae 970 if (track->InRPSelection()) {fNumberOfRPs--;}
971 if (track->InPOISelection()) {fNumberOfPOIs--;}
34b15925 972 track->ResetFlowTags();
1a43aca9 973 }
34b15925 974 }
975}
976
977//_____________________________________________________________________________
978Int_t AliFlowEventSimple::CleanUpDeadTracks()
979{
980 //remove tracks that have no flow tags set and cleanup the container
65ef7d1a 981 //returns number of cleaned tracks
982 Int_t ncleaned=0;
34b15925 983 for (Int_t i=0; i<fNumberOfTracks; i++)
984 {
985 AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
cd62a2a8 986 if (!track) continue;
7d27a354 987 if (track->IsDead()) {delete track;track=NULL;ncleaned++;}
34b15925 988 }
989 fTrackCollection->Compress(); //clean up empty slots
7d27a354 990 fNumberOfTracks-=ncleaned; //update number of tracks
37b19b98 991 delete [] fShuffledIndexes; fShuffledIndexes=NULL;
65ef7d1a 992 return ncleaned;
34b15925 993}
99ff691b 994
995//_____________________________________________________________________________
996TF1* AliFlowEventSimple::SimplePtDepV2()
997{
998 //return a standard pt dependent v2 formula, user has to clean up!
999 return new TF1("StandardPtDepV2","((x<1.0)*(0.05/1.0)*x+(x>=1.0)*0.05)");
1000}
1001
1002//_____________________________________________________________________________
1003TF1* AliFlowEventSimple::SimplePtSpectrum()
1004{
1005 //return a standard pt spectrum, user has to clean up!
1006 return new TF1("StandardPtSpectrum","x*TMath::Exp(-pow(0.13957*0.13957+x*x,0.5)/0.4)",0.1,10.);
1007}
7d27a354 1008
1009//_____________________________________________________________________________
1010void AliFlowEventSimple::ClearFast()
1011{
1012 //clear the counter without deleting allocated objects so they can be reused
1013 fReferenceMultiplicity = 0;
1014 fNumberOfTracks = 0;
1015 fNumberOfRPs = 0;
df23c5ae 1016 fNumberOfPOIs = 0;
7d27a354 1017 fMCReactionPlaneAngle = 0.0;
1018 fMCReactionPlaneAngleIsSet = kFALSE;
1019 fAfterBurnerPrecision = 0.001;
1020 fUserModified = kFALSE;
1021}