AliPhysics  bf041da (bf041da)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliFlowEvent.cxx
Go to the documentation of this file.
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 /*****************************************************************
17  AliFlowEvent: Event container for flow analysis
18 
19  origin: Mikolaj Krzewicki (mikolaj.krzewicki@cern.ch)
20  mods: Redmer A. Bertens (rbertens@cern.ch)
21 *****************************************************************/
22 
23 #include "Riostream.h"
24 #include "TFile.h"
25 #include "TList.h"
26 #include "TH1.h"
27 #include "TH1F.h"
28 #include "TH2F.h"
29 #include "TH3F.h"
30 #include "TArrayD.h"
31 #include "TProfile.h"
32 #include "AliMCEvent.h"
33 #include "AliMCParticle.h"
34 #include "AliCFManager.h"
35 #include "AliESDtrack.h"
36 #include "AliESDPmdTrack.h"
37 #include "AliESDEvent.h"
38 #include "AliAODEvent.h"
39 #include "AliOADBContainer.h"
40 #include "AliGenCocktailEventHeader.h"
41 #include "AliGenEposEventHeader.h"
42 #include "AliGenHijingEventHeader.h"
43 #include "AliGenGeVSimEventHeader.h"
44 #include "AliCollisionGeometry.h"
45 #include "AliMultiplicity.h"
46 #include "AliFlowTrackCuts.h"
47 #include "AliFlowEventSimple.h"
48 #include "AliFlowTrack.h"
49 #include "AliFlowVector.h"
50 #include "AliFlowEvent.h"
51 #include "AliLog.h"
52 
53 using std::endl;
54 using std::cout;
55 ClassImp(AliFlowEvent)
56 
57 //-----------------------------------------------------------------------
58 AliFlowEvent::AliFlowEvent():
59  AliFlowEventSimple(), fApplyRecentering(-1), fDivSigma(kTRUE), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0), fChi2A(0x0), fChi2C(0x0), fChi3A(0x0), fChi3C(0x0)
60 {
61  // constructor
62  for(Int_t i(0); i < 9; i++) {
63  for(Int_t j(0); j < 2; j++) {
64  for(Int_t k(0); k < 2; k++) {
65  fMeanQ[i][j][k] = 0.;
66  fWidthQ[i][j][k] = 0.;
67  fMeanQv3[i][j][k] = 0.;
68  fWidthQv3[i][j][k] = 0.;
69  }
70  }
71  }
72  for(Int_t i(0); i < 5; i++) {
73  fQxavsV0[i] = 0x0;
74  fQyavsV0[i] = 0x0;
75  fQxcvsV0[i] = 0x0;
76  fQycvsV0[i] = 0x0;
77  }
78  //ctor
79  cout << "AliFlowEvent: Default constructor to be used only by root for io" << endl;
80 }
81 
82 //-----------------------------------------------------------------------
83 AliFlowEvent::AliFlowEvent(Int_t n):
84  AliFlowEventSimple(n), fApplyRecentering(-1), fDivSigma(kTRUE), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0), fChi2A(0x0), fChi2C(0x0), fChi3A(0x0), fChi3C(0x0)
85 {
86  // constructor
87  for(Int_t i(0); i < 9; i++) {
88  for(Int_t j(0); j < 2; j++) {
89  for(Int_t k(0); k < 2; k++) {
90  fMeanQ[i][j][k] = 0.;
91  fWidthQ[i][j][k] = 0.;
92  fMeanQv3[i][j][k] = 0.;
93  fWidthQv3[i][j][k] = 0.;
94  }
95  }
96  }
97  for(Int_t i(0); i < 5; i++) {
98  fQxavsV0[i] = 0x0;
99  fQyavsV0[i] = 0x0;
100  fQxcvsV0[i] = 0x0;
101  fQycvsV0[i] = 0x0;
102  }
103 }
104 
105 //-----------------------------------------------------------------------
106 AliFlowEvent::AliFlowEvent(const AliFlowEvent& event):
107  AliFlowEventSimple(event), fApplyRecentering(event.fApplyRecentering), fDivSigma(event.fDivSigma), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0), fChi2A(0x0), fChi2C(0x0), fChi3A(0x0), fChi3C(0x0)
108 {
109  // copy constructor
110  for(Int_t i(0); i < 9; i++) {
111  for(Int_t j(0); j < 2; j++) {
112  for(Int_t k(0); k < 2; k++) {
113  fMeanQ[i][j][k] = 0.;
114  fWidthQ[i][j][k] = 0.;
115  fMeanQv3[i][j][k] = 0.;
116  fWidthQv3[i][j][k] = 0.;
117  }
118  }
119  }
120  for(Int_t i(0); i < 5; i++) {
121  fQxavsV0[i] = 0x0;
122  fQyavsV0[i] = 0x0;
123  fQxcvsV0[i] = 0x0;
124  fQycvsV0[i] = 0x0;
125  }
126 }
127 
128 //-----------------------------------------------------------------------
129 AliFlowEvent& AliFlowEvent::operator=(const AliFlowEvent& event)
130 {
131  //assignment operator
132  if (&event==this) return *this; // check self-assignment
133 
134  fApplyRecentering = event.fApplyRecentering;
135  fCachedRun = event.fCachedRun;
136  fVZEROcentralityBin = event.fVZEROcentralityBin;
137  fDivSigma = event.fDivSigma;
138  fEvent = 0x0; // should never be copied
139  fChi2A = 0x0; // do not clone these; if 0x0 they will be retrieved from the rp cuts object
140  fChi2C = 0x0;
141  fChi3A = 0x0;
142  fChi3C = 0x0;
143  for(Int_t i(0); i < 9; i++) {
144  for(Int_t j(0); j < 2; j++) {
145  for(Int_t k(0); k < 2; k++) {
146  fMeanQ[i][j][k] = event.fMeanQ[i][j][k];
147  fWidthQ[i][j][k] = event.fWidthQ[i][j][k];
148  fMeanQv3[i][j][k] = event.fMeanQv3[i][j][k];
149  fWidthQv3[i][j][k] = event.fWidthQv3[i][j][k];
150  }
151  }
152  }
153  for(Int_t i(0); i < 5; i++) {
154  fQxavsV0[i] = 0x0;
155  fQyavsV0[i] = 0x0;
156  fQxcvsV0[i] = 0x0;
157  fQycvsV0[i] = 0x0;
158  }
159 
160  AliFlowEventSimple::operator=(event);
161  return *this;
162 }
163 //-----------------------------------------------------------------------
164 AliFlowTrack* AliFlowEvent::GetTrack(Int_t i)
165 {
166  //get track i from collection
167  if (i>=fNumberOfTracks) return NULL;
168  AliFlowTrack* pTrack = static_cast<AliFlowTrack*>(fTrackCollection->At(i)) ;
169  return pTrack;
170 }
171 
172 //-----------------------------------------------------------------------
173 void AliFlowEvent::SetMCReactionPlaneAngle(const AliMCEvent* mcEvent)
174 {
175  //sets the event plane angle from the proper header in the MC
176 
177  //COCKTAIL with HIJING
178  if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"Cocktail Header")) //returns 0 if matches
179  {
180  AliGenCocktailEventHeader *headerC = dynamic_cast<AliGenCocktailEventHeader *> (mcEvent-> GenEventHeader());
181  if (headerC)
182  {
183  TList *lhd = headerC->GetHeaders();
184  if (lhd)
185  {
186  AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (lhd->At(0));
187  if (hdh) AliFlowEventSimple::SetMCReactionPlaneAngle( hdh->ReactionPlaneAngle() );
188  }
189  }
190  }
191  //THERMINATOR
192  else if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"Therminator")) //returns 0 if matches
193  {
194  AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
195  if (headerH) AliFlowEventSimple::SetMCReactionPlaneAngle( headerH->ReactionPlaneAngle() );
196  }
197  //GEVSIM
198  else if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"GeVSim header")) //returns 0 if matches
199  {
200  AliGenGeVSimEventHeader* headerG = dynamic_cast<AliGenGeVSimEventHeader*>(mcEvent->GenEventHeader());
201  if (headerG) AliFlowEventSimple::SetMCReactionPlaneAngle( headerG->GetEventPlane() );
202  }
203  //HIJING
204  else if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"Hijing")) //returns 0 if matches
205  {
206  AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
207  if (headerH) AliFlowEventSimple::SetMCReactionPlaneAngle( headerH->ReactionPlaneAngle() );
208  }
209  //AMPT
210  else if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"Ampt")) //returns 0 if matches
211  {
212  AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
213  if (headerH) AliFlowEventSimple::SetMCReactionPlaneAngle( headerH->ReactionPlaneAngle() );
214  }
215  //EPOS
216  else if (!strcmp(mcEvent->GenEventHeader()->GetName(),"EPOS"))
217  {
218  AliGenEposEventHeader* headerE = dynamic_cast<AliGenEposEventHeader*>(mcEvent->GenEventHeader());
219  if (headerE) AliFlowEventSimple::SetMCReactionPlaneAngle( headerE->ReactionPlaneAngle() );
220  }
221  //Hydjet
222  else
223  {
224  AliCollisionGeometry* header = dynamic_cast<AliCollisionGeometry*>(mcEvent->GenEventHeader());
225  if (header) AliFlowEventSimple::SetMCReactionPlaneAngle( header->ReactionPlaneAngle() );
226  }
227 }
228 
229 //-----------------------------------------------------------------------
230 AliFlowEvent::AliFlowEvent( const AliMCEvent* anInput,
231  const AliCFManager* rpCFManager,
232  const AliCFManager* poiCFManager):
233  AliFlowEventSimple(20), fApplyRecentering(-1), fDivSigma(kTRUE), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0), fChi2A(0x0), fChi2C(0x0), fChi3A(0x0), fChi3C(0x0)
234 {
235  // constructor
236  for(Int_t i(0); i < 9; i++) {
237  for(Int_t j(0); j < 2; j++) {
238  for(Int_t k(0); k < 2; k++) {
239  fMeanQ[i][j][k] = 0.;
240  fWidthQ[i][j][k] = 0.;
241  fMeanQv3[i][j][k] = 0.;
242  fWidthQv3[i][j][k] = 0.;
243  }
244  }
245  }
246  for(Int_t i(0); i < 5; i++) {
247  fQxavsV0[i] = 0x0;
248  fQyavsV0[i] = 0x0;
249  fQxcvsV0[i] = 0x0;
250  fQycvsV0[i] = 0x0;
251  }
252 
253  //Fills the event from the MC kinematic information
254  Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
255 
256  //loop over tracks
257  for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
258  {
259  //get input particle
260  AliMCParticle* pParticle = dynamic_cast<AliMCParticle*>(anInput->GetTrack(itrkN));
261  if (!pParticle) continue;
262 
263  //check if pParticle passes the cuts
264  Bool_t rpOK = kTRUE;
265  Bool_t poiOK = kTRUE;
266  if (rpCFManager && poiCFManager)
267  {
268  rpOK = rpCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pParticle);
269  poiOK = poiCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pParticle);
270  }
271  if (!(rpOK||poiOK)) continue;
272 
273  AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
274  pTrack->SetSource(AliFlowTrack::kFromMC);
275 
276  if (rpOK && rpCFManager)
277  {
278  pTrack->SetForRPSelection(kTRUE);
279  IncrementNumberOfPOIs(0);
280  }
281  if (poiOK && poiCFManager)
282  {
283  pTrack->SetForPOISelection(kTRUE);
284  IncrementNumberOfPOIs(1);
285  }
286 
287  AddTrack(pTrack) ;
288  }//for all tracks
289  SetMCReactionPlaneAngle(anInput);
290 }
291 
292 //-----------------------------------------------------------------------
293 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
294  const AliCFManager* rpCFManager,
295  const AliCFManager* poiCFManager ):
296  AliFlowEventSimple(20), fApplyRecentering(-1), fDivSigma(kTRUE), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0), fChi2A(0x0), fChi2C(0x0), fChi3A(0x0), fChi3C(0x0)
297 {
298  // constructor
299  for(Int_t i(0); i < 9; i++) {
300  for(Int_t j(0); j < 2; j++) {
301  for(Int_t k(0); k < 2; k++) {
302  fMeanQ[i][j][k] = 0.;
303  fWidthQ[i][j][k] = 0.;
304  fMeanQv3[i][j][k] = 0.;
305  fWidthQv3[i][j][k] = 0.;
306  }
307  }
308  }
309  for(Int_t i(0); i < 5; i++) {
310  fQxavsV0[i] = 0x0;
311  fQyavsV0[i] = 0x0;
312  fQxcvsV0[i] = 0x0;
313  fQycvsV0[i] = 0x0;
314  }
315  //set run number
316  if(anInput->GetRunNumber()) fRun = anInput->GetRunNumber();
317 
318  //Fills the event from the ESD
319 
320  Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
321 
322  //loop over tracks
323  for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
324  {
325  AliESDtrack* pParticle = anInput->GetTrack(itrkN); //get input particle
326 
327  //check if pParticle passes the cuts
328  Bool_t rpOK = kTRUE;
329  Bool_t poiOK = kTRUE;
330  if (rpCFManager && poiCFManager)
331  {
332  rpOK = ( rpCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
333  rpCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
334  poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
335  poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
336  }
337  if (!(rpOK || poiOK)) continue;
338 
339  //make new AliFLowTrack
340  AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
341  pTrack->SetSource(AliFlowTrack::kFromESD);
342 
343  //marking the particles used for int. flow:
344  if(rpOK && rpCFManager)
345  {
346  pTrack->SetForRPSelection(kTRUE);
347  IncrementNumberOfPOIs(0);
348  }
349  //marking the particles used for diff. flow:
350  if(poiOK && poiCFManager)
351  {
352  pTrack->SetForPOISelection(kTRUE);
353  IncrementNumberOfPOIs(1);
354  }
355 
356  AddTrack(pTrack);
357  }//end of while (itrkN < iNumberOfInputTracks)
358 }
359 
360 //-----------------------------------------------------------------------
361 AliFlowEvent::AliFlowEvent( const AliAODEvent* anInput,
362  const AliCFManager* rpCFManager,
363  const AliCFManager* poiCFManager):
364  AliFlowEventSimple(20), fApplyRecentering(-1), fDivSigma(kTRUE), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0), fChi2A(0x0), fChi2C(0x0), fChi3A(0x0), fChi3C(0x0)
365 {
366  // constructor
367  for(Int_t i(0); i < 9; i++) {
368  for(Int_t j(0); j < 2; j++) {
369  for(Int_t k(0); k < 2; k++) {
370  fMeanQ[i][j][k] = 0.;
371  fWidthQ[i][j][k] = 0.;
372  fMeanQv3[i][j][k] = 0.;
373  fWidthQv3[i][j][k] = 0.;
374  }
375  }
376  }
377  for(Int_t i(0); i < 5; i++) {
378  fQxavsV0[i] = 0x0;
379  fQyavsV0[i] = 0x0;
380  fQxcvsV0[i] = 0x0;
381  fQycvsV0[i] = 0x0;
382  }
383 
384  //set run number
385  if(anInput->GetRunNumber()) fRun = anInput->GetRunNumber();
386 
387  //Fills the event from the AOD
388  Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
389 
390  //loop over tracks
391  for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
392  {
393  AliAODTrack* pParticle = dynamic_cast<AliAODTrack*>(anInput->GetTrack(itrkN));
394  if(!pParticle) AliFatal("Not a standard AOD"); //get input particle
395 
396  //check if pParticle passes the cuts
397  Bool_t rpOK = kTRUE;
398  Bool_t poiOK = kTRUE;
399  if (rpCFManager && poiCFManager)
400  {
401  rpOK = ( rpCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
402  rpCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
403  poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
404  poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
405  }
406  if (!(rpOK || poiOK)) continue;
407 
408  //make new AliFlowTrack
409  AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
410  pTrack->SetSource(AliFlowTrack::kFromAOD);
411 
412  if (rpOK /* && rpCFManager */ ) // to be fixed - with CF managers uncommented only empty events (NULL in header files)
413  {
414  pTrack->SetForRPSelection(kTRUE);
415  IncrementNumberOfPOIs(0);
416  }
417  if (poiOK /* && poiCFManager*/ )
418  {
419  pTrack->SetForPOISelection(kTRUE);
420  IncrementNumberOfPOIs(1);
421  }
422  AddTrack(pTrack);
423  }
424 
425  // if (iSelParticlesRP >= fMinMult && iSelParticlesRP <= fMaxMult)
426  // {
427  // if ( (++fCount % 100) == 0)
428  // {
429  // if (!fMCReactionPlaneAngle == 0) cout<<" MC Reaction Plane Angle = "<< fMCReactionPlaneAngle << endl;
430  // else cout<<" MC Reaction Plane Angle = unknown "<< endl;
431  // cout<<" iGoodTracks = "<<iGoodTracks<<endl;
432  // cout<<" # of RP selected tracks = "<<iSelParticlesRP<<endl;
433  // cout<<" # of POI selected tracks = "<<iSelParticlesPOI<<endl;
434  // cout << "# " << fCount << " events processed" << endl;
435  // }
436  // return pEvent;
437  // }
438  // else
439  // {
440  // cout<<"Not enough tracks in the FlowEventSimple"<<endl;
441  // return 0;
442  // }
443  //}
444  //else
445  //{
446  // cout<<"Event does not pass multiplicity cuts"<<endl;
447  // return 0;
448  //}
449 
450 }
451 
452 //-----------------------------------------------------------------------
453 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
454  const AliMCEvent* anInputMc,
455  KineSource anOption,
456  const AliCFManager* rpCFManager,
457  const AliCFManager* poiCFManager ):
458  AliFlowEventSimple(20), fApplyRecentering(-1), fDivSigma(kTRUE), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0), fChi2A(0x0), fChi2C(0x0), fChi3A(0x0), fChi3C(0x0)
459 {
460  // constructor
461  for(Int_t i(0); i < 9; i++) {
462  for(Int_t j(0); j < 2; j++) {
463  for(Int_t k(0); k < 2; k++) {
464  fMeanQ[i][j][k] = 0.;
465  fWidthQ[i][j][k] = 0.;
466  fMeanQv3[i][j][k] = 0.;
467  fWidthQv3[i][j][k] = 0.;
468  }
469  }
470  }
471  for(Int_t i(0); i < 5; i++) {
472  fQxavsV0[i] = 0x0;
473  fQyavsV0[i] = 0x0;
474  fQxcvsV0[i] = 0x0;
475  fQycvsV0[i] = 0x0;
476  }
477  //fills the event with tracks from the ESD and kinematics from the MC info via the track label
478  if (anOption==kNoKine)
479  {
480  AliFatal("WRONG OPTION IN AliFlowEventMaker::FillTracks(AliESDEvent* anInput, AliMCEvent* anInputMc, KineSource anOption)");
481  exit(1);
482  }
483 
484  Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
485 
486  Int_t iNumberOfInputTracksMC = anInputMc->GetNumberOfTracks() ;
487  if (iNumberOfInputTracksMC==-1)
488  {
489  AliError("Skipping Event -- No MC information available for this event");
490  return;
491  }
492 
493  //loop over ESD tracks
494  for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
495  {
496  AliESDtrack* pParticle = anInput->GetTrack(itrkN); //get input particle
497  //get Label
498  Int_t iLabel = pParticle->GetLabel();
499  //match to mc particle
500  AliMCParticle* pMcParticle = (AliMCParticle*) anInputMc->GetTrack(TMath::Abs(iLabel));
501 
502  //check
503  if (TMath::Abs(pParticle->GetLabel())!=pMcParticle->Label())
504  AliWarning(Form("pParticle->GetLabel()!=pMcParticle->Label(), %i, %i", pParticle->GetLabel(), pMcParticle->Label()));
505 
506  //check if pParticle passes the cuts
507  Bool_t rpOK = kFALSE;
508  Bool_t poiOK = kFALSE;
509  if (rpCFManager && poiCFManager)
510  {
511  if(anOption == kESDkine)
512  {
513  if (rpCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pMcParticle,"mcGenCuts1") &&
514  rpCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle))
515  rpOK=kTRUE;
516  if (poiCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pMcParticle,"mcGenCuts2") &&
517  poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle))
518  poiOK=kTRUE;
519  }
520  else if (anOption == kMCkine)
521  {
522  if (rpCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pMcParticle))
523  rpOK=kTRUE;
524  if (poiCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pMcParticle))
525  poiOK=kTRUE;
526  }
527  }
528 
529  if (!(rpOK || poiOK)) continue;
530 
531  //make new AliFlowTrack
532  AliFlowTrack* pTrack = NULL;
533  if(anOption == kESDkine) //take the PID from the MC & the kinematics from the ESD
534  {
535  pTrack = new AliFlowTrack(pParticle);
536  }
537  else if (anOption == kMCkine) //take the PID and kinematics from the MC
538  {
539  pTrack = new AliFlowTrack(pMcParticle);
540  }
541 
542  if (rpOK && rpCFManager)
543  {
544  IncrementNumberOfPOIs(0);
545  pTrack->SetForRPSelection();
546  }
547  if (poiOK && poiCFManager)
548  {
549  IncrementNumberOfPOIs(1);
550  pTrack->SetForPOISelection();
551  }
552 
553  AddTrack(pTrack);
554  }
555  SetMCReactionPlaneAngle(anInputMc);
556 }
557 
558 //-----------------------------------------------------------------------
559 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
560  const AliMultiplicity* anInputTracklets,
561  const AliCFManager* poiCFManager ):
562  AliFlowEventSimple(20), fApplyRecentering(-1), fDivSigma(kTRUE), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0), fChi2A(0x0), fChi2C(0x0), fChi3A(0x0), fChi3C(0x0)
563 {
564  // constructor
565  for(Int_t i(0); i < 9; i++) {
566  for(Int_t j(0); j < 2; j++) {
567  for(Int_t k(0); k < 2; k++) {
568  fMeanQ[i][j][k] = 0.;
569  fWidthQ[i][j][k] = 0.;
570  fMeanQv3[i][j][k] = 0.;
571  fWidthQv3[i][j][k] = 0.;
572  }
573  }
574  }
575  for(Int_t i(0); i < 5; i++) {
576  fQxavsV0[i] = 0x0;
577  fQyavsV0[i] = 0x0;
578  fQxcvsV0[i] = 0x0;
579  fQycvsV0[i] = 0x0;
580  }
581 
582  //Select the particles of interest from the ESD
583  Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
584 
585  //loop over tracks
586  for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
587  {
588  AliESDtrack* pParticle = anInput->GetTrack(itrkN); //get input particle
589 
590  //check if pParticle passes the cuts
591  Bool_t poiOK = kTRUE;
592  if (poiCFManager)
593  {
594  poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
595  poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
596  }
597  if (!poiOK) continue;
598 
599  //make new AliFLowTrack
600  AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
601 
602  //marking the particles used for the particle of interest (POI) selection:
603  if(poiOK && poiCFManager)
604  {
605  IncrementNumberOfPOIs(1);
606  pTrack->SetForPOISelection(kTRUE);
607  pTrack->SetSource(AliFlowTrack::kFromESD);
608  }
609 
610  AddTrack(pTrack);
611  }//end of while (itrkN < iNumberOfInputTracks)
612 
613  //Select the reference particles from the SPD tracklets
614  anInputTracklets = anInput->GetMultiplicity();
615  Int_t multSPD = anInputTracklets->GetNumberOfTracklets();
616 
617  //loop over tracklets
618  for (Int_t itracklet=0; itracklet<multSPD; ++itracklet) {
619  Float_t thetaTr= anInputTracklets->GetTheta(itracklet);
620  Float_t phiTr= anInputTracklets->GetPhi(itracklet);
621  // calculate eta
622  Float_t etaTr = -TMath::Log(TMath::Tan(thetaTr/2.));
623 
624  //make new AliFLowTrackSimple
625  AliFlowTrack* pTrack = new AliFlowTrack();
626  pTrack->SetPt(0.0);
627  pTrack->SetEta(etaTr);
628  pTrack->SetPhi(phiTr);
629  //marking the particles used for the reference particle (RP) selection:
630  IncrementNumberOfPOIs(0);
631  pTrack->SetForRPSelection(kTRUE);
632  pTrack->SetSource(AliFlowTrack::kFromTracklet);
633 
634  //Add the track to the flowevent
635  AddTrack(pTrack);
636  }
637 
638 }
639 
640 //-----------------------------------------------------------------------
641 AliFlowEvent::AliFlowEvent( const AliESDEvent* esd,
642  const AliCFManager* poiCFManager,
643  Bool_t hybrid):
644  AliFlowEventSimple(20), fApplyRecentering(-1), fDivSigma(kTRUE), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0), fChi2A(0x0), fChi2C(0x0), fChi3A(0x0), fChi3C(0x0)
645 {
646  // constructor
647  for(Int_t i(0); i < 9; i++) {
648  for(Int_t j(0); j < 2; j++) {
649  for(Int_t k(0); k < 2; k++) {
650  fMeanQ[i][j][k] = 0.;
651  fWidthQ[i][j][k] = 0.;
652  fMeanQv3[i][j][k] = 0.;
653  fWidthQv3[i][j][k] = 0.;
654  }
655  }
656  }
657  for(Int_t i(0); i < 5; i++) {
658  fQxavsV0[i] = 0x0;
659  fQyavsV0[i] = 0x0;
660  fQxcvsV0[i] = 0x0;
661  fQycvsV0[i] = 0x0;
662  }
663  //Select the particles of interest from the ESD
664  Int_t iNumberOfInputTracks = esd->GetNumberOfTracks() ;
665 
666  //Double_t gPt = 0.0, gP = 0.0;
667  Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0}; //The impact parameters and their covariance.
668 // Double_t dca3D = 0.0; FIXME unused variable
669 
670  AliESDtrack trackTPC;
671 
672  //loop over tracks
673  for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
674  {
675 
676  if (!esd->GetTrack(itrkN)) continue;
677 
678  Bool_t useTPC = kFALSE;
679 
680  AliESDtrack* pParticle = esd->GetTrack(itrkN); //get input particle
681 
682  //check if pParticle passes the cuts
683  Bool_t poiOK = kTRUE;
684 
685  if (poiCFManager)
686  {
687  poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
688  poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
689  }
690 
691  if (!(poiOK)) continue;
692 
693  AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)pParticle->GetTPCInnerParam();
694 
695  if (tpcTrack)
696  {
697 
698 // gPt = tpcTrack->Pt();
699 // gP = tpcTrack->P();
700 
701  useTPC = kTRUE;
702 
703  const AliESDVertex *vertexSPD = esd->GetPrimaryVertexSPD();
704  const AliESDVertex *vertexTPC = esd->GetPrimaryVertexTPC();
705 
706  AliExternalTrackParam copy(*tpcTrack);
707  if(hybrid)
708  copy.PropagateToDCA(vertexSPD,esd->GetMagneticField(),100.,dca,cov);
709  else
710  copy.PropagateToDCA(vertexTPC,esd->GetMagneticField(),100.,dca,cov);
711 
712 // dca3D = TMath::Sqrt(TMath::Power(dca[0],2)+TMath::Power(dca[1],2)); FIXME unused variable
713 
714  }
715 
716  //make new AliFLowTrack
717  AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
718 
719  pTrack->SetSource(AliFlowTrack::kFromESD);
720 
721  //marking the particles used for diff. flow:
722  if(poiOK && poiCFManager)
723  {
724  pTrack->SetForPOISelection(kTRUE);
725  IncrementNumberOfPOIs(1);
726  }
727 
728  if(useTPC)
729  {
730  pTrack->SetForRPSelection(kTRUE);
731  IncrementNumberOfPOIs(0);
732  }
733 
734  AddTrack(pTrack);
735 
736  }//end of while (itrkN < iNumberOfInputTracks)
737 
738 }
739 
740 //-----------------------------------------------------------------------
741 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
742  const TH2F* anInputFMDhist,
743  const AliCFManager* poiCFManager ):
744  AliFlowEventSimple(20), fApplyRecentering(-1), fDivSigma(kTRUE), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0), fChi2A(0x0), fChi2C(0x0), fChi3A(0x0), fChi3C(0x0)
745 {
746  // constructor
747  for(Int_t i(0); i < 9; i++) {
748  for(Int_t j(0); j < 2; j++) {
749  for(Int_t k(0); k < 2; k++) {
750  fMeanQ[i][j][k] = 0.;
751  fWidthQ[i][j][k] = 0.;
752  fMeanQv3[i][j][k] = 0.;
753  fWidthQv3[i][j][k] = 0.;
754  }
755  }
756  }
757  for(Int_t i(0); i < 5; i++) {
758  fQxavsV0[i] = 0x0;
759  fQyavsV0[i] = 0x0;
760  fQxcvsV0[i] = 0x0;
761  fQycvsV0[i] = 0x0;
762  }
763 
764  //Select the particles of interest from the ESD
765  Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
766 
767  //loop over tracks
768  for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
769  {
770  AliESDtrack* pParticle = anInput->GetTrack(itrkN); //get input particle
771 
772  //check if pParticle passes the cuts
773  Bool_t poiOK = kTRUE;
774  if (poiCFManager)
775  {
776  poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
777  poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
778  }
779  if (!poiOK) continue;
780 
781  //make new AliFLowTrack
782  AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
783 
784  //marking the particles used for the particle of interest (POI) selection:
785  if(poiOK && poiCFManager)
786  {
787  IncrementNumberOfPOIs(1);
788  pTrack->SetForPOISelection(kTRUE);
789  pTrack->SetSource(AliFlowTrack::kFromESD);
790  }
791 
792  AddTrack(pTrack);
793  }//end of while (itrkN < iNumberOfInputTracks)
794 
795  //Select the reference particles from the FMD hits
796  //loop over FMD histogram
797  Int_t iBinsEta = anInputFMDhist->GetNbinsX();
798  Int_t iBinsPhi = anInputFMDhist->GetNbinsY();
799 
800  for (Int_t iEta = 1; iEta <= iBinsEta; iEta++){
801  Double_t etaFMD = anInputFMDhist->GetXaxis()->GetBinCenter(iEta);
802  for (Int_t iPhi = 1; iPhi <= iBinsPhi; iPhi++){
803  Double_t phiFMD = anInputFMDhist->GetYaxis()->GetBinCenter(iPhi);
804  Double_t weightFMD = anInputFMDhist->GetBinContent(iEta,iPhi);
805 
806  if (weightFMD > 0.0) { //do not add empty bins
807  //make new AliFLowTrackSimple
808  AliFlowTrack* pTrack = new AliFlowTrack();
809  pTrack->SetPt(0.0);
810  pTrack->SetEta(etaFMD);
811  pTrack->SetPhi(phiFMD);
812  pTrack->SetWeight(weightFMD);
813  //marking the particles used for the reference particle (RP) selection:
814  pTrack->TagRP();
815  IncrementNumberOfPOIs(0);
816  pTrack->SetSource(AliFlowTrack::kFromFMD);
817 
818  //Add the track to the flowevent
819  AddTrack(pTrack);
820 
821  }
822  }
823  }
824 
825 }
826 
827 //-----------------------------------------------------------------------
828 void AliFlowEvent::FindDaughters(Bool_t keepDaughtersInRPselection)
829 {
830  //each flow track holds it's esd track index as well as its daughters esd index.
831  //fill the array of daughters for every track with the pointers to flow tracks
832  //to associate the mothers with daughters directly
833  for (Int_t iTrack=0; iTrack<fMothersCollection->GetEntriesFast(); iTrack++)
834  {
835  AliFlowTrack* mother = static_cast<AliFlowTrack*>(fMothersCollection->At(iTrack));
836  if (!mother) continue;
837  if (mother->GetNDaughters()<1) continue;
838  for (Int_t iDaughterCandidate=0; iDaughterCandidate<fNumberOfTracks; iDaughterCandidate++)
839  {
840  AliFlowTrack* daughterCandidate = static_cast<AliFlowTrack*>(fTrackCollection->At(iDaughterCandidate));
841  Int_t esdIndexDaughterCandidate = daughterCandidate->GetID();
842  for (Int_t iDaughter=0; iDaughter<mother->GetNDaughters(); iDaughter++)
843  {
844  Int_t esdIndexDaughter = mother->GetIDDaughter(iDaughter);
845  if (esdIndexDaughter==esdIndexDaughterCandidate)
846  {
847  mother->SetDaughter(iDaughter,daughterCandidate);
848  daughterCandidate->SetForRPSelection(keepDaughtersInRPselection);
849  }
850  }
851  }
852  }
853 }
854 
855 //-----------------------------------------------------------------------
856 void AliFlowEvent::Fill( AliFlowTrackCuts* rpCuts,
857  AliFlowTrackCuts* poiCuts )
858 {
859  //Fills the event from a vevent: AliESDEvent,AliAODEvent,AliMCEvent
860  //the input data needs to be attached to the cuts
861  //we have two cases, if we're cutting the same collection of tracks
862  //(same param type) then we can have tracks that are both rp and poi
863  //in the other case we want to have two exclusive sets of rps and pois
864  //e.g. one tracklets, the other PMD or global - USER IS RESPOSIBLE
865  //FOR MAKING SURE THEY DONT OVERLAP OR ELSE THE SAME PARTICLE WILL BE
866  //TAKEN TWICE
867 
868  ClearFast();
869  if (!rpCuts || !poiCuts) return;
870  AliFlowTrackCuts::trackParameterType sourceRP = rpCuts->GetParamType();
871  AliFlowTrackCuts::trackParameterType sourcePOI = poiCuts->GetParamType();
872  AliFlowTrack* pTrack=NULL;
873 
874  //set run
875  if(rpCuts->GetRun()) fRun = rpCuts->GetRun();
876 
877  // if the source for rp's or poi's is the VZERO detector, get the calibration
878  // and set the calibration parameters
879  if (sourceRP == AliFlowTrackCuts::kBetaVZERO) {
880  SetBetaVZEROCalibrationForTrackCuts(rpCuts);
881  fDivSigma = rpCuts->GetDivSigma();
882  if(!rpCuts->GetApplyRecentering()) {
883  // if the user does not want to recenter, switch the flag
884  fApplyRecentering = -1;
885  }
886  // note: this flag is used in the overloaded implementation of Get2Qsub()
887  // and tells the function to use as Qsub vectors the recentered Q-vectors
888  // from the VZERO oadb file or from the event header
889  }
890  if (sourcePOI == AliFlowTrackCuts::kBetaVZERO) {
891  // probably no-one will choose vzero tracks as poi's ...
892  SetBetaVZEROCalibrationForTrackCuts(poiCuts);
893  fDivSigma = poiCuts->GetDivSigma();
894  }
895 
896  if (sourceRP == AliFlowTrackCuts::kDeltaVZERO) {
897  SetDeltaVZEROCalibrationForTrackCuts(rpCuts);
898  fDivSigma = rpCuts->GetDivSigma();
899  if(!rpCuts->GetApplyRecentering()) {
900  // if the user does not want to recenter, switch the flag
901  fApplyRecentering = -1;
902  }
903  }
904  if (sourcePOI == AliFlowTrackCuts::kDeltaVZERO) {
905  // probably no-one will choose vzero tracks as poi's ...
906  SetDeltaVZEROCalibrationForTrackCuts(poiCuts);
907  fDivSigma = poiCuts->GetDivSigma();
908  }
909 
910  if (sourceRP == AliFlowTrackCuts::kVZERO) {
911  SetVZEROCalibrationForTrackCuts(rpCuts);
912  if(!rpCuts->GetApplyRecentering()) {
913  // if the user does not want to recenter, switch the flag
914  fApplyRecentering = -1;
915  }
916  // note: this flag is used in the overloaded implementation of Get2Qsub()
917  // and tells the function to use as Qsub vectors the recentered Q-vectors
918  // from the VZERO oadb file or from the event header
919  }
920  if (sourcePOI == AliFlowTrackCuts::kVZERO) {
921  // probably no-one will choose vzero tracks as poi's ...
922  SetVZEROCalibrationForTrackCuts(poiCuts);
923  }
924 
925  if (sourceRP==sourcePOI)
926  {
927  //loop over tracks
928  Int_t numberOfInputObjects = rpCuts->GetNumberOfInputObjects();
929  for (Int_t i=0; i<numberOfInputObjects; i++)
930  {
931  //get input object (particle)
932  TObject* particle = rpCuts->GetInputObject(i);
933 
934  Bool_t rp = rpCuts->IsSelected(particle,i);
935  Bool_t poi = poiCuts->IsSelected(particle,i);
936 
937  if (!(rp||poi)) continue;
938 
939  //make new AliFlowTrack
940  if (rp)
941  {
942  pTrack = rpCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
943  if (!pTrack) continue;
944  pTrack->Tag(0); IncrementNumberOfPOIs(0);
945  if (poi) {pTrack->Tag(1); IncrementNumberOfPOIs(1);}
946  if (pTrack->GetNDaughters()>0) fMothersCollection->Add(pTrack);
947  }
948  else if (poi)
949  {
950  pTrack = poiCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
951  if (!pTrack) continue;
952  pTrack->Tag(1); IncrementNumberOfPOIs(1);
953  if (pTrack->GetNDaughters()>0) fMothersCollection->Add(pTrack);
954  }
955  fNumberOfTracks++;
956  }//end of while (i < numberOfTracks)
957  }
958  else if (sourceRP!=sourcePOI)
959  {
960  //here we have two different sources of particles, so we fill
961  //them independently
962  //POI
963  for (Int_t i=0; i<poiCuts->GetNumberOfInputObjects(); i++)
964  {
965  TObject* particle = poiCuts->GetInputObject(i);
966  Bool_t poi = poiCuts->IsSelected(particle,i);
967  if (!poi) continue;
968  pTrack = poiCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
969  if (!pTrack) continue;
970  pTrack->Tag(1);
971  IncrementNumberOfPOIs(1);
972  fNumberOfTracks++;
973  if (pTrack->GetNDaughters()>0) fMothersCollection->Add(pTrack);
974  }
975  //RP
976  Int_t numberOfInputObjects = rpCuts->GetNumberOfInputObjects();
977  for (Int_t i=0; i<numberOfInputObjects; i++)
978  {
979  TObject* particle = rpCuts->GetInputObject(i);
980  Bool_t rp = rpCuts->IsSelected(particle,i);
981  if (!rp) continue;
982  pTrack = rpCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
983  if (!pTrack) continue;
984  pTrack->Tag(0);
985  IncrementNumberOfPOIs(0);
986  fNumberOfTracks++;
987  if (pTrack->GetNDaughters()>0) fMothersCollection->Add(pTrack);
988  }
989  }
990 }
991 
992 //-----------------------------------------------------------------------
993 void AliFlowEvent::InsertTrack(AliFlowTrack *track) {
994  // adds a flow track at the end of the container
995  AliFlowTrack *pTrack = ReuseTrack( fNumberOfTracks++ );
996  *pTrack = *track;
997  if (track->GetNDaughters()>0)
998  {
999  fMothersCollection->Add(pTrack);
1000  }
1001  return;
1002 }
1003 
1004 //-----------------------------------------------------------------------
1005 AliFlowTrack* AliFlowEvent::ReuseTrack(Int_t i)
1006 {
1007  //try to reuse an existing track, if empty, make new one
1008  AliFlowTrack* pTrack = static_cast<AliFlowTrack*>(fTrackCollection->At(i));
1009  if (pTrack)
1010  {
1011  pTrack->Clear();
1012  }
1013  else
1014  {
1015  pTrack = new AliFlowTrack();
1016  fTrackCollection->AddAtAndExpand(pTrack,i);
1017  }
1018  return pTrack;
1019 }
1020 
1021 //-----------------------------------------------------------------------
1022 AliFlowEvent::AliFlowEvent( AliFlowTrackCuts* rpCuts,
1023  AliFlowTrackCuts* poiCuts ):
1024  AliFlowEventSimple(20), fApplyRecentering(kFALSE), fDivSigma(kTRUE), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0), fChi2A(0x0), fChi2C(0x0), fChi3A(0x0), fChi3C(0x0)
1025 
1026 {
1027  // constructor
1028  for(Int_t i(0); i < 9; i++) {
1029  for(Int_t j(0); j < 2; j++) {
1030  for(Int_t k(0); k < 2; k++) {
1031  fMeanQ[i][j][k] = 0.;
1032  fWidthQ[i][j][k] = 0.;
1033  fMeanQv3[i][j][k] = 0.;
1034  fWidthQv3[i][j][k] = 0.;
1035  }
1036  }
1037  }
1038  for(Int_t i(0); i < 5; i++) {
1039  fQxavsV0[i] = 0x0;
1040  fQyavsV0[i] = 0x0;
1041  fQxcvsV0[i] = 0x0;
1042  fQycvsV0[i] = 0x0;
1043  }
1044  //Fills the event from a vevent: AliESDEvent,AliAODEvent,AliMCEvent
1045  //the input data needs to be attached to the cuts
1046  Fill(rpCuts,poiCuts);
1047 }
1048 
1049 //-------------------------------------------------------------------//
1050 //---- Including PMD tracks as RP --------------------------//
1051 
1052 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
1053  const AliESDPmdTrack *pmdtracks,
1054  const AliCFManager* poiCFManager ):
1055  AliFlowEventSimple(20), fApplyRecentering(kFALSE), fDivSigma(kTRUE), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0), fChi2A(0x0), fChi2C(0x0), fChi3A(0x0), fChi3C(0x0)
1056 
1057 {
1058  // constructor
1059  for(Int_t i(0); i < 9; i++) {
1060  for(Int_t j(0); j < 2; j++) {
1061  for(Int_t k(0); k < 2; k++) {
1062  fMeanQ[i][j][k] = 0.;
1063  fWidthQ[i][j][k] = 0.;
1064  fMeanQv3[i][j][k] = 0.;
1065  fWidthQv3[i][j][k] = 0.;
1066  }
1067  }
1068  }
1069  for(Int_t i(0); i < 5; i++) {
1070  fQxavsV0[i] = 0x0;
1071  fQyavsV0[i] = 0x0;
1072  fQxcvsV0[i] = 0x0;
1073  fQycvsV0[i] = 0x0;
1074  }
1075 
1076  Float_t GetPmdEta(Float_t xPos, Float_t yPos, Float_t zPos);
1077  Float_t GetPmdPhi(Float_t xPos, Float_t yPos);
1078  //Select the particles of interest from the ESD
1079  Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
1080 
1081  //loop over tracks
1082  for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
1083  {
1084  AliESDtrack* pParticle = anInput->GetTrack(itrkN); //get input particle
1085  //check if pParticle passes the cuts
1086  Bool_t poiOK = kTRUE;
1087  if (poiCFManager)
1088  {
1089  poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
1090  poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
1091  }
1092  if (!poiOK) continue;
1093 
1094  //make new AliFLowTrack
1095  AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
1096 
1097  //marking the particles used for the particle of interest (POI) selection:
1098  if(poiOK && poiCFManager)
1099  {
1100  IncrementNumberOfPOIs(1);
1101  pTrack->SetForPOISelection(kTRUE);
1102  pTrack->SetSource(AliFlowTrack::kFromESD);
1103  }
1104 
1105  AddTrack(pTrack);
1106  }//end of while (itrkN < iNumberOfInputTracks)
1107 
1108  //Select the reference particles from the PMD tracks
1109  Int_t npmdcl = anInput->GetNumberOfPmdTracks();
1110  printf("======There are %d PMD tracks in this event\n-------",npmdcl);
1111  //loop over clusters
1112  for(Int_t iclust=0; iclust < npmdcl; iclust++){
1113  //AliESDPmdTrack *pmdtr = anInput->GetPmdTrack(iclust);
1114  pmdtracks = anInput->GetPmdTrack(iclust);
1115  Int_t det = pmdtracks->GetDetector();
1116  //Int_t smn = pmdtracks->GetSmn();
1117  Float_t clsX = pmdtracks->GetClusterX();
1118  Float_t clsY = pmdtracks->GetClusterY();
1119  Float_t clsZ = pmdtracks->GetClusterZ();
1120  Float_t ncell = pmdtracks->GetClusterCells();
1121  Float_t adc = pmdtracks->GetClusterADC();
1122  //Float_t pid = pmdtracks->GetClusterPID();
1123  Float_t etacls = GetPmdEta(clsX,clsY,clsZ);
1124  Float_t phicls = GetPmdPhi(clsX,clsY);
1125  //make new AliFLowTrackSimple
1126  AliFlowTrack* pTrack = new AliFlowTrack();
1127  //if(det == 0){ //selecting preshower plane only
1128  if(det == 0 && adc > 270 && ncell > 1){ //selecting preshower plane only
1129  //pTrack->SetPt(adc);//cluster adc
1130  pTrack->SetPt(0.0);
1131  pTrack->SetEta(etacls);
1132  pTrack->SetPhi(phicls);
1133  //marking the particles used for the reference particle (RP) selection:
1134  IncrementNumberOfPOIs(0);
1135  pTrack->SetForRPSelection(kTRUE);
1136  pTrack->SetSource(AliFlowTrack::kFromPMD);
1137  //Add the track to the flowevent
1138  AddTrack(pTrack);
1139  }//if det
1140  }
1141 }
1142 //----------------------------------------------------------------------------//
1143 Float_t GetPmdEta(Float_t xPos, Float_t yPos, Float_t zPos)
1144 {
1145  Float_t rpxpy, theta, eta;
1146  rpxpy = TMath::Sqrt(xPos*xPos + yPos*yPos);
1147  theta = TMath::ATan2(rpxpy,zPos);
1148  eta = -TMath::Log(TMath::Tan(0.5*theta));
1149  return eta;
1150 }
1151 //--------------------------------------------------------------------------//
1152 Float_t GetPmdPhi(Float_t xPos, Float_t yPos)
1153 {
1154  Float_t pybypx, phi = 0., phi1;
1155  if(xPos==0)
1156  {
1157  if(yPos>0) phi = 90.;
1158  if(yPos<0) phi = 270.;
1159  }
1160  if(xPos != 0)
1161  {
1162  pybypx = yPos/xPos;
1163  if(pybypx < 0) pybypx = - pybypx;
1164  phi1 = TMath::ATan(pybypx)*180./3.14159;
1165 
1166  if(xPos > 0 && yPos > 0) phi = phi1; // 1st Quadrant
1167  if(xPos < 0 && yPos > 0) phi = 180 - phi1; // 2nd Quadrant
1168  if(xPos < 0 && yPos < 0) phi = 180 + phi1; // 3rd Quadrant
1169  if(xPos > 0 && yPos < 0) phi = 360 - phi1; // 4th Quadrant
1170 
1171  }
1172  phi = phi*3.14159/180.;
1173  return phi;
1174 }
1175 //---------------------------------------------------------------//
1176 AliFlowVector AliFlowEvent::GetQ(Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights)
1177 {
1178  // start with the uncalibrated Q-vector. if we're not looking at vzero data, this will be returned
1179  // so this is completely backward compatible
1180  AliFlowVector vQ = AliFlowEventSimple::GetQ(n, weightsList, usePhiWeights, usePhiWeights, useEtaWeights);
1181  // see if we want to re-center 2010 style
1182  if(fApplyRecentering == 2010) {
1183  AliFlowVector Qarray[2];
1184  AliFlowEvent::Get2Qsub(Qarray, n, weightsList, usePhiWeights, usePtWeights, useEtaWeights);
1185  AliFlowVector vA = Qarray[0];
1186  AliFlowVector vB = Qarray[1];
1187 
1188  // now we have the calibrated sub-event q-vectors, these will be merged using a resolution derived weight
1189  Double_t chiA(1.), chiC(1.), dQX(0.), dQY(0.);
1190  if(n==2) {
1191  chiA = fChi2A->At(fVZEROcentralityBin);
1192  chiC = fChi2C->At(fVZEROcentralityBin);
1193  } else if (n==3) {
1194  chiA = fChi3A->At(fVZEROcentralityBin);
1195  chiC = fChi3C->At(fVZEROcentralityBin);
1196  }
1197 
1198  // combine the vzera and vzeroc signal, note that vzeroa is stored in vB and vzeroc in vA
1199  dQX = chiA*chiA*vB.X()+chiC*chiC*vA.X();
1200  dQY = chiA*chiA*vB.Y()+chiC*chiC*vA.Y();
1201  vQ.Set(dQX, dQY);
1202 
1203  } else if (fApplyRecentering == 2011) { // 11h style recentering
1204  Double_t dQX = 0.;
1205  Double_t dQY = 0.;
1206 
1207  if(fEvent && fEvent->GetEventplane()) {
1208  fEvent->GetEventplane()->CalculateVZEROEventPlane(fEvent, 10, n, dQX, dQY);
1209  // set the new q-vectors (which in this case means REPLACING)
1210  vQ.Set(dQX, dQY);
1211  } // if for some reason the info from the event header is not available, the AliFlowTrackCuts object
1212  // has provided the equalized multiplicity info so this should still be relatively safe
1213  }
1214 
1215  // return the Q-vector
1216  return vQ;
1217 }
1218  //---------------------------------------------------------------//
1219 void AliFlowEvent::Get2Qsub(AliFlowVector* Qarray, Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights)
1220 {
1221  // get q vectors for the subevents. if no recentering is necessary, get the guy from the flow event simple
1222  AliFlowEventSimple::Get2Qsub(Qarray, n, weightsList, usePhiWeights, usePtWeights, useEtaWeights);
1223  AliFlowVector vA = Qarray[0];
1224  AliFlowVector vB = Qarray[1];
1225 
1226  // else get the recentering from the cached info
1227  if (fApplyRecentering == 2010) // 10h style recentering, implemented for n=2 and n=3
1228  {
1229  // first retrieve the q-vectors from the AliFlowEventSimple:: routine
1230  // extract the information form the current flow vectors
1231  Double_t Qxc(vA.X()); // IMPORTANT: user is responsible for the sign of eta
1232  Double_t Qyc(vA.Y()); // vzeroC has negative pseudorapidity and is taken as subevent A
1233  Double_t Qxa(vB.X()); // vzeroA has positive pseudorapidity and is taken as subevent B
1234  Double_t Qya(vB.Y());
1235  // init some values for the corrections
1236 
1237  // default values for vector a (VZEROA)
1238  Double_t Qxamean(0);
1239  Double_t Qxarms(1);
1240  Double_t Qyamean(0);
1241  Double_t Qyarms(1);
1242  // default values for vector b (VZEROC)
1243  Double_t Qxcmean(0);
1244  Double_t Qxcrms(1);
1245  Double_t Qycmean(0);
1246  Double_t Qycrms(1);
1247  // note that defaults are chosen such that for n!=2||n!=3 re-centering is a null-operation
1248 
1249  if( n == 2) { // second order symmetry
1250  Qxamean = fMeanQ[fVZEROcentralityBin][1][0];
1251  Qxarms = fWidthQ[fVZEROcentralityBin][1][0];
1252  Qyamean = fMeanQ[fVZEROcentralityBin][1][1];
1253  Qyarms = fWidthQ[fVZEROcentralityBin][1][1];
1254 
1255  Qxcmean = fMeanQ[fVZEROcentralityBin][0][0];
1256  Qxcrms = fWidthQ[fVZEROcentralityBin][0][0];
1257  Qycmean = fMeanQ[fVZEROcentralityBin][0][1];
1258  Qycrms = fWidthQ[fVZEROcentralityBin][0][1];
1259  } else if (n == 3) { // third order symmetry
1260  Qxamean = fMeanQv3[fVZEROcentralityBin][1][0];
1261  Qxarms = fWidthQv3[fVZEROcentralityBin][1][0];
1262  Qyamean = fMeanQv3[fVZEROcentralityBin][1][1];
1263  Qyarms = fWidthQv3[fVZEROcentralityBin][1][1];
1264 
1265  Qxcmean = fMeanQv3[fVZEROcentralityBin][0][0];
1266  Qxcrms = fWidthQv3[fVZEROcentralityBin][0][0];
1267  Qycmean = fMeanQv3[fVZEROcentralityBin][0][1];
1268  Qycrms = fWidthQv3[fVZEROcentralityBin][0][1];
1269  }
1270  // do the correction
1271  Double_t QxaCor = (Qxa - Qxamean)/Qxarms;
1272  Double_t QyaCor = (Qya - Qyamean)/Qyarms;
1273  Double_t QxcCor = (Qxc - Qxcmean)/Qxcrms;
1274  Double_t QycCor = (Qyc - Qycmean)/Qycrms;
1275  // update the vector
1276  vA.Set(QxcCor, QycCor);
1277  vB.Set(QxaCor, QyaCor);
1278  } else if (fApplyRecentering == 2011) { // 11h style recentering
1279 
1280  Double_t QxaCor = 0.;
1281  Double_t QyaCor = 0.;
1282  Double_t QxcCor = 0.;
1283  Double_t QycCor = 0.;
1284 
1285  if(fEvent && fEvent->GetEventplane()) {
1286  fEvent->GetEventplane()->CalculateVZEROEventPlane(fEvent, 8, n, QxaCor, QyaCor);
1287  fEvent->GetEventplane()->CalculateVZEROEventPlane(fEvent, 9, n, QxcCor, QycCor);
1288 
1289  // set the new q-vectors (which in this case means REPLACING)
1290  vA.Set(QxcCor, QycCor);
1291  vB.Set(QxaCor, QyaCor);
1292  } // if for some reason the info from the event header is not available, the AliFlowTrackCuts object
1293  // has provided the equalized multiplicity info so this should still be relatively safe
1294  } else if (fApplyRecentering == 999) {
1295  // experimental VZERO recentering
1296  // first retrieve the q-vectors from the AliFlowEventSimple:: routine
1297  // extract the information form the current flow vectors
1298  Double_t Qxc(vA.X()); // IMPORTANT: user is responsible for the sign of eta
1299  Double_t Qyc(vA.Y()); // vzeroC has negative pseudorapidity and is taken as subevent A
1300  Double_t Qxa(vB.X()); // vzeroA has positive pseudorapidity and is taken as subevent B
1301  Double_t Qya(vB.Y());
1302  // init some values for the corrections
1303 
1304  Double_t Cen = fEvent->GetCentrality()->GetCentralityPercentile("V0M");
1305  // default values for vector a (VZEROA)
1306  Double_t Qxamean(fQxavsV0[n-1]->GetBinContent(fQxavsV0[n-1]->FindBin(Cen)));
1307  Double_t Qxarms(fQxavsV0[n-1]->GetBinError(fQxavsV0[n-1]->FindBin(Cen)));
1308  Double_t Qyamean(fQyavsV0[n-1]->GetBinContent(fQyavsV0[n-1]->FindBin(Cen)));
1309  Double_t Qyarms(fQyavsV0[n-1]->GetBinError(fQyavsV0[n-1]->FindBin(Cen)));
1310  // default values for vector b (VZEROC)
1311  Double_t Qxcmean(fQxcvsV0[n-1]->GetBinContent(fQxcvsV0[n-1]->FindBin(Cen)));
1312  Double_t Qxcrms(fQxcvsV0[n-1]->GetBinError(fQxcvsV0[n-1]->FindBin(Cen)));
1313  Double_t Qycmean(fQycvsV0[n-1]->GetBinContent(fQycvsV0[n-1]->FindBin(Cen)));
1314  Double_t Qycrms(fQycvsV0[n-1]->GetBinError(fQycvsV0[n-1]->FindBin(Cen)));
1315 
1316  // just a precaution ...
1317  if(n > 5) {
1318  Qxamean = 0;
1319  Qxarms = 1;
1320  Qyamean = 0;
1321  Qyarms = 1;
1322  Qxcmean = 0;
1323  Qxcrms = 1;
1324  Qycmean = 0; // this effectively disables the calibration
1325  Qycrms = 1;
1326  }
1327 
1328  Double_t QxaR = Qxa - Qxamean;
1329  Double_t QyaR = Qya - Qyamean;
1330  Double_t QxcR = Qxc - Qxcmean;
1331  Double_t QycR = Qyc - Qycmean;
1332  if(fDivSigma && Qxarms>0. && Qyarms>0. && Qxcrms>0. && Qycrms>0.) {
1333  QxaR /= Qxarms;
1334  QyaR /= Qyarms;
1335  QxcR /= Qxcrms;
1336  QycR /= Qycrms;
1337  }
1338  // update the vector
1339  vA.Set(QxcR, QycR);
1340  vB.Set(QxaR, QyaR);
1341 
1342  } else if (fApplyRecentering == 666) {
1343  // experimental VZERO recentering for 11h Full TPC Flow data, harmonics 1,2,3
1344  // first retrieve the q-vectors from the AliFlowEventSimple:: routine
1345  // extract the information form the current flow vectors
1346  if(n < 4) {
1347  Double_t Qxc(vA.X()); // IMPORTANT: user is responsible for the sign of eta
1348  Double_t Qyc(vA.Y()); // vzeroC has negative pseudorapidity and is taken as subevent A
1349  Double_t Qxa(vB.X()); // vzeroA has positive pseudorapidity and is taken as subevent B
1350  Double_t Qya(vB.Y());
1351  Double_t MultC = vA.GetMult();
1352  Double_t MultA = vB.GetMult();
1353  // init some values for the corrections
1354 
1355  Double_t Cen = fEvent->GetCentrality()->GetCentralityPercentile("V0M");
1356  // default values for vector a (VZEROA)
1357  Double_t Qxamean(fQxavsV0[n-1]->GetBinContent(fQxavsV0[n-1]->FindBin(Cen)));
1358  Double_t Qxarms(fQxavsV0[n-1]->GetBinError(fQxavsV0[n-1]->FindBin(Cen)));
1359  Double_t Qyamean(fQyavsV0[n-1]->GetBinContent(fQyavsV0[n-1]->FindBin(Cen)));
1360  Double_t Qyarms(fQyavsV0[n-1]->GetBinError(fQyavsV0[n-1]->FindBin(Cen)));
1361  // default values for vector b (VZEROC)
1362  Double_t Qxcmean(fQxcvsV0[n-1]->GetBinContent(fQxcvsV0[n-1]->FindBin(Cen)));
1363  Double_t Qxcrms(fQxcvsV0[n-1]->GetBinError(fQxcvsV0[n-1]->FindBin(Cen)));
1364  Double_t Qycmean(fQycvsV0[n-1]->GetBinContent(fQycvsV0[n-1]->FindBin(Cen)));
1365  Double_t Qycrms(fQycvsV0[n-1]->GetBinError(fQycvsV0[n-1]->FindBin(Cen)));
1366 
1367  Double_t QxaR = Qxa - Qxamean*MultA;
1368  Double_t QyaR = Qya - Qyamean*MultA;
1369  Double_t QxcR = Qxc - Qxcmean*MultC;
1370  Double_t QycR = Qyc - Qycmean*MultC;
1371  if(fDivSigma && Qxarms>0. && Qyarms>0. && Qxcrms>0. && Qycrms>0.) {
1372  QxaR /= Qxarms;
1373  QyaR /= Qyarms;
1374  QxcR /= Qxcrms;
1375  QycR /= Qycrms;
1376  }
1377  // update the vector
1378  vA.Set(QxcR, QycR);
1379  vB.Set(QxaR, QyaR);
1380 
1381  } else {
1382  cout << " WARNING: recentering not possible for harmonic " << n << " (delta calibration) " << endl;
1383  } // end of if(n < 4)
1384  }
1385 
1386  Qarray[0] = vA;
1387  Qarray[1] = vB;
1388 }
1389 //_____________________________________________________________________________
1390 void AliFlowEvent::SetVZEROCalibrationForTrackCuts(AliFlowTrackCuts* cuts) {
1391  // open calibration info, copied from AliAnalyisTaskVnV0.cxx
1392  fEvent = cuts->GetEvent();
1393  if(!fEvent) return; // coverity. we need to know the event to get the runnumber and centrlaity
1394  // get the vzero centrality percentile (cc dependent calibration)
1395  Float_t v0Centr(fEvent->GetCentrality()->GetCentralityPercentile("V0M"));
1396  if(v0Centr < 5) fVZEROcentralityBin = 0;
1397  else if(v0Centr < 10) fVZEROcentralityBin = 1;
1398  else if(v0Centr < 20) fVZEROcentralityBin = 2;
1399  else if(v0Centr < 30) fVZEROcentralityBin = 3;
1400  else if(v0Centr < 40) fVZEROcentralityBin = 4;
1401  else if(v0Centr < 50) fVZEROcentralityBin = 5;
1402  else if(v0Centr < 60) fVZEROcentralityBin = 6;
1403  else if(v0Centr < 70) fVZEROcentralityBin = 7;
1404  else fVZEROcentralityBin = 8;
1405 
1406  // if this event is from the same run as the previous event
1407  // we can use the cached calibration values, no need to re-open the
1408  // aodb file, else cache the new run
1409  Int_t run(fEvent->GetRunNumber());
1410  if(fCachedRun == run) return;
1411  else fCachedRun = run;
1412 
1413  // relevant for 2010 data: check if the proper chi weights for merging vzero a and vzero c ep are present
1414  // if not, use sane defaults. centrality binning is equal to that given in the fVZEROcentralityBin snippet
1415  //
1416  // chi values can be calculated using the static helper function
1417  // AliAnalysisTaskJetV2::CalculateEventPlaneChi(Double_t res) where res is the event plane
1418  // resolution in a given centrality bin
1419  //
1420  // the resolutions that were used for these defaults are
1421  // Double_t R2VZEROA[] = {.35, .40, .48, .50, .48, .45, .38, .26, .16};
1422  // Double_t R2VZEROC[] = {.45, .60, .70, .73, .68, .60, .40, .36, .17};
1423  //
1424  // Double_t R3VZEROA[] = {.22, .23, .22, .19, .15, .12, .08, .00, .00};
1425  // Double_t R3VZEROC[] = {.30, .30, .28, .25, .22, .17, .11, .00, .00};
1426  // this might need a bit of updating as they were read 'by-eye' from a performance plot ..
1427  Double_t chiC2[] = {0.771423, 1.10236, 1.38116, 1.48077, 1.31964, 1.10236, 0.674622, 0.600403, 0.273865};
1428  Double_t chiA2[] = {0.582214, 0.674622, 0.832214, 0.873962, 0.832214, 0.771423, 0.637146, 0.424255, 0.257385};
1429  Double_t chiC3[] = {0.493347, 0.493347, 0.458557, 0.407166, 0.356628, 0.273865, 0.176208, 6.10352e-05, 6.10352e-05};
1430  Double_t chiA3[] = {0.356628, 0.373474, 0.356628, 0.306702, 0.24115, 0.192322, 0.127869, 6.10352e-05, 6.10352e-05};
1431 
1432  // this may seem redundant but in this way the cuts object is owner of the arrays
1433  // even if they're created here (so we won't get into trouble with dtor, assigmnet and copying)
1434  if(!cuts->GetChi2A()) cuts->SetChi2A(new TArrayD(9, chiA2));
1435  if(!cuts->GetChi2C()) cuts->SetChi2C(new TArrayD(9, chiC2));
1436  if(!cuts->GetChi3A()) cuts->SetChi3A(new TArrayD(9, chiA3));
1437  if(!cuts->GetChi3C()) cuts->SetChi3C(new TArrayD(9, chiC3));
1438 
1439  if(!fChi2A) fChi2A = cuts->GetChi2A();
1440  if(!fChi2C) fChi2C = cuts->GetChi2C();
1441  if(!fChi3A) fChi3A = cuts->GetChi3A();
1442  if(!fChi3C) fChi3C = cuts->GetChi3C();
1443 
1444  TFile *foadb = TFile::Open("$ALICE_PHYSICS/OADB/PWGCF/VZERO/VZEROcalibEP.root");
1445  if(!foadb){
1446  printf("OADB file $ALICE_PHYSICS/OADB/PWGCF/VZERO/VZEROcalibEP.root cannot be opened, CALIBRATION FAILED !");
1447  return;
1448  }
1449 
1450  AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
1451  if(!cont){
1452  printf("OADB object hMultV0BefCorr is not available in the file\n");
1453  return;
1454  }
1455  if(!(cont->GetObject(run))){
1456  // if the multiplicity correction cannot be found for the specified run,
1457  // loop over the 11h runs to see if it's 11h data
1458  Int_t runs11h[] = {170593, 170572, 170556, 170552, 170546, 170390, 170389, 170388, 170387, 170315, 170313, 170312, 170311, 170309, 170308, 170306, 170270, 170269, 170268, 170267, 170264, 170230, 170228, 170208, 170207, 170205, 170204, 170203, 170195, 170193, 170163, 170162, 170159, 170155, 170152, 170091, 170089, 170088, 170085, 170084, 170083, 170081, 170040, 170038, 170036, 170027, 169981, 169975, 169969, 169965, 169961, 169956, 169926, 169924, 169923, 169922, 169919, 169918, 169914, 169859, 169858, 169855, 169846, 169838, 169837, 169835, 169683, 169628, 169591, 169590, 169588, 169587, 169586, 169584, 169557, 169555, 169554, 169553, 169550, 169515, 169512, 169506, 169504, 169498, 169475, 169420, 169419, 169418, 169417, 169415, 169411, 169238, 169236, 169167, 169160, 169156, 169148, 169145, 169144, 169143, 169138, 169099, 169094, 169091, 169045, 169044, 169040, 169035, 168992, 168988, 168984, 168826, 168777, 168514, 168512, 168511, 168467, 168464, 168461, 168460, 168458, 168362, 168361, 168356, 168342, 168341, 168325, 168322, 168318, 168311, 168310, 168213, 168212, 168208, 168207, 168206, 168205, 168204, 168203, 168181, 168177, 168175, 168173, 168172, 168171, 168115, 168108, 168107, 168105, 168104, 168103, 168076, 168069, 168068, 168066, 167988, 167987, 167986, 167985, 167921, 167920, 167915, 167909, 167903, 167902, 167818, 167814, 167813, 167808, 167807, 167806, 167713, 167712, 167711, 167706, 167693};
1459  for(Int_t r(0); r < 176; r++) {
1460  if(run == runs11h[r]) {
1461  printf(" > run has been identified as 11h < \n");
1462  if(cuts->GetVZEROgainEqualizationPerRing()) {
1463  // enable or disable rings through the weights, weight 1. is enabled, 0. is disabled
1464  // start with the vzero c rings (segments 0 through 31)
1465  (cuts->GetUseVZERORing(0)) ? cuts->SetVZEROCpol(0, 1.) : cuts->SetVZEROCpol(0, 0.);
1466  (cuts->GetUseVZERORing(1)) ? cuts->SetVZEROCpol(1, 1.) : cuts->SetVZEROCpol(1, 0.);
1467  (cuts->GetUseVZERORing(2)) ? cuts->SetVZEROCpol(2, 1.) : cuts->SetVZEROCpol(2, 0.);
1468  (cuts->GetUseVZERORing(3)) ? cuts->SetVZEROCpol(3, 1.) : cuts->SetVZEROCpol(3, 0.);
1469  // same for vzero a
1470  (cuts->GetUseVZERORing(4)) ? cuts->SetVZEROApol(0, 1.) : cuts->SetVZEROApol(0, 0.);
1471  (cuts->GetUseVZERORing(5)) ? cuts->SetVZEROApol(1, 1.) : cuts->SetVZEROApol(1, 0.);
1472  (cuts->GetUseVZERORing(6)) ? cuts->SetVZEROApol(2, 1.) : cuts->SetVZEROApol(2, 0.);
1473  (cuts->GetUseVZERORing(7)) ? cuts->SetVZEROApol(3, 1.) : cuts->SetVZEROApol(3, 0.);
1474  } else {
1475  // else enable all rings, which is also default
1476  for(Int_t i(0); i < 4; i++) cuts->SetVZEROCpol(i, 1.);
1477  for(Int_t i(0); i < 4; i++) cuts->SetVZEROApol(i, 1.);
1478  }
1479  // pass a NULL pointer to the track cuts object, the NULL pointer will identify 11h runs
1480  cuts->SetVZEROgainEqualisation(NULL);
1481  fApplyRecentering = 2011;
1482  return; // the rest of the steps are not necessary
1483  }
1484  }
1485  // the run has not been identified as lhc11h data, so we assume a template calibration
1486  printf("OADB object hMultVZEROBefCorr is not available for run %i (used default run 137366)\n",run);
1487  run = 137366;
1488  }
1489  printf(" > run has been identified as 10h < \n");
1490  // step 1) get the proper multiplicity weights from the vzero signal
1491  TProfile* fMultVZERO = ((TH2F *) cont->GetObject(run))->ProfileX();
1492 
1493  TF1 *fpol0 = new TF1("fpol0","pol0");
1494  if(cuts->GetVZEROgainEqualizationPerRing()) {
1495  // do the calibration per ring
1496  // start with the vzero c rings (segments 0 through 31)
1497  fMultVZERO->Fit(fpol0, "N0", "", 0, 8);
1498  (cuts->GetUseVZERORing(0)) ? cuts->SetVZEROCpol(0, fpol0->GetParameter(0)) : cuts->SetVZEROCpol(0, 0.);
1499  fMultVZERO->Fit(fpol0, "N0", "", 8, 16);
1500  (cuts->GetUseVZERORing(1)) ? cuts->SetVZEROCpol(1, fpol0->GetParameter(0)) : cuts->SetVZEROCpol(1, 0.);
1501  fMultVZERO->Fit(fpol0, "N0", "", 16, 24);
1502  (cuts->GetUseVZERORing(2)) ? cuts->SetVZEROCpol(2, fpol0->GetParameter(0)) : cuts->SetVZEROCpol(2, 0.);
1503  fMultVZERO->Fit(fpol0, "N0", "", 24, 32);
1504  (cuts->GetUseVZERORing(3)) ? cuts->SetVZEROCpol(3, fpol0->GetParameter(0)) : cuts->SetVZEROCpol(3, 0.);
1505  // same thing for vero A
1506  fMultVZERO->Fit(fpol0, "N0", "", 32, 40);
1507  (cuts->GetUseVZERORing(4)) ? cuts->SetVZEROApol(0, fpol0->GetParameter(0)) : cuts->SetVZEROApol(0, 0.);
1508  fMultVZERO->Fit(fpol0, "N0", "", 40, 48);
1509  (cuts->GetUseVZERORing(5)) ? cuts->SetVZEROApol(1, fpol0->GetParameter(0)) : cuts->SetVZEROApol(1, 0.);
1510  fMultVZERO->Fit(fpol0, "N0", "", 48, 56);
1511  (cuts->GetUseVZERORing(6)) ? cuts->SetVZEROApol(2, fpol0->GetParameter(0)) : cuts->SetVZEROApol(2, 0.);
1512  fMultVZERO->Fit(fpol0, "N0", "", 56, 64);
1513  (cuts->GetUseVZERORing(7)) ? cuts->SetVZEROApol(3, fpol0->GetParameter(0)) : cuts->SetVZEROApol(3, 0.);
1514  } else {
1515  // do the calibration in one go. the calibration will still be
1516  // stored per ring, but each ring has the same weight now
1517  fMultVZERO->Fit(fpol0,"N0","",0,31);
1518  for(Int_t i(0); i < 4; i++) cuts->SetVZEROCpol(i, fpol0->GetParameter(0));
1519  fMultVZERO->Fit(fpol0,"N0","",32,64);
1520  for(Int_t i(0); i < 4; i++) cuts->SetVZEROApol(i, fpol0->GetParameter(0));
1521  }
1522  // the parameters to weigh the vzero track cuts have been extracted now,
1523  // so we can pass them to the current track cuts obect
1524  cuts->SetVZEROgainEqualisation(fMultVZERO); // passed as a TH1
1525 
1526  // step 2) reweight the q-vectors that will be called by flow methods which use
1527  // subevents
1528  // underlying assumption is that subevent a uses VZEROA
1529  // and subevent b uses VZEROC
1530  for(Int_t iside=0;iside<2;iside++){
1531  for(Int_t icoord=0;icoord<2;icoord++){
1532  for(Int_t i=0;i < 9;i++){
1533  char namecont[100];
1534  if(iside==0 && icoord==0)
1535  snprintf(namecont,100,"hQxc2_%i",i);
1536  else if(iside==1 && icoord==0)
1537  snprintf(namecont,100,"hQxa2_%i",i);
1538  else if(iside==0 && icoord==1)
1539  snprintf(namecont,100,"hQyc2_%i",i);
1540  else if(iside==1 && icoord==1)
1541  snprintf(namecont,100,"hQya2_%i",i);
1542 
1543  cont = (AliOADBContainer*) foadb->Get(namecont);
1544  if(!cont){
1545  printf("OADB object %s is not available in the file\n",namecont);
1546  return;
1547  }
1548 
1549  if(!(cont->GetObject(run))){
1550  printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
1551  run = 137366;
1552  }
1553 
1554  // after grabbing all the info, set the CORRECTION TERMS to
1555  // the 2nd and 3rd order qsub-vectors
1556  // we do this here for all centralities, so that subsequent events
1557  // can grab the correction from these cached values
1558  fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
1559  fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
1560 
1561  //for v3
1562  if(iside==0 && icoord==0)
1563  snprintf(namecont,100,"hQxc3_%i",i);
1564  else if(iside==1 && icoord==0)
1565  snprintf(namecont,100,"hQxa3_%i",i);
1566  else if(iside==0 && icoord==1)
1567  snprintf(namecont,100,"hQyc3_%i",i);
1568  else if(iside==1 && icoord==1)
1569  snprintf(namecont,100,"hQya3_%i",i);
1570 
1571  cont = (AliOADBContainer*) foadb->Get(namecont);
1572  if(!cont){
1573  printf("OADB object %s is not available in the file\n",namecont);
1574  return;
1575  }
1576 
1577  if(!(cont->GetObject(run))){
1578  printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
1579  run = 137366;
1580  }
1581  fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
1582  fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
1583 
1584  }
1585  }
1586  }
1587  // set the recentering style (might be switched back to -1 if recentering is disabeled)
1588  fApplyRecentering = 2010;
1589 }
1590 //-----------------------------------------------------------------------
1591 void AliFlowEvent::SetBetaVZEROCalibrationForTrackCuts(AliFlowTrackCuts* cuts) {
1592  // implementation of beta vzero calibration
1593 
1594  fEvent = cuts->GetEvent();
1595  if(!fEvent) return; // coverity. we need to know the event to get the runnumber and centrlaity
1596  // get the vzero centrality percentile (cc dependent calibration)
1597 
1598 
1599  // if this event is from the same run as the previous event
1600  // we can use the cached calibration values, no need to re-open the
1601  // aodb file, else cache the new run
1602  Int_t run(fEvent->GetRunNumber());
1603  if(fCachedRun == run) return;
1604  else fCachedRun = run;
1605 
1606  // check if the proper chi weights for merging vzero a and vzero c ep are present
1607  // if not, use sane defaults. centrality binning is equal to that given in the fVZEROcentralityBin snippet
1608  //
1609  // chi values can be calculated using the static helper function
1610  // AliAnalysisTaskJetV2::CalculateEventPlaneChi(Double_t res) where res is the event plane
1611  // resolution in a given centrality bin
1612  //
1613  // the resolutions that were used for these defaults are
1614  // Double_t R2VZEROA[] = {.35, .40, .48, .50, .48, .45, .38, .26, .16};
1615  // Double_t R2VZEROC[] = {.45, .60, .70, .73, .68, .60, .40, .36, .17};
1616  //
1617  // Double_t R3VZEROA[] = {.22, .23, .22, .19, .15, .12, .08, .00, .00};
1618  // Double_t R3VZEROC[] = {.30, .30, .28, .25, .22, .17, .11, .00, .00};
1619  Double_t chiC2[] = {0.771423, 1.10236, 1.38116, 1.48077, 1.31964, 1.10236, 0.674622, 0.600403, 0.273865};
1620  Double_t chiA2[] = {0.582214, 0.674622, 0.832214, 0.873962, 0.832214, 0.771423, 0.637146, 0.424255, 0.257385};
1621  Double_t chiC3[] = {0.493347, 0.493347, 0.458557, 0.407166, 0.356628, 0.273865, 0.176208, 6.10352e-05, 6.10352e-05};
1622  Double_t chiA3[] = {0.356628, 0.373474, 0.356628, 0.306702, 0.24115, 0.192322, 0.127869, 6.10352e-05, 6.10352e-05};
1623 
1624  // this may seem redundant but in this way the cuts object is owner of the arrays
1625  // even if they're created here (so we won't get into trouble with dtor, assigmnet and copying)
1626  if(!cuts->GetChi2A()) cuts->SetChi2A(new TArrayD(9, chiA2));
1627  if(!cuts->GetChi2C()) cuts->SetChi2C(new TArrayD(9, chiC2));
1628  if(!cuts->GetChi3A()) cuts->SetChi3A(new TArrayD(9, chiA3));
1629  if(!cuts->GetChi3C()) cuts->SetChi3C(new TArrayD(9, chiC3));
1630 
1631  if(!fChi2A) fChi2A = cuts->GetChi2A();
1632  if(!fChi2C) fChi2C = cuts->GetChi2C();
1633  if(!fChi3A) fChi3A = cuts->GetChi3A();
1634  if(!fChi3C) fChi3C = cuts->GetChi3C();
1635 
1636 
1637  TFile *foadb = TFile::Open("$ALICE_PHYSICS/PWGCF/FLOW/database/calibV0_filtered.root");
1638  if(!foadb){
1639  printf("OADB file $ALICE_PHYSICS/PWGCF/FLOW/database/calibV0_filtered.root cannot be opened, CALIBRATION FAILED !");
1640  return;
1641  }
1642 
1643  // first get the mutiplicity of the vzero channels before gain equalization
1644  AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr_filtered");
1645  if(!cont){
1646  printf("OADB object hMultV0BefCorr is not available in the file\n");
1647  return;
1648  }
1649  if(!(cont->GetObject(run))) {
1650  // if the multiplicity correction cannot be found for the specified run,
1651  // loop over the 11h runs to see if it's 11h data
1652  Int_t runs11h[] = {170593, 170572, 170556, 170552, 170546, 170390, 170389, 170388, 170387, 170315, 170313, 170312, 170311, 170309, 170308, 170306, 170270, 170269, 170268, 170267, 170264, 170230, 170228, 170208, 170207, 170205, 170204, 170203, 170195, 170193, 170163, 170162, 170159, 170155, 170152, 170091, 170089, 170088, 170085, 170084, 170083, 170081, 170040, 170038, 170036, 170027, 169981, 169975, 169969, 169965, 169961, 169956, 169926, 169924, 169923, 169922, 169919, 169918, 169914, 169859, 169858, 169855, 169846, 169838, 169837, 169835, 169683, 169628, 169591, 169590, 169588, 169587, 169586, 169584, 169557, 169555, 169554, 169553, 169550, 169515, 169512, 169506, 169504, 169498, 169475, 169420, 169419, 169418, 169417, 169415, 169411, 169238, 169236, 169167, 169160, 169156, 169148, 169145, 169144, 169143, 169138, 169099, 169094, 169091, 169045, 169044, 169040, 169035, 168992, 168988, 168984, 168826, 168777, 168514, 168512, 168511, 168467, 168464, 168461, 168460, 168458, 168362, 168361, 168356, 168342, 168341, 168325, 168322, 168318, 168311, 168310, 168213, 168212, 168208, 168207, 168206, 168205, 168204, 168203, 168181, 168177, 168175, 168173, 168172, 168171, 168115, 168108, 168107, 168105, 168104, 168103, 168076, 168069, 168068, 168066, 167988, 167987, 167986, 167985, 167921, 167920, 167915, 167909, 167903, 167902, 167818, 167814, 167813, 167808, 167807, 167806, 167713, 167712, 167711, 167706, 167693};
1653  for(Int_t r(0); r < 176; r++) {
1654  if(run == runs11h[r]) {
1655  printf(" > run has been identified as 11h < \n");
1656  if(cuts->GetVZEROgainEqualizationPerRing()) {
1657  // enable or disable rings through the weights, weight 1. is enabled, 0. is disabled
1658  // start with the vzero c rings (segments 0 through 31)
1659  (cuts->GetUseVZERORing(0)) ? cuts->SetVZEROCpol(0, 1.) : cuts->SetVZEROCpol(0, 0.);
1660  (cuts->GetUseVZERORing(1)) ? cuts->SetVZEROCpol(1, 1.) : cuts->SetVZEROCpol(1, 0.);
1661  (cuts->GetUseVZERORing(2)) ? cuts->SetVZEROCpol(2, 1.) : cuts->SetVZEROCpol(2, 0.);
1662  (cuts->GetUseVZERORing(3)) ? cuts->SetVZEROCpol(3, 1.) : cuts->SetVZEROCpol(3, 0.);
1663  // same for vzero a
1664  (cuts->GetUseVZERORing(4)) ? cuts->SetVZEROApol(0, 1.) : cuts->SetVZEROApol(0, 0.);
1665  (cuts->GetUseVZERORing(5)) ? cuts->SetVZEROApol(1, 1.) : cuts->SetVZEROApol(1, 0.);
1666  (cuts->GetUseVZERORing(6)) ? cuts->SetVZEROApol(2, 1.) : cuts->SetVZEROApol(2, 0.);
1667  (cuts->GetUseVZERORing(7)) ? cuts->SetVZEROApol(3, 1.) : cuts->SetVZEROApol(3, 0.);
1668  } else {
1669  // else enable all rings, which is also default
1670  for(Int_t i(0); i < 4; i++) cuts->SetVZEROCpol(i, 1.);
1671  for(Int_t i(0); i < 4; i++) cuts->SetVZEROApol(i, 1.);
1672  }
1673  // pass a NULL pointer to the track cuts object, the NULL pointer will identify 11h runs
1674  cuts->SetVZEROgainEqualisation(NULL);
1675  fApplyRecentering = 2011;
1676  return; // the rest of the steps are not necessary
1677  }
1678  }
1679  // the run has not been identified as lhc11h data, so we assume a template calibration
1680  printf("OADB object hMultVZEROBefCorr is not available for run %i (used default run 138275)\n",run);
1681  run = 138275;
1682  }
1683  printf(" > run has been identified as 10h < \n");
1684  // step 0) get the profile which contains average multiplicity per VZERO channel
1685  TProfile* fMultVZERO = static_cast<TProfile*>(cont->GetObject(run));
1686 
1687  TF1 *fpol0 = new TF1("fpol0","pol0");
1688  // step 1) extract the proper weights from the profile. Q-vector recentering relies on
1689  // ring-by-ring gain equalization, so the terms are extracted ring-by-ring here
1690 
1691  // start with the vzero c rings (segments 0 through 31)
1692  fMultVZERO->Fit(fpol0, "N0", "", 0, 8);
1693  (cuts->GetUseVZERORing(0)) ? cuts->SetVZEROCpol(0, fpol0->GetParameter(0)) : cuts->SetVZEROCpol(0, 0.);
1694  fMultVZERO->Fit(fpol0, "N0", "", 8, 16);
1695  (cuts->GetUseVZERORing(1)) ? cuts->SetVZEROCpol(1, fpol0->GetParameter(0)) : cuts->SetVZEROCpol(1, 0.);
1696  fMultVZERO->Fit(fpol0, "N0", "", 16, 24);
1697  (cuts->GetUseVZERORing(2)) ? cuts->SetVZEROCpol(2, fpol0->GetParameter(0)) : cuts->SetVZEROCpol(2, 0.);
1698  fMultVZERO->Fit(fpol0, "N0", "", 24, 32);
1699  (cuts->GetUseVZERORing(3)) ? cuts->SetVZEROCpol(3, fpol0->GetParameter(0)) : cuts->SetVZEROCpol(3, 0.);
1700  // same thing for vero A
1701  fMultVZERO->Fit(fpol0, "N0", "", 32, 40);
1702  (cuts->GetUseVZERORing(4)) ? cuts->SetVZEROApol(0, fpol0->GetParameter(0)) : cuts->SetVZEROApol(0, 0.);
1703  fMultVZERO->Fit(fpol0, "N0", "", 40, 48);
1704  (cuts->GetUseVZERORing(5)) ? cuts->SetVZEROApol(1, fpol0->GetParameter(0)) : cuts->SetVZEROApol(1, 0.);
1705  fMultVZERO->Fit(fpol0, "N0", "", 48, 56);
1706  (cuts->GetUseVZERORing(6)) ? cuts->SetVZEROApol(2, fpol0->GetParameter(0)) : cuts->SetVZEROApol(2, 0.);
1707  fMultVZERO->Fit(fpol0, "N0", "", 56, 64);
1708  (cuts->GetUseVZERORing(7)) ? cuts->SetVZEROApol(3, fpol0->GetParameter(0)) : cuts->SetVZEROApol(3, 0.);
1709 
1710  // the parameters to weigh the vzero track cuts have been extracted now,
1711  // so we can pass them to the current track cuts obect
1712  cuts->SetVZEROgainEqualisation(fMultVZERO); // passed as a TH1
1713 
1714  // step 2) extract the calibration histograms from the database and
1715  // pass them to the cuts object
1716  //
1717  // first index of the oadb array is the harmonic n, the second index is either qax, qay, qcx, qcy
1718  AliOADBContainer* h[5][4];
1719  for(Int_t i(0); i < 5; i++) {
1720  h[i][0] = (AliOADBContainer*)foadb->Get(Form("hQxa%i_filtered", i+1));
1721  if(h[i][0]) fQxavsV0[i] = static_cast<TH1F*>(h[i][0]->GetObject(run));
1722  h[i][1] = (AliOADBContainer*)foadb->Get(Form("hQya%i_filtered", i+1));
1723  if(h[i][1]) fQyavsV0[i] = static_cast<TH1F*>(h[i][1]->GetObject(run));
1724  h[i][2] = (AliOADBContainer*)foadb->Get(Form("hQxc%i_filtered", i+1));
1725  if(h[i][2]) fQxcvsV0[i] = static_cast<TH1F*>(h[i][2]->GetObject(run));
1726  h[i][3] = (AliOADBContainer*)foadb->Get(Form("hQyc%i_filtered", i+1));
1727  if(h[i][3]) fQycvsV0[i] = static_cast<TH1F*>(h[i][3]->GetObject(run));
1728  }
1729 
1730  // set the recentering style (might be switched back to -1 if recentering is disabeled)
1731  // FIXME as an ugly hack, for now I mark this as 999 to denote the experimental nature
1732  // of this and use it transparently without disrupting the existing calbiration
1733  fApplyRecentering = 999;
1734 }
1735 
1736 //-----------------------------------------------------------------------------
1737 
1738 void AliFlowEvent::SetDeltaVZEROCalibrationForTrackCuts(AliFlowTrackCuts* cuts) {
1739  // implementation of delta vzero calibration 2011
1740 
1741  fEvent = cuts->GetEvent();
1742  if(!fEvent) return; // coverity. we need to know the event to get the runnumber and centrlaity
1743  // get the vzero centrality percentile (cc dependent calibration)
1744 
1745 
1746  // if this event is from the same run as the previous event
1747  // we can use the cached calibration values, no need to re-open the
1748  // aodb file, else cache the new run
1749  Int_t run(fEvent->GetRunNumber());
1750  if(fCachedRun == run) return;
1751  else fCachedRun = run;
1752 
1753  // check if the proper chi weights for merging vzero a and vzero c ep are present
1754  // if not, use sane defaults. centrality binning is equal to that given in the fVZEROcentralityBin snippet
1755  //
1756  // chi values can be calculated using the static helper function
1757  // AliAnalysisTaskJetV2::CalculateEventPlaneChi(Double_t res) where res is the event plane
1758  // resolution in a given centrality bin
1759  //
1760  // the resolutions that were used for these defaults are
1761  // Double_t R2VZEROA[] = {.35, .40, .48, .50, .48, .45, .38, .26, .16};
1762  // Double_t R2VZEROC[] = {.45, .60, .70, .73, .68, .60, .40, .36, .17};
1763  //
1764  // Double_t R3VZEROA[] = {.22, .23, .22, .19, .15, .12, .08, .00, .00};
1765  // Double_t R3VZEROC[] = {.30, .30, .28, .25, .22, .17, .11, .00, .00};
1766  Double_t chiC2[] = {0.771423, 1.10236, 1.38116, 1.48077, 1.31964, 1.10236, 0.674622, 0.600403, 0.273865};
1767  Double_t chiA2[] = {0.582214, 0.674622, 0.832214, 0.873962, 0.832214, 0.771423, 0.637146, 0.424255, 0.257385};
1768  Double_t chiC3[] = {0.493347, 0.493347, 0.458557, 0.407166, 0.356628, 0.273865, 0.176208, 6.10352e-05, 6.10352e-05};
1769  Double_t chiA3[] = {0.356628, 0.373474, 0.356628, 0.306702, 0.24115, 0.192322, 0.127869, 6.10352e-05, 6.10352e-05};
1770 
1771  // this may seem redundant but in this way the cuts object is owner of the arrays
1772  // even if they're created here (so we won't get into trouble with dtor, assigmnet and copying)
1773  if(!cuts->GetChi2A()) cuts->SetChi2A(new TArrayD(9, chiA2));
1774  if(!cuts->GetChi2C()) cuts->SetChi2C(new TArrayD(9, chiC2));
1775  if(!cuts->GetChi3A()) cuts->SetChi3A(new TArrayD(9, chiA3));
1776  if(!cuts->GetChi3C()) cuts->SetChi3C(new TArrayD(9, chiC3));
1777 
1778  if(!fChi2A) fChi2A = cuts->GetChi2A();
1779  if(!fChi2C) fChi2C = cuts->GetChi2C();
1780  if(!fChi3A) fChi3A = cuts->GetChi3A();
1781  if(!fChi3C) fChi3C = cuts->GetChi3C();
1782 
1783  // get the calibration file for the gain equalization
1784  TFile *foadb = TFile::Open("alien:///alice/cern.ch/user/j/jmargutt/gainVZERO.LHC11h.root");
1785  if(!foadb){
1786  printf("file alien:///alice/cern.ch/user/j/jmargutt/gainVZERO.LHC11h.root cannot be opened, CALIBRATION FAILED !");
1787  return;
1788  }
1789  TH3F* Weights = dynamic_cast<TH3F*>(foadb->FindObjectAny("LHC11h")->FindObject("gHistVZEROChannelGainEqualizationMap"));
1790  if(!Weights){
1791  printf("gHistVZEROChannelGainEqualizationMap is not available in the file\n");
1792  return;
1793  }
1794 
1795  if(cuts->GetVZEROgainEqualizationPerRing()) {
1796  // enable or disable rings through the weights, weight 1. is enabled, 0. is disabled
1797  // start with the vzero c rings (segments 0 through 31)
1798  (cuts->GetUseVZERORing(0)) ? cuts->SetVZEROCpol(0, 1.) : cuts->SetVZEROCpol(0, 0.);
1799  (cuts->GetUseVZERORing(1)) ? cuts->SetVZEROCpol(1, 1.) : cuts->SetVZEROCpol(1, 0.);
1800  (cuts->GetUseVZERORing(2)) ? cuts->SetVZEROCpol(2, 1.) : cuts->SetVZEROCpol(2, 0.);
1801  (cuts->GetUseVZERORing(3)) ? cuts->SetVZEROCpol(3, 1.) : cuts->SetVZEROCpol(3, 0.);
1802  // same for vzero a
1803  (cuts->GetUseVZERORing(4)) ? cuts->SetVZEROApol(0, 1.) : cuts->SetVZEROApol(0, 0.);
1804  (cuts->GetUseVZERORing(5)) ? cuts->SetVZEROApol(1, 1.) : cuts->SetVZEROApol(1, 0.);
1805  (cuts->GetUseVZERORing(6)) ? cuts->SetVZEROApol(2, 1.) : cuts->SetVZEROApol(2, 0.);
1806  (cuts->GetUseVZERORing(7)) ? cuts->SetVZEROApol(3, 1.) : cuts->SetVZEROApol(3, 0.);
1807  } else {
1808  // else enable all rings, which is also default
1809  for(Int_t i(0); i < 4; i++) cuts->SetVZEROCpol(i, 1.);
1810  for(Int_t i(0); i < 4; i++) cuts->SetVZEROApol(i, 1.);
1811  }
1812 
1813  // loop over the 11h runs to see if it's 11h data (Full TPC Flow)
1814  Int_t runs11h[] = {167915, 168115, 168460, 169035, 169238, 169859, 170228, 167920, 168310, 168464, 169091, 169411, 169923, 170230, 167985, 168311, 168467, 169094, 169415, 170027, 170268, 167987, 168322, 168511, 169138, 169417, 170081, 170269, 167988, 168325, 168512, 169144, 169835, 170155, 170270, 168069, 168341, 168514, 169145, 169837, 170159, 170306, 168076, 168342, 168777, 169148, 169838, 170163, 170308, 168105, 168361, 168826, 169156, 169846, 170193, 170309, 168107, 168362, 168988, 169160, 169855, 170203, 168108, 168458, 168992, 169167, 169858, 170204};
1815  Int_t RunBin = -1;
1816  for(Int_t r(0); r < 68; r++) {
1817  if(run == runs11h[r]) {
1818  RunBin = r+1;
1819  }
1820  }
1821  if(RunBin>0) {
1822  printf(" > run has been identified as 11h Full TPC Flow < \n");
1823  } else {
1824  printf(" > run has NOT been identified as 11h Full TPC Flow, use default for 11h < \n");
1825  // pass a NULL pointer to the track cuts object, the NULL pointer will identify 11h runs
1826  cuts->SetVZEROgainEqualisation(NULL);
1827  fApplyRecentering = 2011;
1828  return; // the rest of the steps are not necessary
1829  }
1830 
1831  // step 0) get the TH2D which contains the correction factor per VZERO channel per centrality bin
1832  // and pass it to the current track cuts obect
1833  Weights->GetXaxis()->SetRange(RunBin,RunBin);
1834  TH2D* fMultVZEROCen = dynamic_cast<TH2D*>(Weights->Project3D("zy"));
1835  cuts->SetVZEROgainEqualisationCen(fMultVZEROCen); // passed as a TH2
1836 
1837  // get the calibration file for the q-vector recentering
1838  TFile *fqvec = TFile::Open("alien:///alice/cern.ch/user/j/jmargutt/recenteringVZERO.LHC11h.root");
1839  if(!fqvec){
1840  printf("file alien:///alice/cern.ch/user/j/jmargutt/recenteringVZERO.LHC11h.root cannot be opened, CALIBRATION FAILED !");
1841  return;
1842  }
1843 
1844  // first index of the oadb array is the harmonic n, the second index is either qax, qay, qcx, qcy
1845  TH2D* h[3][4];
1846  for(Int_t i(0); i < 3; i++) {
1847  h[i][0] = (TH2D*)(fqvec->FindObjectAny("LHC11h")->FindObject(Form("gHistVZEROAQ%ixRecenteringMap",i+1)));
1848  if(h[i][0]) fQxavsV0[i] = static_cast<TH1D*>(h[i][0]->ProjectionY(Form("fQxavsV0[%i]",i),RunBin,RunBin));
1849  h[i][1] = (TH2D*)(fqvec->FindObjectAny("LHC11h")->FindObject(Form("gHistVZEROAQ%iyRecenteringMap",i+1)));
1850  if(h[i][1]) fQyavsV0[i] = static_cast<TH1D*>(h[i][1]->ProjectionY(Form("fQyavsV0[%i]",i),RunBin,RunBin));
1851  h[i][2] = (TH2D*)(fqvec->FindObjectAny("LHC11h")->FindObject(Form("gHistVZEROCQ%ixRecenteringMap",i+1)));
1852  if(h[i][2]) fQxcvsV0[i] = static_cast<TH1D*>(h[i][2]->ProjectionY(Form("fQxcvsV0[%i]",i),RunBin,RunBin));
1853  h[i][3] = (TH2D*)(fqvec->FindObjectAny("LHC11h")->FindObject(Form("gHistVZEROCQ%iyRecenteringMap",i+1)));
1854  if(h[i][3]) fQycvsV0[i] = static_cast<TH1D*>(h[i][3]->ProjectionY(Form("fQycvsV0[%i]",i),RunBin,RunBin));
1855  }
1856 
1857  // set the recentering style (might be switched back to -1 if recentering is disabeled)
1858  // FIXME as an ugly hack, for now I mark this as 666 to denote the experimental nature
1859  // of this and use it transparently without disrupting the existing calbiration
1860  fApplyRecentering = 666;
1861 }
1862 
1863 //-----------------------------------------------------------------------------
1864 
1865 void AliFlowEvent::ClearFast()
1866 {
1867  //clear the event without releasing any memory
1868  //note that cached run number of recentering settigns are not clear
1869  //(see AliFlowEvent::ClearCachedRun() )
1870  AliFlowEventSimple::ClearFast();
1871 }
1872 
1873 //_____________________________________________________________________________
1874 
1875 void AliFlowEvent::ClearCachedRun()
1876 {
1877  //clear the cached run (not in clear fast as cache needs to be persistent in most cases )
1878  fCachedRun=0;
1879  fApplyRecentering=0;
1880 }
AliFlowEvent::GetQ
virtual AliFlowVector GetQ(Int_t n=2, TList *weightsList=NULL, Bool_t usePhiWeights=kFALSE, Bool_t usePtWeights=kFALSE, Bool_t useEtaWeights=kFALSE)
Definition: AliFlowEvent.cxx:1176
AliFlowTrack.h
GetPmdEta
Float_t GetPmdEta(Float_t xPos, Float_t yPos, Float_t zPos)
Definition: AliFlowEvent.cxx:1143
AliFlowEvent.h
AliFlowEvent::fChi2A
TArrayD * fChi2A
current event
Definition: AliFlowEvent.h:117
AliFlowEvent::ClearCachedRun
virtual void ClearCachedRun()
Definition: AliFlowEvent.cxx:1875
AliFlowEvent::FindDaughters
void FindDaughters(Bool_t keepDaughtersInRPselection=kFALSE)
Definition: AliFlowEvent.cxx:828
AliFlowEvent::SetVZEROCalibrationForTrackCuts
void SetVZEROCalibrationForTrackCuts(AliFlowTrackCuts *cuts)
Definition: AliFlowEvent.cxx:1390
AliFlowTrackCuts::SetChi3C
void SetChi3C(TArrayD *Chi3C)
Definition: AliFlowTrackCuts.h:328
AliFlowTrackSimple::SetDaughter
virtual void SetDaughter(Int_t, AliFlowTrackSimple *)
Definition: AliFlowTrackSimple.h:131
AliFlowEventSimple::GetQ
virtual AliFlowVector GetQ(Int_t n=2, TList *weightsList=NULL, Bool_t usePhiWeights=kFALSE, Bool_t usePtWeights=kFALSE, Bool_t useEtaWeights=kFALSE)
Definition: AliFlowEventSimple.cxx:461
AliFlowTrackSimple::SetEta
void SetEta(Double_t eta)
Definition: AliFlowTrackSimple.h:174
AliFlowTrack::Clear
virtual void Clear(Option_t *o="")
Definition: AliFlowTrack.h:42
AliFlowTrackCuts::GetInputObject
TObject * GetInputObject(Int_t i)
Definition: AliFlowTrackCuts.cxx:2672
AliFlowTrackSimple::GetNDaughters
virtual Int_t GetNDaughters() const
Definition: AliFlowTrackSimple.h:128
AliFlowTrackCuts.h
AliFlowTrackCuts::SetChi3A
void SetChi3A(TArrayD *Chi3A)
Definition: AliFlowTrackCuts.h:326
AliFlowEvent::fChi3C
TArrayD * fChi3C
chi vs cent for vzero A ep_3
Definition: AliFlowEvent.h:120
AliFlowEvent::fChi2C
TArrayD * fChi2C
chi vs cent for vzero A ep_2
Definition: AliFlowEvent.h:118
AliFlowVector::GetMult
Double_t GetMult() const
Definition: AliFlowVector.h:46
AliFlowEvent::SetMCReactionPlaneAngle
void SetMCReactionPlaneAngle(const AliMCEvent *mcEvent)
Definition: AliFlowEvent.cxx:173
AliFlowEventSimple::operator=
AliFlowEventSimple & operator=(const AliFlowEventSimple &anEvent)
Definition: AliFlowEventSimple.cxx:234
AliFlowTrackCuts::GetRun
Int_t GetRun() const
Definition: AliFlowTrackCuts.h:339
AliFlowEvent::operator=
AliFlowEvent & operator=(const AliFlowEvent &event)
Definition: AliFlowEvent.cxx:129
AliFlowTrackCuts::GetChi2A
TArrayD * GetChi2A()
Definition: AliFlowTrackCuts.h:330
AliFlowTrackCuts::SetChi2A
void SetChi2A(TArrayD *Chi2A)
Definition: AliFlowTrackCuts.h:325
AliFlowTrackSimple::GetIDDaughter
virtual Int_t GetIDDaughter(Int_t) const
Definition: AliFlowTrackSimple.h:130
AliFlowEventSimple::AddTrack
void AddTrack(AliFlowTrackSimple *track)
Definition: AliFlowEventSimple.cxx:420
AliFlowTrackCuts::IsSelected
virtual Bool_t IsSelected(TObject *obj, Int_t id=-666)
Definition: AliFlowTrackCuts.cxx:729
AliFlowEvent::fEvent
AliVEvent * fEvent
recentering
Definition: AliFlowEvent.h:116
AliFlowTrackCuts::SetVZEROgainEqualisationCen
void SetVZEROgainEqualisationCen(TH2 *g)
Definition: AliFlowTrackCuts.h:307
ClassImp
ClassImp(AliFlowEvent) AliFlowEvent
Definition: AliFlowEvent.cxx:55
AliFlowTrackCuts::GetChi2C
TArrayD * GetChi2C()
Definition: AliFlowTrackCuts.h:332
AliFlowEventSimple.h
AliFlowTrackCuts::GetParamType
trackParameterType GetParamType() const
Definition: AliFlowTrackCuts.h:226
AliFlowTrackSimple::TagRP
void TagRP(Bool_t b=kTRUE)
Definition: AliFlowTrackSimple.h:53
AliFlowEvent::fDivSigma
Bool_t fDivSigma
Definition: AliFlowEvent.h:103
AliFlowEventSimple::IncrementNumberOfPOIs
void IncrementNumberOfPOIs(Int_t poiType=1)
Definition: AliFlowEventSimple.cxx:226
AliFlowEvent::GetTrack
AliFlowTrack * GetTrack(Int_t i)
Definition: AliFlowEvent.cxx:164
AliFlowTrackCuts::GetVZEROgainEqualizationPerRing
Bool_t GetVZEROgainEqualizationPerRing() const
Definition: AliFlowTrackCuts.h:314
AliFlowTrackSimple::SetForRPSelection
void SetForRPSelection(Bool_t b=kTRUE)
Definition: AliFlowTrackSimple.h:188
AliFlowTrackCuts::GetApplyRecentering
Bool_t GetApplyRecentering() const
Definition: AliFlowTrackCuts.h:312
AliFlowTrackCuts::GetNumberOfInputObjects
Int_t GetNumberOfInputObjects() const
Definition: AliFlowTrackCuts.cxx:2622
AliFlowEvent::SetDeltaVZEROCalibrationForTrackCuts
void SetDeltaVZEROCalibrationForTrackCuts(AliFlowTrackCuts *cuts)
Definition: AliFlowEvent.cxx:1738
AliFlowEvent::Fill
void Fill(AliFlowTrackCuts *rpCuts, AliFlowTrackCuts *poiCuts)
Definition: AliFlowEvent.cxx:856
AliFlowTrackCuts::FillFlowTrack
Bool_t FillFlowTrack(AliFlowTrack *track) const
Definition: AliFlowTrackCuts.cxx:2130
AliFlowEvent::fQxcvsV0
TH1 * fQxcvsV0[5]
recentering
Definition: AliFlowEvent.h:113
AliFlowTrackCuts::SetVZEROCpol
void SetVZEROCpol(Int_t ring, Float_t f)
Definition: AliFlowTrackCuts.h:309
AliFlowTrackCuts::GetChi3A
TArrayD * GetChi3A()
Definition: AliFlowTrackCuts.h:331
AliFlowTrackCuts::GetDivSigma
Bool_t GetDivSigma() const
Definition: AliFlowTrackCuts.h:316
AliFlowEvent::fWidthQ
Float_t fWidthQ[9][2][2]
recentering
Definition: AliFlowEvent.h:107
AliFlowTrackCuts::SetVZEROgainEqualisation
void SetVZEROgainEqualisation(TH1 *g)
Definition: AliFlowTrackCuts.h:306
GetPmdPhi
Float_t GetPmdPhi(Float_t xPos, Float_t yPos)
Definition: AliFlowEvent.cxx:1152
AliFlowTrackSimple::Tag
void Tag(Int_t n, Bool_t b=kTRUE)
Definition: AliFlowTrackSimple.h:55
AliFlowTrackCuts::GetEvent
AliVEvent * GetEvent() const
Definition: AliFlowTrackCuts.h:245
AliFlowEventSimple::fTrackCollection
TObjArray * fTrackCollection
Definition: AliFlowEventSimple.h:147
AliFlowTrack::SetSource
void SetSource(trackSource s)
Definition: AliFlowTrack.h:37
AliFlowEvent::fMeanQv3
Float_t fMeanQv3[9][2][2]
recentering
Definition: AliFlowEvent.h:108
AliFlowEventSimple::Get2Qsub
virtual void Get2Qsub(AliFlowVector *Qarray, Int_t n=2, TList *weightsList=NULL, Bool_t usePhiWeights=kFALSE, Bool_t usePtWeights=kFALSE, Bool_t useEtaWeights=kFALSE)
Definition: AliFlowEventSimple.cxx:578
AliFlowTrackSimple::SetPhi
void SetPhi(Double_t phi)
Definition: AliFlowTrackSimple.h:178
AliFlowEvent::fVZEROcentralityBin
Int_t fVZEROcentralityBin
cached calibration info for vzero
Definition: AliFlowEvent.h:105
AliFlowTrackSimple::SetPt
void SetPt(Double_t pt)
Definition: AliFlowTrackSimple.h:176
AliFlowTrackCuts::GetUseVZERORing
Bool_t GetUseVZERORing(Int_t i) const
Definition: AliFlowTrackCuts.h:324
AliFlowTrackCuts::GetChi3C
TArrayD * GetChi3C()
Definition: AliFlowTrackCuts.h:333
AliFlowTrackSimple::SetWeight
void SetWeight(Double_t weight)
Definition: AliFlowTrackSimple.h:180
AliFlowEvent::SetBetaVZEROCalibrationForTrackCuts
void SetBetaVZEROCalibrationForTrackCuts(AliFlowTrackCuts *cuts)
Definition: AliFlowEvent.cxx:1591
AliFlowTrackSimple::SetForPOISelection
void SetForPOISelection(Bool_t b=kTRUE)
Definition: AliFlowTrackSimple.h:190
AliFlowEventSimple::fMothersCollection
TObjArray * fMothersCollection
Definition: AliFlowEventSimple.h:170
AliFlowEventSimple::SetMCReactionPlaneAngle
void SetMCReactionPlaneAngle(Double_t fPhiRP)
Definition: AliFlowEventSimple.h:81
AliFlowEvent::fMeanQ
Float_t fMeanQ[9][2][2]
centrality bin for the current event
Definition: AliFlowEvent.h:106
AliFlowVector.h
AliFlowEvent::fQycvsV0
TH1 * fQycvsV0[5]
recentering
Definition: AliFlowEvent.h:114
AliFlowEvent::ReuseTrack
AliFlowTrack * ReuseTrack(Int_t i)
Definition: AliFlowEvent.cxx:1005
AliFlowEvent::Get2Qsub
virtual void Get2Qsub(AliFlowVector *Qarray, Int_t n=2, TList *weightsList=0x0, Bool_t usePhiWeights=0x0, Bool_t usePtWeights=0x0, Bool_t useEtaWeights=0x0)
Definition: AliFlowEvent.cxx:1219
AliFlowEvent::fQyavsV0
TH1 * fQyavsV0[5]
recentering
Definition: AliFlowEvent.h:112
AliFlowEvent::fCachedRun
Int_t fCachedRun
Definition: AliFlowEvent.h:104
AliFlowEvent::InsertTrack
void InsertTrack(AliFlowTrack *)
Definition: AliFlowEvent.cxx:993
AliFlowEvent::ClearFast
virtual void ClearFast()
Definition: AliFlowEvent.cxx:1865
AliFlowTrackCuts::SetChi2C
void SetChi2C(TArrayD *Chi2C)
Definition: AliFlowTrackCuts.h:327
AliFlowEvent::fQxavsV0
TH1 * fQxavsV0[5]
recentering
Definition: AliFlowEvent.h:111
AliFlowEvent::fChi3A
TArrayD * fChi3A
chi vs cent for vzero C ep_2
Definition: AliFlowEvent.h:119
AliFlowEvent::fWidthQv3
Float_t fWidthQv3[9][2][2]
recentering
Definition: AliFlowEvent.h:109
AliFlowTrackCuts::SetVZEROApol
void SetVZEROApol(Int_t ring, Float_t f)
Definition: AliFlowTrackCuts.h:308
AliFlowEvent::fApplyRecentering
Int_t fApplyRecentering
Definition: AliFlowEvent.h:102