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