more efficient use of vzero calibration code, remove the false impression that the...
[u/mrichter/AliRoot.git] / PWG / FLOW / Tasks / AliFlowEvent.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /*****************************************************************
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 "TH2F.h"
28 #include "TProfile.h"
29 #include "AliMCEvent.h"
30 #include "AliMCParticle.h"
31 #include "AliCFManager.h"
32 #include "AliESDtrack.h"
33 #include "AliESDPmdTrack.h"
34 #include "AliESDEvent.h"
35 #include "AliAODEvent.h"
36 #include "AliOADBContainer.h"
37 #include "AliGenCocktailEventHeader.h"
38 #include "AliGenEposEventHeader.h"
39 #include "AliGenHijingEventHeader.h"
40 #include "AliGenGeVSimEventHeader.h"
41 #include "AliCollisionGeometry.h"
42 #include "AliMultiplicity.h"
43 #include "AliFlowTrackCuts.h"
44 #include "AliFlowEventSimple.h"
45 #include "AliFlowTrack.h"
46 #include "AliFlowVector.h"
47 #include "AliFlowEvent.h"
48 #include "AliLog.h"
49
50 using std::endl;
51 using std::cout;
52 ClassImp(AliFlowEvent)
53
54 //-----------------------------------------------------------------------
55 AliFlowEvent::AliFlowEvent():
56   AliFlowEventSimple(), fApplyRecentering(-1), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0)
57 {
58     // constructor
59     for(Int_t i(0); i < 9; i++) {
60         for(Int_t j(0); j < 2; j++) {
61             for(Int_t k(0); k < 2; k++) {
62                 fMeanQ[i][j][k] = 0.; 
63                 fWidthQ[i][j][k] = 0.;  
64                 fMeanQv3[i][j][k] = 0.; 
65                 fWidthQv3[i][j][k] = 0.;
66             }
67         }
68     }
69   //ctor
70   cout << "AliFlowEvent: Default constructor to be used only by root for io" << endl;
71 }
72
73 //-----------------------------------------------------------------------
74 AliFlowEvent::AliFlowEvent(Int_t n):
75   AliFlowEventSimple(n), fApplyRecentering(-1), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0)
76 {
77     // constructor
78     for(Int_t i(0); i < 9; i++) {
79         for(Int_t j(0); j < 2; j++) {
80             for(Int_t k(0); k < 2; k++) {
81                 fMeanQ[i][j][k] = 0.; 
82                 fWidthQ[i][j][k] = 0.;  
83                 fMeanQv3[i][j][k] = 0.; 
84                 fWidthQv3[i][j][k] = 0.;
85             }
86         }
87     }
88 }
89
90 //-----------------------------------------------------------------------
91 AliFlowEvent::AliFlowEvent(const AliFlowEvent& event):
92   AliFlowEventSimple(event), fApplyRecentering(event.fApplyRecentering), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0)
93 {
94     // copy constructor 
95     for(Int_t i(0); i < 9; i++) {
96         for(Int_t j(0); j < 2; j++) {
97             for(Int_t k(0); k < 2; k++) {
98                 fMeanQ[i][j][k] = 0.; 
99                 fWidthQ[i][j][k] = 0.;  
100                 fMeanQv3[i][j][k] = 0.; 
101                 fWidthQv3[i][j][k] = 0.;
102             }
103         }
104     }
105 }
106
107 //-----------------------------------------------------------------------
108 AliFlowEvent& AliFlowEvent::operator=(const AliFlowEvent& event)
109 {
110   //assignment operator
111   if (&event==this) return *this;       // check self-assignment
112
113   fApplyRecentering = event.fApplyRecentering; 
114   fCachedRun = event.fCachedRun; 
115   fVZEROcentralityBin = event.fVZEROcentralityBin;
116   fEvent = 0x0;                         // should never be copied
117   for(Int_t i(0); i < 9; i++) {
118       for(Int_t j(0); j < 2; j++) {
119           for(Int_t k(0); k < 2; k++) {
120               fMeanQ[i][j][k] = event.fMeanQ[i][j][k]; 
121               fWidthQ[i][j][k] = event.fWidthQ[i][j][k];  
122               fMeanQv3[i][j][k] = event.fMeanQv3[i][j][k]; 
123               fWidthQv3[i][j][k] = event.fWidthQv3[i][j][k];
124           }
125       }
126   }
127   AliFlowEventSimple::operator=(event);
128   return *this;
129 }
130
131 //-----------------------------------------------------------------------
132 AliFlowTrack* AliFlowEvent::GetTrack(Int_t i)
133 {
134   //get track i from collection
135   if (i>=fNumberOfTracks) return NULL;
136   AliFlowTrack* pTrack = static_cast<AliFlowTrack*>(fTrackCollection->At(i)) ;
137   return pTrack;
138 }
139
140 //-----------------------------------------------------------------------
141 void AliFlowEvent::SetMCReactionPlaneAngle(const AliMCEvent* mcEvent)
142 {
143   //sets the event plane angle from the proper header in the MC
144
145   //COCKTAIL with HIJING
146   if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"Cocktail Header"))   //returns 0 if matches
147   {
148     AliGenCocktailEventHeader *headerC = dynamic_cast<AliGenCocktailEventHeader *> (mcEvent-> GenEventHeader());
149     if (headerC)
150     {
151       TList *lhd = headerC->GetHeaders();
152       if (lhd)
153       {
154         AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (lhd->At(0));
155         if (hdh) AliFlowEventSimple::SetMCReactionPlaneAngle( hdh->ReactionPlaneAngle() );
156       }
157     }
158   }
159   //THERMINATOR
160   else if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"Therminator"))   //returns 0 if matches
161   {
162     AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
163     if (headerH) AliFlowEventSimple::SetMCReactionPlaneAngle( headerH->ReactionPlaneAngle() );
164   }
165   //GEVSIM
166   else if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"GeVSim header"))   //returns 0 if matches
167   {
168     AliGenGeVSimEventHeader* headerG = dynamic_cast<AliGenGeVSimEventHeader*>(mcEvent->GenEventHeader());
169     if (headerG) AliFlowEventSimple::SetMCReactionPlaneAngle( headerG->GetEventPlane() );
170   }
171   //HIJING
172   else if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"Hijing"))   //returns 0 if matches
173   {
174     AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
175     if (headerH) AliFlowEventSimple::SetMCReactionPlaneAngle( headerH->ReactionPlaneAngle() );
176   }
177   //AMPT
178   else if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"Ampt"))   //returns 0 if matches
179   {
180     AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
181     if (headerH) AliFlowEventSimple::SetMCReactionPlaneAngle( headerH->ReactionPlaneAngle() );
182   }
183   //EPOS
184   else if (!strcmp(mcEvent->GenEventHeader()->GetName(),"EPOS"))
185   {
186     AliGenEposEventHeader* headerE = dynamic_cast<AliGenEposEventHeader*>(mcEvent->GenEventHeader());
187     if (headerE) AliFlowEventSimple::SetMCReactionPlaneAngle( headerE->ReactionPlaneAngle() );
188   }
189   //Hydjet
190   else
191   {
192     AliCollisionGeometry* header = dynamic_cast<AliCollisionGeometry*>(mcEvent->GenEventHeader());
193     if (header) AliFlowEventSimple::SetMCReactionPlaneAngle( header->ReactionPlaneAngle() );
194   }
195 }
196
197 //-----------------------------------------------------------------------
198 AliFlowEvent::AliFlowEvent( const AliMCEvent* anInput,
199                             const AliCFManager* rpCFManager,
200                             const AliCFManager* poiCFManager):
201   AliFlowEventSimple(20), fApplyRecentering(-1), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0)
202 {
203     // constructor
204     for(Int_t i(0); i < 9; i++) {
205         for(Int_t j(0); j < 2; j++) {
206             for(Int_t k(0); k < 2; k++) {
207                 fMeanQ[i][j][k] = 0.; 
208                 fWidthQ[i][j][k] = 0.;  
209                 fMeanQv3[i][j][k] = 0.; 
210                 fWidthQv3[i][j][k] = 0.;
211             }
212         }
213     }
214   //Fills the event from the MC kinematic information
215
216   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
217
218   //loop over tracks
219   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
220   {
221     //get input particle
222     AliMCParticle* pParticle = dynamic_cast<AliMCParticle*>(anInput->GetTrack(itrkN));
223     if (!pParticle) continue;
224
225     //check if pParticle passes the cuts
226     Bool_t rpOK = kTRUE;
227     Bool_t poiOK = kTRUE;
228     if (rpCFManager && poiCFManager)
229     {
230       rpOK = rpCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pParticle);
231       poiOK = poiCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pParticle);
232     }
233     if (!(rpOK||poiOK)) continue;
234
235     AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
236     pTrack->SetSource(AliFlowTrack::kFromMC);
237
238     if (rpOK && rpCFManager)
239     {
240       pTrack->SetForRPSelection(kTRUE);
241       IncrementNumberOfPOIs(0);
242     }
243     if (poiOK && poiCFManager)
244     {
245       pTrack->SetForPOISelection(kTRUE);
246       IncrementNumberOfPOIs(1);
247     }
248
249     AddTrack(pTrack) ;
250   }//for all tracks
251   SetMCReactionPlaneAngle(anInput);
252 }
253
254 //-----------------------------------------------------------------------
255 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
256                             const AliCFManager* rpCFManager,
257                             const AliCFManager* poiCFManager ):
258   AliFlowEventSimple(20),  fApplyRecentering(-1), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0)
259 {
260     // constructor
261     for(Int_t i(0); i < 9; i++) {
262         for(Int_t j(0); j < 2; j++) {
263             for(Int_t k(0); k < 2; k++) {
264                 fMeanQ[i][j][k] = 0.; 
265                 fWidthQ[i][j][k] = 0.;  
266                 fMeanQv3[i][j][k] = 0.; 
267                 fWidthQv3[i][j][k] = 0.;
268             }
269         }
270     }
271    
272   //Fills the event from the ESD
273
274   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
275
276   //loop over tracks
277   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
278   {
279     AliESDtrack* pParticle = anInput->GetTrack(itrkN);   //get input particle
280
281     //check if pParticle passes the cuts
282     Bool_t rpOK = kTRUE;
283     Bool_t poiOK = kTRUE;
284     if (rpCFManager && poiCFManager)
285     {
286       rpOK = ( rpCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
287                rpCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
288       poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
289                 poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
290     }
291     if (!(rpOK || poiOK)) continue;
292
293     //make new AliFLowTrack
294     AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
295     pTrack->SetSource(AliFlowTrack::kFromESD);
296
297     //marking the particles used for int. flow:
298     if(rpOK && rpCFManager)
299     {
300       pTrack->SetForRPSelection(kTRUE);
301       IncrementNumberOfPOIs(0);
302     }
303     //marking the particles used for diff. flow:
304     if(poiOK && poiCFManager)
305     {
306       pTrack->SetForPOISelection(kTRUE);
307       IncrementNumberOfPOIs(1);
308     }
309
310     AddTrack(pTrack);
311   }//end of while (itrkN < iNumberOfInputTracks)
312 }
313
314 //-----------------------------------------------------------------------
315 AliFlowEvent::AliFlowEvent( const AliAODEvent* anInput,
316                             const AliCFManager* rpCFManager,
317                             const AliCFManager* poiCFManager):
318   AliFlowEventSimple(20), fApplyRecentering(-1), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0)
319 {
320     // constructor
321     for(Int_t i(0); i < 9; i++) {
322         for(Int_t j(0); j < 2; j++) {
323             for(Int_t k(0); k < 2; k++) {
324                 fMeanQ[i][j][k] = 0.; 
325                 fWidthQ[i][j][k] = 0.;  
326                 fMeanQv3[i][j][k] = 0.; 
327                 fWidthQv3[i][j][k] = 0.;
328             }
329         }
330     }
331   //Fills the event from the AOD
332   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
333
334   //loop over tracks
335   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
336   {
337     AliAODTrack* pParticle = anInput->GetTrack(itrkN);   //get input particle
338
339     //check if pParticle passes the cuts
340     Bool_t rpOK = kTRUE;
341     Bool_t poiOK = kTRUE;
342     if (rpCFManager && poiCFManager)
343     {
344       rpOK = ( rpCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
345                rpCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
346       poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
347                 poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
348     }
349     if (!(rpOK || poiOK)) continue;
350
351     //make new AliFlowTrack
352     AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
353     pTrack->SetSource(AliFlowTrack::kFromAOD);
354
355     if (rpOK /* && rpCFManager */ ) // to be fixed - with CF managers uncommented only empty events (NULL in header files)
356     {
357       pTrack->SetForRPSelection(kTRUE);
358       IncrementNumberOfPOIs(0);
359     }
360     if (poiOK /* && poiCFManager*/ )
361     {
362       pTrack->SetForPOISelection(kTRUE);
363       IncrementNumberOfPOIs(1);
364     }
365     AddTrack(pTrack);
366   }
367
368   //  if (iSelParticlesRP >= fMinMult && iSelParticlesRP <= fMaxMult)
369   //  {
370   //    if ( (++fCount % 100) == 0)
371   //    {
372   //      if (!fMCReactionPlaneAngle == 0) cout<<" MC Reaction Plane Angle = "<<  fMCReactionPlaneAngle << endl;
373   //      else cout<<" MC Reaction Plane Angle = unknown "<< endl;
374   //      cout<<" iGoodTracks = "<<iGoodTracks<<endl;
375   //      cout<<" # of RP selected tracks = "<<iSelParticlesRP<<endl;
376   //      cout<<" # of POI selected tracks = "<<iSelParticlesPOI<<endl;
377   //      cout << "# " << fCount << " events processed" << endl;
378   //    }
379   //    return pEvent;
380   //  }
381   //  else
382   //  {
383   //    cout<<"Not enough tracks in the FlowEventSimple"<<endl;
384   //    return 0;
385   //  }
386   //}
387   //else
388   //{
389   //  cout<<"Event does not pass multiplicity cuts"<<endl;
390   //  return 0;
391   //}
392
393 }
394
395 //-----------------------------------------------------------------------
396 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
397                             const AliMCEvent* anInputMc,
398                             KineSource anOption,
399                             const AliCFManager* rpCFManager,
400                             const AliCFManager* poiCFManager ):
401   AliFlowEventSimple(20), fApplyRecentering(-1), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0)
402 {
403     // constructor
404     for(Int_t i(0); i < 9; i++) {
405         for(Int_t j(0); j < 2; j++) {
406             for(Int_t k(0); k < 2; k++) {
407                 fMeanQ[i][j][k] = 0.; 
408                 fWidthQ[i][j][k] = 0.;  
409                 fMeanQv3[i][j][k] = 0.; 
410                 fWidthQv3[i][j][k] = 0.;
411             }
412         }
413     }
414   //fills the event with tracks from the ESD and kinematics from the MC info via the track label
415   if (anOption==kNoKine)
416   {
417     AliFatal("WRONG OPTION IN AliFlowEventMaker::FillTracks(AliESDEvent* anInput, AliMCEvent* anInputMc, KineSource anOption)");
418     exit(1);
419   }
420
421   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
422
423   Int_t iNumberOfInputTracksMC = anInputMc->GetNumberOfTracks() ;
424   if (iNumberOfInputTracksMC==-1)
425   {
426     AliError("Skipping Event -- No MC information available for this event");
427     return;
428   }
429
430   //loop over ESD tracks
431   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
432   {
433     AliESDtrack* pParticle = anInput->GetTrack(itrkN);   //get input particle
434     //get Label
435     Int_t iLabel = pParticle->GetLabel();
436     //match to mc particle
437     AliMCParticle* pMcParticle = (AliMCParticle*) anInputMc->GetTrack(TMath::Abs(iLabel));
438
439     //check
440     if (TMath::Abs(pParticle->GetLabel())!=pMcParticle->Label())
441       AliWarning(Form("pParticle->GetLabel()!=pMcParticle->Label(), %i, %i", pParticle->GetLabel(), pMcParticle->Label()));
442
443     //check if pParticle passes the cuts
444     Bool_t rpOK = kFALSE;
445     Bool_t poiOK = kFALSE;
446     if (rpCFManager && poiCFManager)
447     {
448       if(anOption == kESDkine)
449       {
450         if (rpCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pMcParticle,"mcGenCuts1") &&
451             rpCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle))
452           rpOK=kTRUE;
453         if (poiCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pMcParticle,"mcGenCuts2") &&
454             poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle))
455           poiOK=kTRUE;
456       }
457       else if (anOption == kMCkine)
458       {
459         if (rpCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pMcParticle))
460           rpOK=kTRUE;
461         if (poiCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pMcParticle))
462           poiOK=kTRUE;
463       }
464     }
465
466     if (!(rpOK || poiOK)) continue;
467
468     //make new AliFlowTrack
469     AliFlowTrack* pTrack = NULL;
470     if(anOption == kESDkine)   //take the PID from the MC & the kinematics from the ESD
471     {
472       pTrack = new AliFlowTrack(pParticle);
473     }
474     else if (anOption == kMCkine)   //take the PID and kinematics from the MC
475     {
476       pTrack = new AliFlowTrack(pMcParticle);
477     }
478
479     if (rpOK && rpCFManager)
480     {
481       IncrementNumberOfPOIs(0);
482       pTrack->SetForRPSelection();
483     }
484     if (poiOK && poiCFManager) 
485     { 
486       IncrementNumberOfPOIs(1);
487       pTrack->SetForPOISelection();
488     }
489
490     AddTrack(pTrack);
491   }
492   SetMCReactionPlaneAngle(anInputMc);
493 }
494
495 //-----------------------------------------------------------------------
496 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
497                             const AliMultiplicity* anInputTracklets,
498                             const AliCFManager* poiCFManager ):
499   AliFlowEventSimple(20), fApplyRecentering(-1), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0)
500 {
501     // constructor
502     for(Int_t i(0); i < 9; i++) {
503         for(Int_t j(0); j < 2; j++) {
504             for(Int_t k(0); k < 2; k++) {
505                 fMeanQ[i][j][k] = 0.; 
506                 fWidthQ[i][j][k] = 0.;  
507                 fMeanQv3[i][j][k] = 0.; 
508                 fWidthQv3[i][j][k] = 0.;
509             }
510         }
511     }
512
513   //Select the particles of interest from the ESD
514   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
515
516   //loop over tracks
517   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
518     {
519       AliESDtrack* pParticle = anInput->GetTrack(itrkN);   //get input particle
520
521       //check if pParticle passes the cuts
522       Bool_t poiOK = kTRUE;
523       if (poiCFManager)
524         {
525           poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
526                     poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
527         }
528       if (!poiOK) continue;
529       
530       //make new AliFLowTrack
531       AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
532           
533       //marking the particles used for the particle of interest (POI) selection:
534       if(poiOK && poiCFManager)
535         {
536           IncrementNumberOfPOIs(1);
537           pTrack->SetForPOISelection(kTRUE);
538           pTrack->SetSource(AliFlowTrack::kFromESD);
539         }
540
541       AddTrack(pTrack);
542     }//end of while (itrkN < iNumberOfInputTracks)
543
544   //Select the reference particles from the SPD tracklets
545   anInputTracklets = anInput->GetMultiplicity();
546   Int_t multSPD = anInputTracklets->GetNumberOfTracklets();
547   
548   //loop over tracklets
549   for (Int_t itracklet=0; itracklet<multSPD; ++itracklet) {
550     Float_t thetaTr= anInputTracklets->GetTheta(itracklet);
551     Float_t phiTr= anInputTracklets->GetPhi(itracklet);
552     // calculate eta
553     Float_t etaTr = -TMath::Log(TMath::Tan(thetaTr/2.));
554     
555     //make new AliFLowTrackSimple
556     AliFlowTrack* pTrack = new AliFlowTrack();
557     pTrack->SetPt(0.0);
558     pTrack->SetEta(etaTr);
559     pTrack->SetPhi(phiTr);
560     //marking the particles used for the reference particle (RP) selection:
561     IncrementNumberOfPOIs(0);
562     pTrack->SetForRPSelection(kTRUE);
563     pTrack->SetSource(AliFlowTrack::kFromTracklet);
564
565     //Add the track to the flowevent
566     AddTrack(pTrack);
567   }
568
569 }
570
571 //-----------------------------------------------------------------------
572 AliFlowEvent::AliFlowEvent( const AliESDEvent* esd,
573                             const AliCFManager* poiCFManager,
574                             Bool_t hybrid):
575   AliFlowEventSimple(20), fApplyRecentering(-1), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0)
576 {
577     // constructor
578     for(Int_t i(0); i < 9; i++) {
579         for(Int_t j(0); j < 2; j++) {
580             for(Int_t k(0); k < 2; k++) {
581                 fMeanQ[i][j][k] = 0.; 
582                 fWidthQ[i][j][k] = 0.;  
583                 fMeanQv3[i][j][k] = 0.; 
584                 fWidthQv3[i][j][k] = 0.;
585             }
586         }
587     }
588
589   //Select the particles of interest from the ESD
590   Int_t iNumberOfInputTracks = esd->GetNumberOfTracks() ;
591
592   //Double_t gPt = 0.0, gP = 0.0;
593   Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};  //The impact parameters and their covariance.
594 //  Double_t dca3D = 0.0;       FIXME unused variable
595
596   AliESDtrack trackTPC;
597
598   //loop over tracks
599   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
600     {
601
602       if (!esd->GetTrack(itrkN)) continue;
603
604       Bool_t useTPC = kFALSE;
605
606       AliESDtrack* pParticle = esd->GetTrack(itrkN);   //get input particle
607
608       //check if pParticle passes the cuts
609       Bool_t poiOK = kTRUE;
610
611       if (poiCFManager)
612       {
613         poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
614                   poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
615       }
616
617       if (!(poiOK)) continue;
618
619       AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)pParticle->GetTPCInnerParam();
620
621       if (tpcTrack)
622       {
623
624 //      gPt = tpcTrack->Pt();
625 //      gP = tpcTrack->P();
626
627         useTPC = kTRUE;
628
629         const AliESDVertex *vertexSPD = esd->GetPrimaryVertexSPD();
630         const AliESDVertex *vertexTPC = esd->GetPrimaryVertexTPC();
631
632         AliExternalTrackParam copy(*tpcTrack);
633         if(hybrid)
634           copy.PropagateToDCA(vertexSPD,esd->GetMagneticField(),100.,dca,cov);
635         else
636           copy.PropagateToDCA(vertexTPC,esd->GetMagneticField(),100.,dca,cov);
637
638 //        dca3D = TMath::Sqrt(TMath::Power(dca[0],2)+TMath::Power(dca[1],2));   FIXME unused variable
639
640       }
641
642       //make new AliFLowTrack
643       AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
644
645       pTrack->SetSource(AliFlowTrack::kFromESD);
646
647       //marking the particles used for diff. flow:
648       if(poiOK && poiCFManager)
649       {
650         pTrack->SetForPOISelection(kTRUE);
651         IncrementNumberOfPOIs(1);
652       }
653
654       if(useTPC)
655       {
656         pTrack->SetForRPSelection(kTRUE);
657         IncrementNumberOfPOIs(0);
658       }
659
660       AddTrack(pTrack);
661
662     }//end of while (itrkN < iNumberOfInputTracks)
663
664 }
665
666 //-----------------------------------------------------------------------
667 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
668                             const TH2F* anInputFMDhist,
669                             const AliCFManager* poiCFManager ):
670   AliFlowEventSimple(20), fApplyRecentering(-1), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0)
671 {
672     // constructor
673     for(Int_t i(0); i < 9; i++) {
674         for(Int_t j(0); j < 2; j++) {
675             for(Int_t k(0); k < 2; k++) {
676                 fMeanQ[i][j][k] = 0.; 
677                 fWidthQ[i][j][k] = 0.;  
678                 fMeanQv3[i][j][k] = 0.; 
679                 fWidthQv3[i][j][k] = 0.;
680             }
681         }
682     }
683
684   //Select the particles of interest from the ESD
685   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
686
687   //loop over tracks
688   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
689     {
690       AliESDtrack* pParticle = anInput->GetTrack(itrkN);   //get input particle
691
692       //check if pParticle passes the cuts
693       Bool_t poiOK = kTRUE;
694       if (poiCFManager)
695         {
696           poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
697                     poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
698         }
699       if (!poiOK) continue;
700  
701       //make new AliFLowTrack
702       AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
703           
704       //marking the particles used for the particle of interest (POI) selection:
705       if(poiOK && poiCFManager)
706         {
707     IncrementNumberOfPOIs(1);
708           pTrack->SetForPOISelection(kTRUE);
709           pTrack->SetSource(AliFlowTrack::kFromESD);
710         }
711
712       AddTrack(pTrack);
713     }//end of while (itrkN < iNumberOfInputTracks)
714
715   //Select the reference particles from the FMD hits
716   //loop over FMD histogram
717   Int_t iBinsEta = anInputFMDhist->GetNbinsX();
718   Int_t iBinsPhi = anInputFMDhist->GetNbinsY();
719   
720   for (Int_t iEta = 1; iEta <= iBinsEta; iEta++){
721     Double_t etaFMD = anInputFMDhist->GetXaxis()->GetBinCenter(iEta);
722     for (Int_t iPhi = 1; iPhi <= iBinsPhi; iPhi++){
723       Double_t phiFMD = anInputFMDhist->GetYaxis()->GetBinCenter(iPhi);
724       Double_t weightFMD = anInputFMDhist->GetBinContent(iEta,iPhi);
725     
726       if (weightFMD > 0.0) { //do not add empty bins
727         //make new AliFLowTrackSimple
728         AliFlowTrack* pTrack = new AliFlowTrack();
729         pTrack->SetPt(0.0);
730         pTrack->SetEta(etaFMD);
731         pTrack->SetPhi(phiFMD);
732         pTrack->SetWeight(weightFMD);
733         //marking the particles used for the reference particle (RP) selection:
734         pTrack->TagRP();
735         IncrementNumberOfPOIs(0);
736         pTrack->SetSource(AliFlowTrack::kFromFMD);
737
738         //Add the track to the flowevent
739         AddTrack(pTrack);
740         
741       }
742     }
743   }
744
745 }
746
747 //-----------------------------------------------------------------------
748 void AliFlowEvent::FindDaughters(Bool_t keepDaughtersInRPselection)
749 {
750   //each flow track holds it's esd track index as well as its daughters esd index.
751   //fill the array of daughters for every track with the pointers to flow tracks
752   //to associate the mothers with daughters directly
753   for (Int_t iTrack=0; iTrack<fMothersCollection->GetEntriesFast(); iTrack++)
754   {
755     AliFlowTrack* mother = static_cast<AliFlowTrack*>(fMothersCollection->At(iTrack));
756     if (!mother) continue;
757     if (mother->GetNDaughters()<1) continue;
758     for (Int_t iDaughterCandidate=0; iDaughterCandidate<fNumberOfTracks; iDaughterCandidate++)
759     {
760       AliFlowTrack* daughterCandidate = static_cast<AliFlowTrack*>(fTrackCollection->At(iDaughterCandidate));
761       Int_t esdIndexDaughterCandidate = daughterCandidate->GetID();
762       for (Int_t iDaughter=0; iDaughter<mother->GetNDaughters(); iDaughter++)
763       {
764         Int_t esdIndexDaughter = mother->GetIDDaughter(iDaughter);
765         if (esdIndexDaughter==esdIndexDaughterCandidate)
766         {
767           mother->SetDaughter(iDaughter,daughterCandidate);
768           daughterCandidate->SetForRPSelection(keepDaughtersInRPselection);
769         }
770       }
771     }
772   }
773 }
774
775 //-----------------------------------------------------------------------
776 void AliFlowEvent::Fill( AliFlowTrackCuts* rpCuts,
777                          AliFlowTrackCuts* poiCuts )
778 {
779   //Fills the event from a vevent: AliESDEvent,AliAODEvent,AliMCEvent
780   //the input data needs to be attached to the cuts
781   //we have two cases, if we're cutting the same collection of tracks
782   //(same param type) then we can have tracks that are both rp and poi
783   //in the other case we want to have two exclusive sets of rps and pois
784   //e.g. one tracklets, the other PMD or global - USER IS RESPOSIBLE
785   //FOR MAKING SURE THEY DONT OVERLAP OR ELSE THE SAME PARTICLE WILL BE
786   //TAKEN TWICE
787
788   ClearFast();
789   if (!rpCuts || !poiCuts) return;
790   AliFlowTrackCuts::trackParameterType sourceRP = rpCuts->GetParamType();
791   AliFlowTrackCuts::trackParameterType sourcePOI = poiCuts->GetParamType();
792   AliFlowTrack* pTrack=NULL;
793  
794   // if the source for rp's or poi's is the VZERO detector, get the calibration 
795   // and set the calibration parameters
796   if (sourceRP == AliFlowTrackCuts::kVZERO) {
797       SetVZEROCalibrationForTrackCuts(rpCuts);
798       if(!rpCuts->GetApplyRecentering()) {
799           // if the user does not want to recenter, switch the flag
800           fApplyRecentering = -1;
801       }
802       // note: this flag is used in the overloaded implementation of Get2Qsub()
803       // and tells the function to use as Qsub vectors the recentered Q-vectors
804       // from the VZERO oadb file or from the event header
805   }
806   if (sourcePOI == AliFlowTrackCuts::kVZERO) {
807       // probably no-one will choose vzero tracks as poi's ...
808       SetVZEROCalibrationForTrackCuts(poiCuts); 
809   }
810   
811
812   if (sourceRP==sourcePOI)
813   {
814     //loop over tracks
815     Int_t numberOfInputObjects = rpCuts->GetNumberOfInputObjects();
816     for (Int_t i=0; i<numberOfInputObjects; i++)
817     {
818       //get input object (particle)
819       TObject* particle = rpCuts->GetInputObject(i);
820
821       Bool_t rp = rpCuts->IsSelected(particle,i);
822       Bool_t poi = poiCuts->IsSelected(particle,i);
823
824       if (!(rp||poi)) continue;
825
826       //make new AliFlowTrack
827       if (rp)
828       {
829         pTrack = rpCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
830         if (!pTrack) continue;
831         pTrack->Tag(0); IncrementNumberOfPOIs(0);
832         if (poi) {pTrack->Tag(1); IncrementNumberOfPOIs(1);}
833         if (pTrack->GetNDaughters()>0) fMothersCollection->Add(pTrack);
834       }
835       else if (poi)
836       {
837         pTrack = poiCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
838         if (!pTrack) continue;
839         pTrack->Tag(1); IncrementNumberOfPOIs(1);
840         if (pTrack->GetNDaughters()>0) fMothersCollection->Add(pTrack);
841       }
842       fNumberOfTracks++;
843     }//end of while (i < numberOfTracks)
844   }
845   else if (sourceRP!=sourcePOI)
846   {
847     //here we have two different sources of particles, so we fill
848     //them independently
849     //POI
850     for (Int_t i=0; i<poiCuts->GetNumberOfInputObjects(); i++)
851     {
852       TObject* particle = poiCuts->GetInputObject(i);
853       Bool_t poi = poiCuts->IsSelected(particle,i);
854       if (!poi) continue;
855       pTrack = poiCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
856       if (!pTrack) continue;
857       pTrack->Tag(1);
858       IncrementNumberOfPOIs(1);
859       fNumberOfTracks++;
860       if (pTrack->GetNDaughters()>0) fMothersCollection->Add(pTrack);
861     }
862     //RP
863     Int_t numberOfInputObjects = rpCuts->GetNumberOfInputObjects();
864     for (Int_t i=0; i<numberOfInputObjects; i++)
865       {
866       TObject* particle = rpCuts->GetInputObject(i);
867       Bool_t rp = rpCuts->IsSelected(particle,i);
868       if (!rp) continue;
869       pTrack = rpCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
870       if (!pTrack) continue;
871       pTrack->Tag(0);
872       IncrementNumberOfPOIs(0);
873       fNumberOfTracks++;
874       if (pTrack->GetNDaughters()>0) fMothersCollection->Add(pTrack);
875     }
876   }
877 }
878
879 //-----------------------------------------------------------------------
880 void AliFlowEvent::InsertTrack(AliFlowTrack *track) {
881   // adds a flow track at the end of the container
882   AliFlowTrack *pTrack = ReuseTrack( fNumberOfTracks++ );
883   *pTrack = *track;
884   if (track->GetNDaughters()>0)
885   {
886     fMothersCollection->Add(pTrack);
887   }
888   return;
889 }
890
891 //-----------------------------------------------------------------------
892 AliFlowTrack* AliFlowEvent::ReuseTrack(Int_t i)
893 {
894   //try to reuse an existing track, if empty, make new one
895   AliFlowTrack* pTrack = static_cast<AliFlowTrack*>(fTrackCollection->At(i));
896   if (pTrack)
897   {
898     pTrack->Clear();
899   }
900   else 
901   {
902     pTrack = new AliFlowTrack();
903     fTrackCollection->AddAtAndExpand(pTrack,i);
904   }
905   return pTrack;
906 }
907
908 //-----------------------------------------------------------------------
909 AliFlowEvent::AliFlowEvent( AliFlowTrackCuts* rpCuts,
910                             AliFlowTrackCuts* poiCuts ):
911   AliFlowEventSimple(20), fApplyRecentering(kFALSE), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0)
912 {
913     // constructor
914     for(Int_t i(0); i < 9; i++) {
915         for(Int_t j(0); j < 2; j++) {
916             for(Int_t k(0); k < 2; k++) {
917                 fMeanQ[i][j][k] = 0.; 
918                 fWidthQ[i][j][k] = 0.;  
919                 fMeanQv3[i][j][k] = 0.; 
920                 fWidthQv3[i][j][k] = 0.;
921             }
922         }
923     }
924   //Fills the event from a vevent: AliESDEvent,AliAODEvent,AliMCEvent
925   //the input data needs to be attached to the cuts
926   Fill(rpCuts,poiCuts);
927 }
928
929 //-------------------------------------------------------------------//
930 //---- Including PMD tracks as RP --------------------------//
931
932 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
933                             const AliESDPmdTrack *pmdtracks,
934                             const AliCFManager* poiCFManager ):
935   AliFlowEventSimple(20), fApplyRecentering(kFALSE), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0)
936 {
937     // constructor
938     for(Int_t i(0); i < 9; i++) {
939         for(Int_t j(0); j < 2; j++) {
940             for(Int_t k(0); k < 2; k++) {
941                 fMeanQ[i][j][k] = 0.; 
942                 fWidthQ[i][j][k] = 0.;  
943                 fMeanQv3[i][j][k] = 0.; 
944                 fWidthQv3[i][j][k] = 0.;
945             }
946         }
947     }
948   Float_t GetPmdEta(Float_t xPos, Float_t yPos, Float_t zPos);
949   Float_t GetPmdPhi(Float_t xPos, Float_t yPos);
950   //Select the particles of interest from the ESD
951   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
952   
953   //loop over tracks
954   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
955     {
956       AliESDtrack* pParticle = anInput->GetTrack(itrkN);   //get input particle
957       //check if pParticle passes the cuts
958       Bool_t poiOK = kTRUE;
959       if (poiCFManager)
960         {
961           poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
962                     poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
963         }
964       if (!poiOK) continue;
965       
966       //make new AliFLowTrack
967       AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
968       
969       //marking the particles used for the particle of interest (POI) selection:
970       if(poiOK && poiCFManager)
971         {
972           IncrementNumberOfPOIs(1);
973           pTrack->SetForPOISelection(kTRUE);
974           pTrack->SetSource(AliFlowTrack::kFromESD);
975         }
976       
977       AddTrack(pTrack);
978     }//end of while (itrkN < iNumberOfInputTracks)
979   
980   //Select the reference particles from the PMD tracks
981   Int_t npmdcl = anInput->GetNumberOfPmdTracks();
982   printf("======There are %d PMD tracks in this event\n-------",npmdcl);
983   //loop over clusters 
984   for(Int_t iclust=0; iclust < npmdcl; iclust++){
985     //AliESDPmdTrack *pmdtr = anInput->GetPmdTrack(iclust);
986     pmdtracks = anInput->GetPmdTrack(iclust);
987     Int_t   det   = pmdtracks->GetDetector();
988     //Int_t   smn   = pmdtracks->GetSmn();
989     Float_t clsX  = pmdtracks->GetClusterX();
990     Float_t clsY  = pmdtracks->GetClusterY();
991     Float_t clsZ  = pmdtracks->GetClusterZ();
992     Float_t ncell = pmdtracks->GetClusterCells();
993     Float_t adc   = pmdtracks->GetClusterADC();
994     //Float_t pid   = pmdtracks->GetClusterPID();
995     Float_t etacls = GetPmdEta(clsX,clsY,clsZ);
996     Float_t phicls = GetPmdPhi(clsX,clsY);
997     //make new AliFLowTrackSimple
998     AliFlowTrack* pTrack = new AliFlowTrack();
999     //if(det == 0){ //selecting preshower plane only
1000     if(det == 0 && adc > 270 && ncell > 1){ //selecting preshower plane only
1001       //pTrack->SetPt(adc);//cluster adc
1002       pTrack->SetPt(0.0);
1003       pTrack->SetEta(etacls);
1004       pTrack->SetPhi(phicls);
1005       //marking the particles used for the reference particle (RP) selection:
1006       IncrementNumberOfPOIs(0);
1007       pTrack->SetForRPSelection(kTRUE);
1008       pTrack->SetSource(AliFlowTrack::kFromPMD);
1009       //Add the track to the flowevent
1010       AddTrack(pTrack);
1011     }//if det
1012   }
1013 }
1014 //----------------------------------------------------------------------------//
1015 Float_t GetPmdEta(Float_t xPos, Float_t yPos, Float_t zPos)
1016 {
1017   Float_t rpxpy, theta, eta;
1018   rpxpy  = TMath::Sqrt(xPos*xPos + yPos*yPos);
1019   theta  = TMath::ATan2(rpxpy,zPos);
1020   eta    = -TMath::Log(TMath::Tan(0.5*theta));
1021   return eta;
1022 }
1023 //--------------------------------------------------------------------------//
1024 Float_t GetPmdPhi(Float_t xPos, Float_t yPos)
1025 {
1026   Float_t pybypx, phi = 0., phi1;
1027   if(xPos==0)
1028     {
1029       if(yPos>0) phi = 90.;
1030       if(yPos<0) phi = 270.;
1031     }
1032   if(xPos != 0)
1033     {
1034       pybypx = yPos/xPos;
1035       if(pybypx < 0) pybypx = - pybypx;
1036       phi1 = TMath::ATan(pybypx)*180./3.14159;
1037       
1038       if(xPos > 0 && yPos > 0) phi = phi1;        // 1st Quadrant
1039       if(xPos < 0 && yPos > 0) phi = 180 - phi1;  // 2nd Quadrant
1040       if(xPos < 0 && yPos < 0) phi = 180 + phi1;  // 3rd Quadrant
1041       if(xPos > 0 && yPos < 0) phi = 360 - phi1;  // 4th Quadrant
1042       
1043     }
1044   phi = phi*3.14159/180.;
1045   return   phi;
1046 }
1047 //---------------------------------------------------------------//
1048
1049 void AliFlowEvent::Get2Qsub(AliFlowVector* Qarray, Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights)
1050 {
1051   // get q vectors for the subevents. if no recentering is necessary, get the guy from the flow event simple
1052   AliFlowEventSimple::Get2Qsub(Qarray, n, weightsList, usePhiWeights, usePtWeights, useEtaWeights);
1053   AliFlowVector vA = Qarray[0];
1054   AliFlowVector vB = Qarray[1];
1055  
1056   // else get the recentering from the cached info
1057   if (fApplyRecentering == 2010)        // 10h style recentering, implemented for n=2 and n=3
1058   {     
1059     // first retrieve the q-vectors from the AliFlowEventSimple:: routine
1060     // extract the information form the current flow vectors
1061     Double_t Qxc(vA.X());       // IMPORTANT: user is responsible for the sign of eta
1062     Double_t Qyc(vA.Y());       // vzeroC has negative pseudorapidity and is taken as subevent A
1063     Double_t Qxa(vB.X());       // vzeroA has positive pseudorapidity and is taken as subevent B
1064     Double_t Qya(vB.Y());
1065     // init some values for the corrections
1066     
1067     // default values for vector a (VZEROA)
1068     Double_t Qxamean(0);
1069     Double_t Qxarms(1);
1070     Double_t Qyamean(0);
1071     Double_t Qyarms(1);
1072     // default values for vector b (VZEROC)
1073     Double_t Qxcmean(0);
1074     Double_t Qxcrms(1);
1075     Double_t Qycmean(0);
1076     Double_t Qycrms(1); 
1077     // note that defaults are chosen such that for n!=2||n!=3 re-centering is a null-operation
1078     
1079     if( n == 2) {       // second order symmetry
1080         Qxamean = fMeanQ[fVZEROcentralityBin][1][0];
1081         Qxarms  = fWidthQ[fVZEROcentralityBin][1][0];
1082         Qyamean = fMeanQ[fVZEROcentralityBin][1][1];
1083         Qyarms  = fWidthQ[fVZEROcentralityBin][1][1];
1084
1085         Qxcmean = fMeanQ[fVZEROcentralityBin][0][0];
1086         Qxcrms  = fWidthQ[fVZEROcentralityBin][0][0];
1087         Qycmean = fMeanQ[fVZEROcentralityBin][0][1];
1088         Qycrms  = fWidthQ[fVZEROcentralityBin][0][1];   
1089     } else if (n == 3) {        // third order symmetry
1090         Qxamean = fMeanQv3[fVZEROcentralityBin][1][0];
1091         Qxarms  = fWidthQv3[fVZEROcentralityBin][1][0];
1092         Qyamean = fMeanQv3[fVZEROcentralityBin][1][1];
1093         Qyarms  = fWidthQv3[fVZEROcentralityBin][1][1];
1094   
1095         Qxcmean = fMeanQv3[fVZEROcentralityBin][0][0];
1096         Qxcrms  = fWidthQv3[fVZEROcentralityBin][0][0];
1097         Qycmean = fMeanQv3[fVZEROcentralityBin][0][1];
1098         Qycrms  = fWidthQv3[fVZEROcentralityBin][0][1]; 
1099     }
1100 //    printf(" > n %i Qx %.2f Qxrms %.2f \n", n, Qxamean, Qxarms);
1101     // do the correction    
1102     Double_t QxaCor = (Qxa - Qxamean)/Qxarms;
1103     Double_t QyaCor = (Qya - Qyamean)/Qyarms;
1104     Double_t QxcCor = (Qxc - Qxcmean)/Qxcrms;
1105     Double_t QycCor = (Qyc - Qycmean)/Qycrms;
1106     // update the vector
1107     vA.Set(QxcCor, QycCor);
1108     vB.Set(QxaCor, QyaCor);
1109   } else if (fApplyRecentering == 2011) { // 11h style recentering
1110
1111     Double_t QxaCor = 0.;
1112     Double_t QyaCor = 0.;
1113     Double_t QxcCor = 0.;
1114     Double_t QycCor = 0.;
1115
1116     if(fEvent && fEvent->GetEventplane()) {
1117       fEvent->GetEventplane()->CalculateVZEROEventPlane(fEvent, 8, n, QxaCor, QyaCor);
1118       fEvent->GetEventplane()->CalculateVZEROEventPlane(fEvent, 9, n, QxcCor, QycCor);
1119  
1120       // set the new q-vectors (which in this case means REPLACING) 
1121       vA.Set(QxcCor, QycCor);
1122       vB.Set(QxaCor, QyaCor);
1123     } // if for some reason the info from the event header is not available, the AliFlowTrackCuts object
1124       // has provided the equalized multiplicity info so this should still be relatively safe
1125   }
1126 }
1127 //_____________________________________________________________________________
1128 void AliFlowEvent::SetVZEROCalibrationForTrackCuts(AliFlowTrackCuts* cuts) {
1129     // open calibration info, copied from AliAnalyisTaskVnV0.cxx
1130     fEvent = cuts->GetEvent();
1131     if(!fEvent) return; // coverity. we need to know the event to get the runnumber and centrlaity
1132     // get the vzero centrality percentile (cc dependent calibration)
1133     Float_t v0Centr(fEvent->GetCentrality()->GetCentralityPercentile("V0M"));
1134     if(v0Centr < 5) fVZEROcentralityBin = 0;
1135     else if(v0Centr < 10) fVZEROcentralityBin = 1;
1136     else if(v0Centr < 20) fVZEROcentralityBin = 2;
1137     else if(v0Centr < 30) fVZEROcentralityBin = 3;
1138     else if(v0Centr < 40) fVZEROcentralityBin = 4;
1139     else if(v0Centr < 50) fVZEROcentralityBin = 5;
1140     else if(v0Centr < 60) fVZEROcentralityBin = 6;
1141     else if(v0Centr < 70) fVZEROcentralityBin = 7;
1142     else fVZEROcentralityBin = 8;
1143
1144     // if this event is from the same run as the previous event
1145     // we can use the cached calibration values, no need to re-open the 
1146     // aodb file, else cache the new run
1147     Int_t run(fEvent->GetRunNumber());
1148     if(fCachedRun == run) return;
1149     else fCachedRun = run;
1150     
1151     TFile *foadb = TFile::Open("$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root");
1152     if(!foadb){
1153         printf("OADB file $ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root cannot be opened, CALIBRATION FAILED !");
1154         return;
1155     }
1156
1157     AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
1158     if(!cont){
1159         printf("OADB object hMultV0BefCorr is not available in the file\n");
1160         return; 
1161     }
1162     if(!(cont->GetObject(run))){
1163         // if the multiplicity correction cannot be found for the specified run, 
1164         // loop over the 11h runs to see if it's 11h data
1165         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};
1166         for(Int_t r(0); r < 176; r++) {
1167             if(run == runs11h[r]) {
1168                 printf(" > run has been identified as 11h < \n");
1169  if(cuts->GetVZEROgainEqualizationPerRing()) {
1170                     // enable or disable rings through the weights, weight 1. is enabled, 0. is disabled
1171                     // start with the vzero c rings (segments 0 through 31)
1172                     (cuts->GetUseVZERORing(0)) ? cuts->SetVZEROCpol(0, 1.) : cuts->SetVZEROCpol(0, 0.);
1173                     (cuts->GetUseVZERORing(1)) ? cuts->SetVZEROCpol(1, 1.) : cuts->SetVZEROCpol(1, 0.);
1174                     (cuts->GetUseVZERORing(2)) ? cuts->SetVZEROCpol(2, 1.) : cuts->SetVZEROCpol(2, 0.);
1175                     (cuts->GetUseVZERORing(3)) ? cuts->SetVZEROCpol(3, 1.) : cuts->SetVZEROCpol(3, 0.);
1176                     // same for vzero a
1177                     (cuts->GetUseVZERORing(4)) ? cuts->SetVZEROApol(0, 1.) : cuts->SetVZEROApol(0, 0.);
1178                     (cuts->GetUseVZERORing(5)) ? cuts->SetVZEROApol(1, 1.) : cuts->SetVZEROApol(1, 0.);
1179                     (cuts->GetUseVZERORing(6)) ? cuts->SetVZEROApol(2, 1.) : cuts->SetVZEROApol(2, 0.);
1180                     (cuts->GetUseVZERORing(7)) ? cuts->SetVZEROApol(3, 1.) : cuts->SetVZEROApol(3, 0.);
1181                 } else {
1182                     // else enable all rings, which is also default
1183                     for(Int_t i(0); i < 4; i++) cuts->SetVZEROCpol(i, 1.);
1184                     for(Int_t i(0); i < 4; i++) cuts->SetVZEROApol(i, 1.);
1185                 }
1186                 // pass a NULL pointer to the track cuts object, the NULL pointer will identify 11h runs
1187                 cuts->SetVZEROgainEqualisation(NULL);
1188                 fApplyRecentering = 2011;
1189                 return; // the rest of the steps are not necessary
1190             }
1191         }
1192         // the run has not been identified as lhc11h data, so we assume a template calibration
1193         printf("OADB object hMultVZEROBefCorr is not available for run %i (used default run 137366)\n",run);
1194         run = 137366;
1195     }
1196     printf(" > run has been identified as 10h < \n");
1197     // step 1) get the proper multiplicity weights from the vzero signal
1198     TProfile* fMultVZERO = ((TH2F *) cont->GetObject(run))->ProfileX();
1199
1200     TF1 *fpol0 = new TF1("fpol0","pol0"); 
1201     if(cuts->GetVZEROgainEqualizationPerRing()) {
1202         // do the calibration per ring
1203         // start with the vzero c rings (segments 0 through 31)
1204         fMultVZERO->Fit(fpol0, "", "", 0, 8);
1205         (cuts->GetUseVZERORing(0)) ? cuts->SetVZEROCpol(0, fpol0->GetParameter(0)) : cuts->SetVZEROCpol(0, 0.);
1206         fMultVZERO->Fit(fpol0, "", "", 8, 16);
1207         (cuts->GetUseVZERORing(1)) ? cuts->SetVZEROCpol(1, fpol0->GetParameter(0)) : cuts->SetVZEROCpol(1, 0.);
1208         fMultVZERO->Fit(fpol0, "", "", 16, 24);
1209         (cuts->GetUseVZERORing(2)) ? cuts->SetVZEROCpol(2, fpol0->GetParameter(0)) : cuts->SetVZEROCpol(2, 0.);
1210         fMultVZERO->Fit(fpol0, "", "", 24, 32);
1211         (cuts->GetUseVZERORing(3)) ? cuts->SetVZEROCpol(3, fpol0->GetParameter(0)) : cuts->SetVZEROCpol(3, 0.);
1212         // same thing for vero A
1213         fMultVZERO->Fit(fpol0, "", "", 32, 40);
1214         (cuts->GetUseVZERORing(4)) ? cuts->SetVZEROApol(0, fpol0->GetParameter(0)) : cuts->SetVZEROApol(0, 0.);
1215         fMultVZERO->Fit(fpol0, "", "", 40, 48);
1216         (cuts->GetUseVZERORing(5)) ? cuts->SetVZEROApol(1, fpol0->GetParameter(0)) : cuts->SetVZEROApol(1, 0.);
1217         fMultVZERO->Fit(fpol0, "", "", 48, 56);
1218         (cuts->GetUseVZERORing(6)) ? cuts->SetVZEROApol(2, fpol0->GetParameter(0)) : cuts->SetVZEROApol(2, 0.);
1219         fMultVZERO->Fit(fpol0, "", "", 56, 64);
1220         (cuts->GetUseVZERORing(7)) ? cuts->SetVZEROApol(3, fpol0->GetParameter(0)) : cuts->SetVZEROApol(3, 0.);
1221     } else {
1222         // do the calibration in one go. the calibration will still be 
1223         // stored per ring, but each ring has the same weight now
1224        fMultVZERO->Fit(fpol0,"","",0,31);
1225        for(Int_t i(0); i < 4; i++) cuts->SetVZEROCpol(i, fpol0->GetParameter(0));
1226        fMultVZERO->Fit(fpol0,"","",32,64);
1227        for(Int_t i(0); i < 4; i++) cuts->SetVZEROApol(i, fpol0->GetParameter(0));
1228     }
1229     // the parameters to weigh the vzero track cuts have been extracted now, 
1230     // so we can pass them to the current track cuts obect
1231     cuts->SetVZEROgainEqualisation(fMultVZERO);       // passed as a TH1
1232
1233     // step 2) reweight the q-vectors that will be  called by flow methods which use
1234     // subevents
1235     // underlying assumption is that subevent a uses VZEROA
1236     // and subevent b uses VZEROC
1237     for(Int_t iside=0;iside<2;iside++){
1238         for(Int_t icoord=0;icoord<2;icoord++){
1239             for(Int_t i=0;i  < 9;i++){
1240                 char namecont[100];
1241                 if(iside==0 && icoord==0)
1242                   snprintf(namecont,100,"hQxc2_%i",i);
1243                 else if(iside==1 && icoord==0)
1244                   snprintf(namecont,100,"hQxa2_%i",i);
1245                 else if(iside==0 && icoord==1)
1246                   snprintf(namecont,100,"hQyc2_%i",i);
1247                 else if(iside==1 && icoord==1)
1248                   snprintf(namecont,100,"hQya2_%i",i);
1249
1250                 cont = (AliOADBContainer*) foadb->Get(namecont);
1251                 if(!cont){
1252                     printf("OADB object %s is not available in the file\n",namecont);
1253                     return;     
1254                 }
1255         
1256                 if(!(cont->GetObject(run))){
1257                     printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
1258                     run = 137366;
1259                 }
1260
1261                 // after grabbing all the info, set the CORRECTION TERMS to
1262                 // the 2nd and 3rd order qsub-vectors
1263                 // we do this here for all centralities, so that subsequent events
1264                 // can grab the correction from these cached values
1265                 fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
1266                 fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
1267
1268                 //for v3
1269                 if(iside==0 && icoord==0)
1270                   snprintf(namecont,100,"hQxc3_%i",i);
1271                 else if(iside==1 && icoord==0)
1272                   snprintf(namecont,100,"hQxa3_%i",i);
1273                 else if(iside==0 && icoord==1)
1274                   snprintf(namecont,100,"hQyc3_%i",i);
1275                 else if(iside==1 && icoord==1)
1276                   snprintf(namecont,100,"hQya3_%i",i);
1277
1278                 cont = (AliOADBContainer*) foadb->Get(namecont);
1279                 if(!cont){
1280                     printf("OADB object %s is not available in the file\n",namecont);
1281                     return;     
1282                 }
1283                 
1284                 if(!(cont->GetObject(run))){
1285                     printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
1286                     run = 137366;
1287                 }
1288                 fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
1289                 fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
1290
1291             }
1292         }
1293     }
1294     // set the recentering style (might be switched back to -1 if recentering is disabeled)
1295     fApplyRecentering = 2010;
1296 }
1297 //_____________________________________________________________________________
1298
1299 void AliFlowEvent::ClearFast()
1300 {
1301   //clear the event without releasing any memory
1302   //note that cached run number of recentering settigns are not clear 
1303   //(see AliFlowEvent::ClearCachedRun() )
1304   AliFlowEventSimple::ClearFast();
1305 }
1306
1307 //_____________________________________________________________________________
1308
1309 void AliFlowEvent::ClearCachedRun()
1310 {
1311     //clear the cached run (not in clear fast as cache needs to be persistent in most cases )
1312   fCachedRun=0;  
1313   fApplyRecentering=0;
1314 }