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