1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /*****************************************************************
17 AliFlowEventSimple: A simple event
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 *****************************************************************/
26 #include "Riostream.h"
27 #include "TObjArray.h"
31 #include "TParticle.h"
36 #include "TParameter.h"
38 #include "AliFlowVector.h"
39 #include "AliFlowTrackSimple.h"
40 #include "AliFlowEventSimple.h"
41 #include "AliFlowTrackSimpleCuts.h"
44 ClassImp(AliFlowEventSimple)
46 //-----------------------------------------------------------------------
47 AliFlowEventSimple::AliFlowEventSimple():
48 fTrackCollection(NULL),
50 fEventNSelTracksRP(0),
51 fMCReactionPlaneAngle(0.),
52 fMCReactionPlaneAngleIsSet(kFALSE),
53 fNumberOfTracksWrap(NULL),
54 fEventNSelTracksRPWrap(NULL),
55 fMCReactionPlaneAngleWrap(NULL)
57 cout << "AliFlowEventSimple: Default constructor to be used only by root for io" << endl;
60 //-----------------------------------------------------------------------
61 AliFlowEventSimple::AliFlowEventSimple(Int_t aLenght):
62 fTrackCollection(NULL),
64 fEventNSelTracksRP(0),
65 fMCReactionPlaneAngle(0.),
66 fMCReactionPlaneAngleIsSet(kFALSE),
67 fNumberOfTracksWrap(NULL),
68 fEventNSelTracksRPWrap(NULL),
69 fMCReactionPlaneAngleWrap(NULL)
72 fTrackCollection = new TObjArray(aLenght);
75 //-----------------------------------------------------------------------
76 AliFlowEventSimple::AliFlowEventSimple(const AliFlowEventSimple& anEvent):
78 fTrackCollection(NULL),
79 fNumberOfTracks(anEvent.fNumberOfTracks),
80 fEventNSelTracksRP(anEvent.fEventNSelTracksRP),
81 fMCReactionPlaneAngle(anEvent.fMCReactionPlaneAngle),
82 fMCReactionPlaneAngleIsSet(anEvent.fMCReactionPlaneAngleIsSet),
83 fNumberOfTracksWrap(anEvent.fNumberOfTracksWrap),
84 fEventNSelTracksRPWrap(anEvent.fEventNSelTracksRPWrap),
85 fMCReactionPlaneAngleWrap(anEvent.fMCReactionPlaneAngleWrap)
88 fTrackCollection = (TObjArray*)anEvent.Clone(); //deep copy
91 //-----------------------------------------------------------------------
92 AliFlowEventSimple& AliFlowEventSimple::operator=(const AliFlowEventSimple& anEvent)
95 delete fTrackCollection;
96 fTrackCollection = (TObjArray*)anEvent.Clone(); //deep copy
97 fNumberOfTracks = anEvent.fNumberOfTracks;
98 fEventNSelTracksRP = anEvent.fEventNSelTracksRP;
99 fMCReactionPlaneAngle = anEvent.fMCReactionPlaneAngle;
100 fMCReactionPlaneAngleIsSet = anEvent.fMCReactionPlaneAngleIsSet;
101 fNumberOfTracksWrap = anEvent.fNumberOfTracksWrap;
102 fEventNSelTracksRPWrap = anEvent.fEventNSelTracksRPWrap;
103 fMCReactionPlaneAngleWrap=anEvent.fMCReactionPlaneAngleWrap;
107 //-----------------------------------------------------------------------
108 AliFlowEventSimple::~AliFlowEventSimple()
111 if (fTrackCollection) fTrackCollection->Delete();
112 delete fTrackCollection;
113 if (fNumberOfTracksWrap) delete fNumberOfTracksWrap;
114 if (fEventNSelTracksRPWrap) delete fEventNSelTracksRPWrap;
115 if (fMCReactionPlaneAngleWrap) delete fMCReactionPlaneAngleWrap;
118 //-----------------------------------------------------------------------
119 AliFlowTrackSimple* AliFlowEventSimple::GetTrack(Int_t i)
121 //get track i from collection
122 if (i>=fNumberOfTracks) return NULL;
123 AliFlowTrackSimple* pTrack = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i)) ;
127 //-----------------------------------------------------------------------
128 void AliFlowEventSimple::AddTrack( AliFlowTrackSimple* track )
131 if (!fTrackCollection) return;
132 fTrackCollection->AddLast(track);
136 //-----------------------------------------------------------------------
137 AliFlowVector AliFlowEventSimple::GetQ(Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights)
139 // calculate Q-vector in harmonic n without weights (default harmonic n=2)
146 Double_t iUsedTracks = 0;
151 AliFlowTrackSimple* pTrack = NULL;
154 Double_t dBinWidthPt=0.;
156 Double_t dBinWidthEta=0.;
159 Double_t wPhi=1.; // weight Phi
160 Double_t wPt=1.; // weight Pt
161 Double_t wEta=1.; // weight Eta
163 TH1F *phiWeights = NULL;
164 TH1D *ptWeights = NULL;
165 TH1D *etaWeights = NULL;
171 phiWeights = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights"));
172 if(phiWeights) nBinsPhi = phiWeights->GetNbinsX();
176 ptWeights = dynamic_cast<TH1D *>(weightsList->FindObject("pt_weights"));
179 dBinWidthPt = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
180 dPtMin = (ptWeights->GetXaxis())->GetXmin();
185 etaWeights = dynamic_cast<TH1D *>(weightsList->FindObject("eta_weights"));
188 dBinWidthEta = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
189 dEtaMin = (etaWeights->GetXaxis())->GetXmin();
192 } // end of if(weightsList)
195 for(Int_t i=0; i<fNumberOfTracks; i++)
197 pTrack = (AliFlowTrackSimple*)fTrackCollection->At(i);
200 if(pTrack->InRPSelection())
202 dPhi = pTrack->Phi();
204 dEta = pTrack->Eta();
206 // determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
207 if(phiWeights && nBinsPhi)
209 wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi())));
211 // determine v'(pt) weight:
212 if(ptWeights && dBinWidthPt)
214 wPt=ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-dPtMin)/dBinWidthPt)));
216 // determine v'(eta) weight:
217 if(etaWeights && dBinWidthEta)
219 wEta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-dEtaMin)/dBinWidthEta)));
222 // building up the weighted Q-vector:
223 dQX += wPhi*wPt*wEta*TMath::Cos(iOrder*dPhi);
224 dQY += wPhi*wPt*wEta*TMath::Sin(iOrder*dPhi);
226 // weighted multiplicity:
227 iUsedTracks+=wPhi*wPt*wEta;
229 } // end of if (pTrack->InRPSelection())
230 } // end of if (pTrack)
233 cerr << "no particle!!!"<<endl;
235 } // loop over particles
238 vQ.SetMult(iUsedTracks);
244 //-----------------------------------------------------------------------
245 void AliFlowEventSimple::Get2Qsub(AliFlowVector* Qarray, Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights)
248 // calculate Q-vector in harmonic n without weights (default harmonic n=2)
253 Double_t iUsedTracks = 0;
258 AliFlowTrackSimple* pTrack = NULL;
261 Double_t dBinWidthPt = 0.;
262 Double_t dPtMin = 0.;
263 Double_t dBinWidthEta= 0.;
264 Double_t dEtaMin = 0.;
266 Double_t dWphi = 1.; // weight Phi
267 Double_t dWpt = 1.; // weight Pt
268 Double_t dWeta = 1.; // weight Eta
270 TH1F* phiWeights = NULL;
271 TH1D* ptWeights = NULL;
272 TH1D* etaWeights = NULL;
278 phiWeights = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights"));
281 iNbinsPhi = phiWeights->GetNbinsX();
286 ptWeights = dynamic_cast<TH1D *>(weightsList->FindObject("pt_weights"));
289 dBinWidthPt = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
290 dPtMin = (ptWeights->GetXaxis())->GetXmin();
295 etaWeights = dynamic_cast<TH1D *>(weightsList->FindObject("eta_weights"));
298 dBinWidthEta = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
299 dEtaMin = (etaWeights->GetXaxis())->GetXmin();
302 } // end of if(weightsList)
304 //loop over the two subevents
305 for (Int_t s=0; s<2; s++)
308 for(Int_t i=0; i<fNumberOfTracks; i++)
310 pTrack = (AliFlowTrackSimple*)fTrackCollection->At(i);
313 if(pTrack->InRPSelection())
315 if (pTrack->InSubevent(s))
317 dPhi = pTrack->Phi();
319 dEta = pTrack->Eta();
321 // determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
322 if(phiWeights && iNbinsPhi)
324 dWphi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*iNbinsPhi/TMath::TwoPi())));
326 // determine v'(pt) weight:
327 if(ptWeights && dBinWidthPt)
329 dWpt=ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-dPtMin)/dBinWidthPt)));
331 // determine v'(eta) weight:
332 if(etaWeights && dBinWidthEta)
334 dWeta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-dEtaMin)/dBinWidthEta)));
337 // building up the weighted Q-vector:
338 dQX += dWphi*dWpt*dWeta*TMath::Cos(iOrder*dPhi);
339 dQY += dWphi*dWpt*dWeta*TMath::Sin(iOrder*dPhi);
341 // weighted multiplicity:
342 iUsedTracks+=dWphi*dWpt*dWeta;
345 } // end of if (pTrack->InRPSelection())
346 } // end of if (pTrack)
349 cerr << "no particle!!!"<<endl;
351 } // loop over particles
352 Qarray[s].Set(dQX,dQY);
353 Qarray[s].SetMult(iUsedTracks);
363 //-----------------------------------------------------------------------
364 void AliFlowEventSimple::Print(Option_t *option) const
366 // -*-*-*-*-*Print some global quantities for this histogram collection class *-*-*-*-*-*-*-*
367 // ===============================================
368 // printf( "TH1.Print Name = %s, Entries= %d, Total sum= %g\n",GetName(),Int_t(fEntries),GetSumOfWeights());
369 printf( "Class.Print Name = %s, Total number of tracks= %d, Number of selected tracks= %d, MC EventPlaneAngle= %f\n",
370 GetName(),fNumberOfTracks, fEventNSelTracksRP, fMCReactionPlaneAngle );
372 if (fTrackCollection)
374 fTrackCollection->Print(option);
378 printf( "Empty track collection \n");
382 //-----------------------------------------------------------------------
383 void AliFlowEventSimple::Browse(TBrowser *b)
386 if (!fNumberOfTracksWrap)
388 fNumberOfTracksWrap = new TParameter<int>("fNumberOfTracks", fNumberOfTracks);
389 b->Add(fNumberOfTracksWrap);
391 if (!fEventNSelTracksRPWrap)
393 fEventNSelTracksRPWrap = new TParameter<int>("fEventNSelTracksRP", fEventNSelTracksRP);
394 b->Add(fEventNSelTracksRPWrap);
396 if (!fMCReactionPlaneAngleWrap)
398 fMCReactionPlaneAngleWrap = new TParameter<double>(" fMCReactionPlaneAngle", fMCReactionPlaneAngle);
399 b->Add( fMCReactionPlaneAngleWrap);
401 if (fTrackCollection) b->Add(fTrackCollection,"AliFlowTracksSimple");
404 //-----------------------------------------------------------------------
405 AliFlowEventSimple::AliFlowEventSimple( TTree* inputTree,
406 const AliFlowTrackSimpleCuts* rpCuts,
407 const AliFlowTrackSimpleCuts* poiCuts):
408 fTrackCollection(NULL),
410 fEventNSelTracksRP(0),
411 fMCReactionPlaneAngle(0.),
412 fMCReactionPlaneAngleIsSet(kFALSE),
413 fNumberOfTracksWrap(NULL),
414 fEventNSelTracksRPWrap(NULL),
415 fMCReactionPlaneAngleWrap(NULL)
417 //constructor, fills the event from a TTree of kinematic.root files
418 //applies RP and POI cuts, tags the tracks
420 Int_t numberOfInputTracks = inputTree->GetEntries() ;
421 fTrackCollection = new TObjArray(numberOfInputTracks/2);
423 TParticle* pParticle = new TParticle();
424 inputTree->SetBranchAddress("Particles",&pParticle);
426 Int_t iSelParticlesPOI = 0;
428 for (Int_t i=0; i<numberOfInputTracks; i++)
430 inputTree->GetEntry(i); //get input particle
432 if (!pParticle) continue; //no particle
433 if (!pParticle->IsPrimary()) continue;
435 Bool_t rpOK = rpCuts->PassesCuts(pParticle);
436 Bool_t poiOK = poiCuts->PassesCuts(pParticle);
440 AliFlowTrackSimple* pTrack = new AliFlowTrackSimple(pParticle);
442 //marking the particles used for int. flow:
445 pTrack->SetForRPSelection(kTRUE);
446 fEventNSelTracksRP++;
448 //marking the particles used for diff. flow:
451 pTrack->SetForPOISelection(kTRUE);
454 //adding a particles which were used either for int. or diff. flow to the list
461 //_____________________________________________________________________________
462 void AliFlowEventSimple::CloneTracks(Int_t n)
464 //clone every track n times to add non-flow
465 for (Int_t i=1; i<n; i++)
467 for (Int_t itrack=0; itrack<fNumberOfTracks; itrack++)
469 AliFlowTrackSimple* track = dynamic_cast<AliFlowTrackSimple*>(fTrackCollection->At(itrack));
470 if (!track) continue;
471 AddTrack(new AliFlowTrackSimple(*track));
476 //_____________________________________________________________________________
477 void AliFlowEventSimple::ResolutionPt(Double_t res)
479 //smear pt of all tracks by gaussian with sigma=res
480 for (Int_t i=0; i<fNumberOfTracks; i++)
482 AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
483 if (track) track->ResolutionPt(res);
487 //_____________________________________________________________________________
488 void AliFlowEventSimple::TagSubeventsInEta(Double_t etaMinA, Double_t etaMaxA, Double_t etaMinB, Double_t etaMaxB )
490 //Flag two subevents in given eta ranges
491 for (Int_t i=0; i<fNumberOfTracks; i++)
493 AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
494 Double_t eta=track->Eta();
495 if (eta >= etaMinA && eta <= etaMaxA) track->SetForSubevent(0);
496 if (eta >= etaMinB && eta <= etaMaxB) track->SetForSubevent(1);
500 //_____________________________________________________________________________
501 void AliFlowEventSimple::AddFlow(Double_t flow)
503 //add flow to all tracks wrt the reaction plane angle
504 for (Int_t i=0; i<fNumberOfTracks; i++)
506 AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
507 if (track) track->AddFlow(flow, fMCReactionPlaneAngle);