updates from redmer
[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), fCurrentCentrality(-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), fCurrentCentrality(-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), fCurrentCentrality(-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   fCurrentCentrality = event.fCurrentCentrality;
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), fCurrentCentrality(-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), fCurrentCentrality(-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), fCurrentCentrality(-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), fCurrentCentrality(-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), fCurrentCentrality(-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), fCurrentCentrality(-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         if(hybrid)
632           tpcTrack->PropagateToDCA(vertexSPD,esd->GetMagneticField(),100.,dca,cov);
633         else
634           tpcTrack->PropagateToDCA(vertexTPC,esd->GetMagneticField(),100.,dca,cov);
635
636 //        dca3D = TMath::Sqrt(TMath::Power(dca[0],2)+TMath::Power(dca[1],2));   FIXME unused variable
637
638       }
639
640       //make new AliFLowTrack
641       AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
642
643       pTrack->SetSource(AliFlowTrack::kFromESD);
644
645       //marking the particles used for diff. flow:
646       if(poiOK && poiCFManager)
647       {
648         pTrack->SetForPOISelection(kTRUE);
649         IncrementNumberOfPOIs(1);
650       }
651
652       if(useTPC)
653       {
654         pTrack->SetForRPSelection(kTRUE);
655         IncrementNumberOfPOIs(0);
656       }
657
658       AddTrack(pTrack);
659
660     }//end of while (itrkN < iNumberOfInputTracks)
661
662 }
663
664 //-----------------------------------------------------------------------
665 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
666                             const TH2F* anInputFMDhist,
667                             const AliCFManager* poiCFManager ):
668   AliFlowEventSimple(20), fApplyRecentering(-1), fCachedRun(-1), fCurrentCentrality(-1)
669 {
670     // constructor
671     for(Int_t i(0); i < 9; i++) {
672         for(Int_t j(0); j < 2; j++) {
673             for(Int_t k(0); k < 2; k++) {
674                 fMeanQ[i][j][k] = 0.; 
675                 fWidthQ[i][j][k] = 0.;  
676                 fMeanQv3[i][j][k] = 0.; 
677                 fWidthQv3[i][j][k] = 0.;
678             }
679         }
680     }
681
682   //Select the particles of interest from the ESD
683   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
684
685   //loop over tracks
686   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
687     {
688       AliESDtrack* pParticle = anInput->GetTrack(itrkN);   //get input particle
689
690       //check if pParticle passes the cuts
691       Bool_t poiOK = kTRUE;
692       if (poiCFManager)
693         {
694           poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
695                     poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
696         }
697       if (!poiOK) continue;
698  
699       //make new AliFLowTrack
700       AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
701           
702       //marking the particles used for the particle of interest (POI) selection:
703       if(poiOK && poiCFManager)
704         {
705           IncrementNumberOfPOIs(1);
706           pTrack->SetForPOISelection(kTRUE);
707           pTrack->SetSource(AliFlowTrack::kFromESD);
708         }
709
710       AddTrack(pTrack);
711     }//end of while (itrkN < iNumberOfInputTracks)
712
713   //Select the reference particles from the FMD hits
714   //loop over FMD histogram
715   Int_t iBinsEta = anInputFMDhist->GetNbinsX();
716   Int_t iBinsPhi = anInputFMDhist->GetNbinsY();
717   
718   for (Int_t iEta = 1; iEta <= iBinsEta; iEta++){
719     Double_t etaFMD = anInputFMDhist->GetXaxis()->GetBinCenter(iEta);
720     for (Int_t iPhi = 1; iPhi <= iBinsPhi; iPhi++){
721       Double_t phiFMD = anInputFMDhist->GetYaxis()->GetBinCenter(iPhi);
722       Double_t weightFMD = anInputFMDhist->GetBinContent(iEta,iPhi);
723     
724       if (weightFMD > 0.0) { //do not add empty bins
725         //make new AliFLowTrackSimple
726         AliFlowTrack* pTrack = new AliFlowTrack();
727         pTrack->SetPt(0.0);
728         pTrack->SetEta(etaFMD);
729         pTrack->SetPhi(phiFMD);
730         pTrack->SetWeight(weightFMD);
731         //marking the particles used for the reference particle (RP) selection:
732         pTrack->TagRP();
733         IncrementNumberOfPOIs(0);
734         pTrack->SetSource(AliFlowTrack::kFromFMD);
735
736         //Add the track to the flowevent
737         AddTrack(pTrack);
738         
739       }
740     }
741   }
742
743 }
744
745 //-----------------------------------------------------------------------
746 void AliFlowEvent::FindDaughters(Bool_t keepDaughtersInRPselection)
747 {
748   //each flow track holds it's esd track index as well as its daughters esd index.
749   //fill the array of daughters for every track with the pointers to flow tracks
750   //to associate the mothers with daughters directly
751   for (Int_t iTrack=0; iTrack<fMothersCollection->GetEntriesFast(); iTrack++)
752   {
753     AliFlowTrack* mother = static_cast<AliFlowTrack*>(fMothersCollection->At(iTrack));
754     if (!mother) continue;
755     if (mother->GetNDaughters()<1) continue;
756     for (Int_t iDaughterCandidate=0; iDaughterCandidate<fNumberOfTracks; iDaughterCandidate++)
757     {
758       AliFlowTrack* daughterCandidate = static_cast<AliFlowTrack*>(fTrackCollection->At(iDaughterCandidate));
759       Int_t esdIndexDaughterCandidate = daughterCandidate->GetID();
760       for (Int_t iDaughter=0; iDaughter<mother->GetNDaughters(); iDaughter++)
761       {
762         Int_t esdIndexDaughter = mother->GetIDDaughter(iDaughter);
763         if (esdIndexDaughter==esdIndexDaughterCandidate)
764         {
765           mother->SetDaughter(iDaughter,daughterCandidate);
766           daughterCandidate->SetForRPSelection(keepDaughtersInRPselection);
767         }
768       }
769     }
770   }
771 }
772
773 //-----------------------------------------------------------------------
774 void AliFlowEvent::Fill( AliFlowTrackCuts* rpCuts,
775                          AliFlowTrackCuts* poiCuts )
776 {
777   //Fills the event from a vevent: AliESDEvent,AliAODEvent,AliMCEvent
778   //the input data needs to be attached to the cuts
779   //we have two cases, if we're cutting the same collection of tracks
780   //(same param type) then we can have tracks that are both rp and poi
781   //in the other case we want to have two exclusive sets of rps and pois
782   //e.g. one tracklets, the other PMD or global - USER IS RESPOSIBLE
783   //FOR MAKING SURE THEY DONT OVERLAP OR ELSE THE SAME PARTICLE WILL BE
784   //TAKEN TWICE
785
786   ClearFast();
787
788   if (!rpCuts || !poiCuts) return;
789   AliFlowTrackCuts::trackParameterType sourceRP = rpCuts->GetParamType();
790   AliFlowTrackCuts::trackParameterType sourcePOI = poiCuts->GetParamType();
791   AliFlowTrack* pTrack=NULL;
792
793   if (sourceRP==sourcePOI)
794   {
795     //loop over tracks
796     Int_t numberOfInputObjects = rpCuts->GetNumberOfInputObjects();
797     for (Int_t i=0; i<numberOfInputObjects; i++)
798     {
799       //get input object (particle)
800       TObject* particle = rpCuts->GetInputObject(i);
801
802       Bool_t rp = rpCuts->IsSelected(particle,i);
803       Int_t poiClass = poiCuts->IsSelected(particle,i);
804
805       if (!(rp||poiClass>0)) continue;
806
807       //make new AliFlowTrack
808       if (rp)
809       {
810         pTrack = rpCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
811         if (!pTrack) continue;
812         pTrack->TagRP(); IncrementNumberOfPOIs(0);
813         if (poiClass>0) {pTrack->Tag(poiClass); IncrementNumberOfPOIs(poiClass);}
814         if (pTrack->GetNDaughters()>0) fMothersCollection->Add(pTrack);
815       }
816       else if (poiClass>0)
817       {
818         pTrack = poiCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
819         if (!pTrack) continue;
820         pTrack->Tag(poiClass); IncrementNumberOfPOIs(poiClass);
821         if (pTrack->GetNDaughters()>0) fMothersCollection->Add(pTrack);
822       }
823       fNumberOfTracks++;
824     }//end of while (i < numberOfTracks)
825   }
826   else if (sourceRP!=sourcePOI)
827   {
828     //here we have two different sources of particles, so we fill
829     //them independently
830     //POI
831     for (Int_t i=0; i<poiCuts->GetNumberOfInputObjects(); i++)
832     {
833       TObject* particle = poiCuts->GetInputObject(i);
834       Int_t poiClass = poiCuts->IsSelected(particle,i);
835       if (poiClass<=0) continue;
836       pTrack = poiCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
837       if (!pTrack) continue;
838       pTrack->Tag(poiClass);
839       IncrementNumberOfPOIs(poiClass);
840       fNumberOfTracks++;
841       if (pTrack->GetNDaughters()>0) fMothersCollection->Add(pTrack);
842     }
843     //RP
844     Int_t numberOfInputObjects = rpCuts->GetNumberOfInputObjects();
845     for (Int_t i=0; i<numberOfInputObjects; i++)
846       {
847       TObject* particle = rpCuts->GetInputObject(i);
848       Int_t rp = rpCuts->IsSelected(particle,i);
849       if (rp<1) continue;
850       pTrack = rpCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
851       if (!pTrack) continue;
852       pTrack->TagRP();
853       IncrementNumberOfPOIs(0);
854       fNumberOfTracks++;
855       if (pTrack->GetNDaughters()>0) fMothersCollection->Add(pTrack);
856     }
857   }
858 }
859
860 //-----------------------------------------------------------------------
861 void AliFlowEvent::InsertTrack(AliFlowTrack *track) {
862   // adds a flow track at the end of the container
863   AliFlowTrack *pTrack = ReuseTrack( fNumberOfTracks++ );
864   *pTrack = *track;
865   if (track->GetNDaughters()>0)
866   {
867     fMothersCollection->Add(track);
868   }
869   //pTrack->SetPt( track->Pt() );
870   //pTrack->SetPhi( track->Phi() );
871   //pTrack->SetEta( track->Eta() );
872   //pTrack->SetWeight( track->Weight() );
873   //pTrack->SetCharge( track->Charge() );
874   //pTrack->SetMass( track->Mass() );
875   //pTrack->SetForRPSelection( track->InRPSelection() );
876   //pTrack->SetForPOISelection( track->InPOISelection() );
877   //if(track->InSubevent(0)) pTrack->SetForSubevent(0);
878   //if(track->InSubevent(1)) pTrack->SetForSubevent(1);
879   //pTrack->SetID( track->GetID() );
880   return;
881 }
882
883 //-----------------------------------------------------------------------
884 AliFlowTrack* AliFlowEvent::ReuseTrack(Int_t i)
885 {
886   //try to reuse an existing track, if empty, make new one
887   AliFlowTrack* pTrack = static_cast<AliFlowTrack*>(fTrackCollection->At(i));
888   if (pTrack)
889   {
890     pTrack->Clear();
891   }
892   else 
893   {
894     pTrack = new AliFlowTrack();
895     fTrackCollection->AddAtAndExpand(pTrack,i);
896   }
897   return pTrack;
898 }
899
900 //-----------------------------------------------------------------------
901 AliFlowEvent::AliFlowEvent( AliFlowTrackCuts* rpCuts,
902                             AliFlowTrackCuts* poiCuts ):
903   AliFlowEventSimple(20), fApplyRecentering(kFALSE), fCachedRun(-1), fCurrentCentrality(-1)
904 {
905     // constructor
906     for(Int_t i(0); i < 9; i++) {
907         for(Int_t j(0); j < 2; j++) {
908             for(Int_t k(0); k < 2; k++) {
909                 fMeanQ[i][j][k] = 0.; 
910                 fWidthQ[i][j][k] = 0.;  
911                 fMeanQv3[i][j][k] = 0.; 
912                 fWidthQv3[i][j][k] = 0.;
913             }
914         }
915     }
916   //Fills the event from a vevent: AliESDEvent,AliAODEvent,AliMCEvent
917   //the input data needs to be attached to the cuts
918   //we have two cases, if we're cutting the same collection of tracks
919   //(same param type) then we can have tracks that are both rp and poi
920   //in the other case we want to have two exclusive sets of rps and pois
921   //e.g. one tracklets, the other PMD or global - USER IS RESPOSIBLE
922   //FOR MAKING SURE THEY DONT OVERLAP OR ELSE THE SAME PARTICLE WILL BE
923   //TAKEN TWICE
924
925   if (!rpCuts || !poiCuts) return;
926   AliFlowTrackCuts::trackParameterType sourceRP = rpCuts->GetParamType();
927   AliFlowTrackCuts::trackParameterType sourcePOI = poiCuts->GetParamType();
928
929   if (sourceRP==sourcePOI)
930   {
931     //loop over tracks
932     Int_t numberOfInputObjects = rpCuts->GetNumberOfInputObjects();
933     for (Int_t i=0; i<numberOfInputObjects; i++)
934     {
935       //get input object (particle)
936       TObject* particle = rpCuts->GetInputObject(i);
937
938       Bool_t rp = rpCuts->IsSelected(particle,i);
939       Int_t poiClass = poiCuts->IsSelected(particle,i);
940       
941       if (!(rp||poiClass>0)) continue;
942
943       //make new AliFLowTrack
944       AliFlowTrack* pTrack = NULL;
945       if (rp)
946       {
947         pTrack = rpCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
948         if (!pTrack) continue;
949         pTrack->TagRP(); IncrementNumberOfPOIs(0);
950         if (poiClass>0) {pTrack->Tag(poiClass); IncrementNumberOfPOIs(poiClass);}
951       }
952       else
953       if (poiClass>0)
954       {
955         pTrack = poiCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
956         if (!pTrack) continue;
957         pTrack->Tag(poiClass); IncrementNumberOfPOIs(poiClass);
958       }
959       TrackAdded();
960     }//end of while (i < numberOfTracks)
961   }
962   else if (sourceRP!=sourcePOI)
963   {
964     //here we have two different sources of particles, so we fill
965     //them independently
966     AliFlowTrack* pTrack = NULL;
967     //RP
968     Int_t numberOfInputObjects = rpCuts->GetNumberOfInputObjects();
969     for (Int_t i=0; i<numberOfInputObjects; i++)
970     {
971       TObject* particle = rpCuts->GetInputObject(i);
972       Bool_t rp = rpCuts->IsSelected(particle,i);
973       if (!rp) continue;
974       pTrack = rpCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
975       if (!pTrack) continue;
976       pTrack->TagRP(); IncrementNumberOfPOIs(0);
977       TrackAdded();
978     }
979     //POI
980     numberOfInputObjects = poiCuts->GetNumberOfInputObjects();
981     for (Int_t i=0; i<numberOfInputObjects; i++)
982     {
983       TObject* particle = poiCuts->GetInputObject(i);
984       Int_t poiClass = poiCuts->IsSelected(particle,i);
985       if (poiClass<=0) continue;
986       pTrack = poiCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
987       if (!pTrack) continue;
988       pTrack->Tag(poiClass); IncrementNumberOfPOIs(poiClass);
989       TrackAdded();
990     }
991   }
992 }
993
994 //-------------------------------------------------------------------//
995 //---- Including PMD tracks as RP --------------------------//
996
997 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
998                             const AliESDPmdTrack *pmdtracks,
999                             const AliCFManager* poiCFManager ):
1000   AliFlowEventSimple(20), fApplyRecentering(kFALSE), fCachedRun(-1), fCurrentCentrality(-1)
1001 {
1002     // constructor
1003     for(Int_t i(0); i < 9; i++) {
1004         for(Int_t j(0); j < 2; j++) {
1005             for(Int_t k(0); k < 2; k++) {
1006                 fMeanQ[i][j][k] = 0.; 
1007                 fWidthQ[i][j][k] = 0.;  
1008                 fMeanQv3[i][j][k] = 0.; 
1009                 fWidthQv3[i][j][k] = 0.;
1010             }
1011         }
1012     }
1013   Float_t GetPmdEta(Float_t xPos, Float_t yPos, Float_t zPos);
1014   Float_t GetPmdPhi(Float_t xPos, Float_t yPos);
1015   //Select the particles of interest from the ESD
1016   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
1017   
1018   //loop over tracks
1019   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
1020     {
1021       AliESDtrack* pParticle = anInput->GetTrack(itrkN);   //get input particle
1022       //check if pParticle passes the cuts
1023       Bool_t poiOK = kTRUE;
1024       if (poiCFManager)
1025         {
1026           poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
1027                     poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
1028         }
1029       if (!poiOK) continue;
1030       
1031       //make new AliFLowTrack
1032       AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
1033       
1034       //marking the particles used for the particle of interest (POI) selection:
1035       if(poiOK && poiCFManager)
1036         {
1037           IncrementNumberOfPOIs(1);
1038           pTrack->SetForPOISelection(kTRUE);
1039           pTrack->SetSource(AliFlowTrack::kFromESD);
1040         }
1041       
1042       AddTrack(pTrack);
1043     }//end of while (itrkN < iNumberOfInputTracks)
1044   
1045   //Select the reference particles from the PMD tracks
1046   Int_t npmdcl = anInput->GetNumberOfPmdTracks();
1047   printf("======There are %d PMD tracks in this event\n-------",npmdcl);
1048   //loop over clusters 
1049   for(Int_t iclust=0; iclust < npmdcl; iclust++){
1050     //AliESDPmdTrack *pmdtr = anInput->GetPmdTrack(iclust);
1051     pmdtracks = anInput->GetPmdTrack(iclust);
1052     Int_t   det   = pmdtracks->GetDetector();
1053     //Int_t   smn   = pmdtracks->GetSmn();
1054     Float_t clsX  = pmdtracks->GetClusterX();
1055     Float_t clsY  = pmdtracks->GetClusterY();
1056     Float_t clsZ  = pmdtracks->GetClusterZ();
1057     Float_t ncell = pmdtracks->GetClusterCells();
1058     Float_t adc   = pmdtracks->GetClusterADC();
1059     //Float_t pid   = pmdtracks->GetClusterPID();
1060     Float_t etacls = GetPmdEta(clsX,clsY,clsZ);
1061     Float_t phicls = GetPmdPhi(clsX,clsY);
1062     //make new AliFLowTrackSimple
1063     AliFlowTrack* pTrack = new AliFlowTrack();
1064     //if(det == 0){ //selecting preshower plane only
1065     if(det == 0 && adc > 270 && ncell > 1){ //selecting preshower plane only
1066       //pTrack->SetPt(adc);//cluster adc
1067       pTrack->SetPt(0.0);
1068       pTrack->SetEta(etacls);
1069       pTrack->SetPhi(phicls);
1070       //marking the particles used for the reference particle (RP) selection:
1071       IncrementNumberOfPOIs(0);
1072       pTrack->SetForRPSelection(kTRUE);
1073       pTrack->SetSource(AliFlowTrack::kFromPMD);
1074       //Add the track to the flowevent
1075       AddTrack(pTrack);
1076     }//if det
1077   }
1078 }
1079 //----------------------------------------------------------------------------//
1080 Float_t GetPmdEta(Float_t xPos, Float_t yPos, Float_t zPos)
1081 {
1082   Float_t rpxpy, theta, eta;
1083   rpxpy  = TMath::Sqrt(xPos*xPos + yPos*yPos);
1084   theta  = TMath::ATan2(rpxpy,zPos);
1085   eta    = -TMath::Log(TMath::Tan(0.5*theta));
1086   return eta;
1087 }
1088 //--------------------------------------------------------------------------//
1089 Float_t GetPmdPhi(Float_t xPos, Float_t yPos)
1090 {
1091   Float_t pybypx, phi = 0., phi1;
1092   if(xPos==0)
1093     {
1094       if(yPos>0) phi = 90.;
1095       if(yPos<0) phi = 270.;
1096     }
1097   if(xPos != 0)
1098     {
1099       pybypx = yPos/xPos;
1100       if(pybypx < 0) pybypx = - pybypx;
1101       phi1 = TMath::ATan(pybypx)*180./3.14159;
1102       
1103       if(xPos > 0 && yPos > 0) phi = phi1;        // 1st Quadrant
1104       if(xPos < 0 && yPos > 0) phi = 180 - phi1;  // 2nd Quadrant
1105       if(xPos < 0 && yPos < 0) phi = 180 + phi1;  // 3rd Quadrant
1106       if(xPos > 0 && yPos < 0) phi = 360 - phi1;  // 4th Quadrant
1107       
1108     }
1109   phi = phi*3.14159/180.;
1110   return   phi;
1111 }
1112 //---------------------------------------------------------------//
1113
1114 //_____________________________________________________________________________
1115 void AliFlowEvent::SetVZEROCalibrationForTrackCuts(AliFlowTrackCuts* cuts) {
1116     // open calibration info, copied from AliAnalyisTaskVnV0.cxx
1117     if(!cuts->GetEvent()) return; // coverity. we need to know the event to get the runnumber and centrlaity
1118     // get the vzero centrality percentile (cc dependent calibration)
1119     Float_t v0Centr(cuts->GetEvent()->GetCentrality()->GetCentralityPercentile("V0M"));
1120     if(v0Centr < 5) fCurrentCentrality = 0;
1121     else if(v0Centr < 10) fCurrentCentrality = 1;
1122     else if(v0Centr < 20) fCurrentCentrality = 2;
1123     else if(v0Centr < 30) fCurrentCentrality = 3;
1124     else if(v0Centr < 40) fCurrentCentrality = 4;
1125     else if(v0Centr < 50) fCurrentCentrality = 5;
1126     else if(v0Centr < 60) fCurrentCentrality = 6;
1127     else if(v0Centr < 70) fCurrentCentrality = 7;
1128     else fCurrentCentrality = 8;
1129
1130     // if this event is from the same run as the previous event
1131     // we can use the cached calibration values, no need to re-open the 
1132     // aodb file
1133     Int_t run(cuts->GetEvent()->GetRunNumber());
1134 //    printf ( " > run number is %i \n", run);
1135     if(fCachedRun == run) {
1136         // the runnumber did not change, no need to open the database again
1137         // in case of 11h style recentering, update the q-sub vectors
1138         if(fApplyRecentering == 2011) SetVZEROCalibrationForTrackCuts2011(cuts); 
1139         return;
1140     }
1141     // set the chached run number
1142     fCachedRun = run;
1143     
1144     TString oadbfilename = "$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root";
1145     TFile *foadb = TFile::Open(oadbfilename.Data());
1146
1147     if(!foadb){
1148         printf("OADB file %s cannot be opened\n",oadbfilename.Data());
1149         return;
1150     }
1151
1152     AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
1153     if(!cont){
1154         printf("OADB object hMultV0BefCorr is not available in the file\n");
1155         return; 
1156     }
1157     if(!(cont->GetObject(run))){
1158         // if the multiplicity correction cannot be found for the specified run, 
1159         // loop over the 11h runs to see if it's 11h data
1160         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};
1161         for(Int_t r(0); r < 176; r++) {
1162             if(run == runs11h[r]) {
1163                 printf(" > run has been identified as 11h < \n");
1164                 if(cuts->GetV0gainEqualizationPerRing()) {
1165                     // enable or disable rings through the weights, weight 1. is enabled, 0. is disabled
1166                     // start with the vzero c rings (segments 0 through 31)
1167                     (cuts->GetUseVZERORing(0)) ? cuts->SetV0Cpol(0, 1.) : cuts->SetV0Cpol(0, 0.);
1168                     (cuts->GetUseVZERORing(1)) ? cuts->SetV0Cpol(1, 1.) : cuts->SetV0Cpol(1, 0.);
1169                     (cuts->GetUseVZERORing(2)) ? cuts->SetV0Cpol(2, 1.) : cuts->SetV0Cpol(2, 0.);
1170                     (cuts->GetUseVZERORing(3)) ? cuts->SetV0Cpol(3, 1.) : cuts->SetV0Cpol(3, 0.);
1171                     // same for vzero a
1172                     (cuts->GetUseVZERORing(4)) ? cuts->SetV0Apol(0, 1.) : cuts->SetV0Apol(0, 0.);
1173                     (cuts->GetUseVZERORing(5)) ? cuts->SetV0Apol(1, 1.) : cuts->SetV0Apol(1, 0.);
1174                     (cuts->GetUseVZERORing(6)) ? cuts->SetV0Apol(2, 1.) : cuts->SetV0Apol(2, 0.);
1175                     (cuts->GetUseVZERORing(7)) ? cuts->SetV0Apol(3, 1.) : cuts->SetV0Apol(3, 0.);
1176                 } else {
1177                     // else enable all rings
1178                     for(Int_t i(0); i < 4; i++) cuts->SetV0Cpol(i, 1.);
1179                     for(Int_t i(0); i < 4; i++) cuts->SetV0Apol(i, 1.);
1180                 }
1181                 // pass a NULL pointer to the track cuts object
1182                 // the NULL pointer will identify 11h runs
1183                 cuts->SetV0gainEqualisation(NULL);
1184                 // this will identify the recentering style that is required. flight might be changed if recenetering is disabled
1185                 fApplyRecentering = 2011;
1186                 SetVZEROCalibrationForTrackCuts2011(cuts); 
1187                 return; // the rest of the steps are not necessary
1188             }
1189         }
1190         // the run has not been identified as lhc11h data, so we assume a template calibration
1191         printf("OADB object hMultV0BefCorr is not available for run %i (used run 137366)\n",run);
1192         run = 137366;
1193     }
1194     printf(" > run has been identified as 10h < \n");
1195     // step 1) get the proper multiplicity weights from the vzero signal
1196     TProfile* fMultV0 = ((TH2F *) cont->GetObject(run))->ProfileX();
1197
1198     TF1 *fpol0 = new TF1("fpol0","pol0"); 
1199     if(cuts->GetV0gainEqualizationPerRing()) {
1200         // do the calibration per ring
1201         // start with the vzero c rings (segments 0 through 31)
1202         fMultV0->Fit(fpol0, "", "", 0, 8);
1203         (cuts->GetUseVZERORing(0)) ? cuts->SetV0Cpol(0, fpol0->GetParameter(0)) : cuts->SetV0Cpol(0, 0.);
1204         fMultV0->Fit(fpol0, "", "", 8, 16);
1205         (cuts->GetUseVZERORing(1)) ? cuts->SetV0Cpol(1, fpol0->GetParameter(0)) : cuts->SetV0Cpol(1, 0.);
1206         fMultV0->Fit(fpol0, "", "", 16, 24);
1207         (cuts->GetUseVZERORing(2)) ? cuts->SetV0Cpol(2, fpol0->GetParameter(0)) : cuts->SetV0Cpol(2, 0.);
1208         fMultV0->Fit(fpol0, "", "", 24, 32);
1209         (cuts->GetUseVZERORing(3)) ? cuts->SetV0Cpol(3, fpol0->GetParameter(0)) : cuts->SetV0Cpol(3, 0.);
1210         // same thing for vero A
1211         fMultV0->Fit(fpol0, "", "", 32, 40);
1212         (cuts->GetUseVZERORing(4)) ? cuts->SetV0Apol(0, fpol0->GetParameter(0)) : cuts->SetV0Apol(0, 0.);
1213         fMultV0->Fit(fpol0, "", "", 40, 48);
1214         (cuts->GetUseVZERORing(5)) ? cuts->SetV0Apol(1, fpol0->GetParameter(0)) : cuts->SetV0Apol(1, 0.);
1215         fMultV0->Fit(fpol0, "", "", 48, 56);
1216         (cuts->GetUseVZERORing(6)) ? cuts->SetV0Apol(2, fpol0->GetParameter(0)) : cuts->SetV0Apol(2, 0.);
1217         fMultV0->Fit(fpol0, "", "", 56, 64);
1218         (cuts->GetUseVZERORing(7)) ? cuts->SetV0Apol(3, fpol0->GetParameter(0)) : cuts->SetV0Apol(3, 0.);
1219     } else {
1220         // do the calibration in one go. the calibration will still be 
1221         // stored per ring, but each ring has the same weight now
1222        fMultV0->Fit(fpol0,"","",0,31);
1223        for(Int_t i(0); i < 4; i++) cuts->SetV0Cpol(i, fpol0->GetParameter(0));
1224        fMultV0->Fit(fpol0,"","",32,64);
1225        for(Int_t i(0); i < 4; i++) cuts->SetV0Apol(i, fpol0->GetParameter(0));
1226     }
1227     // the parameters to weigh the vzero track cuts have been extracted now, 
1228     // so we can pass them to the current track cuts obect
1229     cuts->SetV0gainEqualisation(fMultV0);       // passed as a TH1
1230
1231     // step 2) reweight the q-vectors that will be  called by flow methods which use
1232     // subevents
1233     // underlying assumption is that subevent a uses VZEROA
1234     // and subevent b uses VZEROC
1235     for(Int_t iside=0;iside<2;iside++){
1236         for(Int_t icoord=0;icoord<2;icoord++){
1237             for(Int_t i=0;i  < 9;i++){
1238                 char namecont[100];
1239                 if(iside==0 && icoord==0)
1240                   snprintf(namecont,100,"hQxc2_%i",i);
1241                 else if(iside==1 && icoord==0)
1242                   snprintf(namecont,100,"hQxa2_%i",i);
1243                 else if(iside==0 && icoord==1)
1244                   snprintf(namecont,100,"hQyc2_%i",i);
1245                 else if(iside==1 && icoord==1)
1246                   snprintf(namecont,100,"hQya2_%i",i);
1247
1248                 cont = (AliOADBContainer*) foadb->Get(namecont);
1249                 if(!cont){
1250                     printf("OADB object %s is not available in the file\n",namecont);
1251                     return;     
1252                 }
1253         
1254                 if(!(cont->GetObject(run))){
1255                     printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
1256                     run = 137366;
1257                 }
1258
1259                 // after grabbing all the info, set the CORRECTION TERMS to
1260                 // the 2nd and 3rd order qsub-vectors
1261                 // we do this here for all centralities, so that subsequent events
1262                 // can grab the correction from these cached values
1263                 fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
1264                 fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
1265
1266                 //for v3
1267                 if(iside==0 && icoord==0)
1268                   snprintf(namecont,100,"hQxc3_%i",i);
1269                 else if(iside==1 && icoord==0)
1270                   snprintf(namecont,100,"hQxa3_%i",i);
1271                 else if(iside==0 && icoord==1)
1272                   snprintf(namecont,100,"hQyc3_%i",i);
1273                 else if(iside==1 && icoord==1)
1274                   snprintf(namecont,100,"hQya3_%i",i);
1275
1276                 cont = (AliOADBContainer*) foadb->Get(namecont);
1277                 if(!cont){
1278                     printf("OADB object %s is not available in the file\n",namecont);
1279                     return;     
1280                 }
1281                 
1282                 if(!(cont->GetObject(run))){
1283                     printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
1284                     run = 137366;
1285                 }
1286                 fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
1287                 fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
1288
1289             }
1290         }
1291     }
1292     // set the recentering style (might be switched back to -1 if recentering is disabeled)
1293     fApplyRecentering = 2010;
1294 }
1295 //_____________________________________________________________________________
1296 void AliFlowEvent::SetVZEROCalibrationForTrackCuts2011(AliFlowTrackCuts* cuts)
1297 {
1298     // load the vzero q-sub vectors
1299     if(!cuts->GetEvent() || !cuts->GetEvent()->GetEventplane()) return;       // coverity
1300     Double_t qxEPa = 0, qyEPa = 0;
1301     Double_t qxEPc = 0, qyEPc = 0;
1302     Double_t qxEPa3 = 0, qyEPa3 = 0;
1303     Double_t qxEPc3 = 0, qyEPc3 = 0;
1304
1305     // get the q-vectors from the header 
1306     cuts->GetEvent()->GetEventplane()->CalculateVZEROEventPlane(cuts->GetEvent(), 8, 2, qxEPa, qyEPa);
1307     cuts->GetEvent()->GetEventplane()->CalculateVZEROEventPlane(cuts->GetEvent(), 9, 2, qxEPc, qyEPc);
1308     cuts->GetEvent()->GetEventplane()->CalculateVZEROEventPlane(cuts->GetEvent(), 8, 3, qxEPa3, qyEPa3);
1309     cuts->GetEvent()->GetEventplane()->CalculateVZEROEventPlane(cuts->GetEvent(), 9, 3, qxEPc3, qyEPc3);
1310  
1311     // store the values temporarily. this may seem
1312     // inelegant, but we don't want to include
1313     // aliflowtrackcuts or alivevnet in get2qsub
1314
1315     // qx and qy for vzero a, second harmonc
1316     fMeanQ[0][1][0] = qxEPa;
1317     fMeanQ[0][1][1] = qyEPa;
1318     // qx and qx for vzero c, second harmonic
1319     fMeanQ[0][0][0] = qxEPc;
1320     fMeanQ[0][0][1] = qyEPc;
1321     // qx and qy for vzero a, third harmonic
1322     fMeanQv3[0][1][0] = qxEPa3;
1323     fMeanQv3[0][1][1] = qyEPa3;
1324     // qx and qy for vzero c, third harmonic
1325     fMeanQv3[0][0][0] = qxEPc3;
1326     fMeanQv3[0][0][1] = qyEPc3;
1327
1328 //_____________________________________________________________________________
1329