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