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