AliPhysics  d497547 (d497547)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskCRCZDC.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2008, 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 thce 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  * analysis task for CRC with ZDC *
18  * *
19  * author: Jacopo Margutti *
20  * (margutti@nikhef.nl) *
21  **********************************/
22 
23 #include "Riostream.h" //needed as include
24 #include "TChain.h"
25 #include "TList.h"
26 #include "TTree.h"
27 #include "TRandom3.h"
28 #include "TTimeStamp.h"
29 #include "TStopwatch.h"
30 #include "TProfile.h"
31 #include <TList.h>
32 #include <TH1F.h>
33 #include <TH2F.h>
34 #include <TH3D.h>
35 #include <TFile.h>
36 #include <TString.h>
37 #include <TCanvas.h>
38 #include <TParticle.h>
39 #include "TProfile2D.h"
40 #include "TProfile3D.h"
41 
42 #include "AliAnalysisManager.h"
43 #include "AliInputEventHandler.h"
44 #include "AliVEvent.h"
45 #include "AliESD.h"
46 #include "AliESDEvent.h"
47 #include "AliESDHeader.h"
48 #include "AliESDInputHandler.h"
49 #include "AliESDZDC.h"
50 #include "AliMultiplicity.h"
51 #include "AliAnalysisUtils.h"
52 #include "AliAODHandler.h"
53 #include "AliAODTrack.h"
54 #include "AliAODEvent.h"
55 #include "AliAODHeader.h"
56 #include "AliAODVertex.h"
57 #include "AliAODVZERO.h"
58 #include "AliAODZDC.h"
59 #include "AliAODMCHeader.h"
60 #include "AliMCEventHandler.h"
61 #include "AliMCEvent.h"
62 #include "AliHeader.h"
63 #include "AliVParticle.h"
64 #include "AliStack.h"
65 #include "AliAODMCParticle.h"
66 #include "AliAnalysisTaskSE.h"
67 #include "AliGenEventHeader.h"
68 #include "AliPhysicsSelectionTask.h"
69 #include "AliPhysicsSelection.h"
70 #include "AliBackgroundSelection.h"
71 #include "AliTriggerAnalysis.h"
72 #include "AliCentrality.h"
73 #include "AliAnalysisTaskCRCZDC.h"
74 #include "AliMultSelection.h"
75 #include "AliLumiTools.h"
76 
77 // ALICE Correction Framework
78 #include "AliCFManager.h"
79 
80 // Interface to Event generators to get Reaction Plane Angle
81 #include "AliGenCocktailEventHeader.h"
82 #include "AliGenPythiaEventHeader.h"
83 #include "AliGenHijingEventHeader.h"
84 #include "AliGenGeVSimEventHeader.h"
85 #include "AliGenEposEventHeader.h"
86 
87 // Interface to Load short life particles
88 #include "TObjArray.h"
89 #include "AliFlowCandidateTrack.h"
90 
91 // Interface to make the Flow Event Simple used in the flow analysis methods
92 #include "AliFlowEvent.h"
93 #include "AliFlowTrackCuts.h"
94 #include "AliFlowEventCuts.h"
95 #include "AliFlowCommonConstants.h"
96 
97 ClassImp(AliAnalysisTaskCRCZDC)
98 
99 //________________________________________________________________________
102 fAnalysisType(kAUTOMATIC),
103 fRPType(""),
104 fCFManager1(NULL),
105 fCFManager2(NULL),
106 fCutsEvent(NULL),
107 fCutsRP(NULL),
108 fCutsPOI(NULL),
109 fCutContainer(new TList()),
110 fQAList(NULL),
111 fMinMult(0),
112 fMaxMult(10000000),
113 fMinA(-1.0),
114 fMaxA(-0.01),
115 fMinB(0.01),
116 fMaxB(1.0),
117 fGenHeader(NULL),
118 fPythiaGenHeader(NULL),
119 fHijingGenHeader(NULL),
120 fFlowTrack(NULL),
121 fAnalysisUtil(NULL),
122 fQAon(kFALSE),
123 fLoadCandidates(kFALSE),
124 fNbinsMult(10000),
125 fNbinsPt(100),
126 fNbinsPhi(100),
127 fNbinsEta(200),
128 fNbinsQ(500),
129 fNbinsMass(1),
130 fMultMin(0.),
131 fMultMax(10000.),
132 fPtMin(0.),
133 fPtMax(10.),
134 fPhiMin(0.),
135 fPhiMax(TMath::TwoPi()),
136 fEtaMin(-5.),
137 fEtaMax(5.),
138 fQMin(0.),
139 fQMax(3.),
140 fMassMin(-1.),
141 fMassMax(0.),
142 fHistWeightvsPhiMin(0.),
143 fHistWeightvsPhiMax(3.),
144 fExcludedEtaMin(0.),
145 fExcludedEtaMax(0.),
146 fExcludedPhiMin(0.),
147 fExcludedPhiMax(0.),
148 fAfterburnerOn(kFALSE),
149 fNonFlowNumberOfTrackClones(0),
150 fV1(0.),
151 fV2(0.),
152 fV3(0.),
153 fV4(0.),
154 fV5(0.),
155 fDifferentialV2(0),
156 fFlowEvent(NULL),
157 fShuffleTracks(kFALSE),
158 fMyTRandom3(NULL),
159 fAnalysisInput(kAOD),
160 fIsMCInput(kFALSE),
161 fUseMCCen(kTRUE),
162 fRejectPileUp(kTRUE),
163 fRejectPileUpTight(kFALSE),
164 fResetNegativeZDC(kFALSE),
165 fCentrLowLim(0.),
166 fCentrUpLim(100.),
167 fCentrEstimator(kV0M),
168 fOutput(0x0),
169 fhZNCvsZNA(0x0),
170 fhZDCCvsZDCCA(0x0),
171 fhZNCvsZPC(0x0),
172 fhZNAvsZPA(0x0),
173 fhZNvsZP(0x0),
174 fhZNvsVZERO(0x0),
175 fhZDCvsVZERO(0x0),
176 fhAsymm(0x0),
177 fhZNAvsAsymm(0x0),
178 fhZNCvsAsymm(0x0),
179 fhZNCvscentrality(0x0),
180 fhZNAvscentrality(0x0),
181 fhZPCvscentrality(0x0),
182 fhZPAvscentrality(0x0),
183 fCRCnRun(0),
184 fZDCGainAlpha(0.395),
185 fDataSet(kAny),
186 fStack(0x0),
187 fCutTPC(kFALSE),
188 fCenDis(0x0),
189 fVZEROMult(0x0),
190 fMultSelection(0x0),
191 fPileUpCount(0x0),
192 fPileUpMultSelCount(0x0),
193 fMultTOFLowCut(0x0),
194 fMultTOFHighCut(0x0),
195 fUseTowerEq(kFALSE),
196 fTowerEqList(NULL),
197 fUseBadTowerCalib(kFALSE),
198 fBadTowerCalibList(NULL),
199 fVZEROGainEqList(NULL),
200 fVZEROQVecRecList(NULL),
201 fUseZDCSpectraCorr(kFALSE),
202 fCorrectPhiTracklets(kFALSE),
203 fZDCSpectraCorrList(NULL),
204 fSpectraMCList(NULL),
205 fTrackQAList(NULL),
206 fBadTowerStuffList(NULL),
207 fVZEROStuffList(NULL),
208 fVZEROGainEqHist(NULL),
209 fMinRingVZC(1),
210 fMaxRingVZC(4),
211 fMinRingVZA(5),
212 fMaxRingVZA(8),
213 fCachedRunNum(0),
214 fhZNSpectra(0x0),
215 fhZNSpectraCor(0x0),
216 fhZNSpectraPow(0x0),
217 fhZNBCCorr(0x0)
218 {
219  for(int i=0; i<5; i++){
220  fhZNCPM[i] = 0x0;
221  fhZNAPM[i] = 0x0;
222  }
223  for(int i=0; i<4; i++){
224  fhZNCPMQiPMC[i] = 0x0;
225  fhZNAPMQiPMC[i] = 0x0;
226  }
227  for(Int_t r=0; r<fCRCMaxnRun; r++) {
228  fRunList[r] = 0;
229  }
230  for(Int_t c=0; c<2; c++) {
231  for(Int_t i=0; i<5; i++) {
232  fTowerGainEq[c][i] = NULL;
233  }
234  }
235  for(Int_t c=0; c<100; c++) {
236  fBadTowerCalibHist[c] = NULL;
237  }
238  for (Int_t k=0; k<fkVZEROnHar; k++) {
239 // fVZEROQVectorRecQx[k] = NULL;
240 // fVZEROQVectorRecQy[k] = NULL;
241  fVZEROQVectorRecQxStored[k] = NULL;
242  fVZEROQVectorRecQyStored[k] = NULL;
243  for (Int_t t=0; t<fkVZEROnQAplots; t++) {
244  fVZEROQVectorRecFinal[k][t] = NULL;
245  }
246  }
247  for(Int_t i=0; i<8; i++) {
248  SpecCorMu1[i] = NULL;
249  SpecCorMu2[i] = NULL;
250  SpecCorSi[i] = NULL;
251  SpecCorAv[i] = NULL;
252  }
253  this->InitializeRunArrays();
254  fMyTRandom3 = new TRandom3(1);
255  gRandom->SetSeed(fMyTRandom3->Integer(65539));
256  for(Int_t j=0; j<2; j++) {
257  for(Int_t c=0; c<10; c++) {
258  fPtSpecGen[j][c] = NULL;
259  fPtSpecFB32[j][c] = NULL;
260  fPtSpecFB96[j][c] = NULL;
261  fPtSpecFB128[j][c] = NULL;
262  fPtSpecFB768[j][c] = NULL;
263  }
264  }
265  for (Int_t c=0; c<2; c++) {
266  fhZNCenDis[c] = NULL;
267  }
268  for(Int_t fb=0; fb<fKNFBs; fb++){
269  for (Int_t i=0; i<4; i++) {
270  fTrackQADCAxy[fb][i] = NULL;
271  fTrackQADCAz[fb][i] = NULL;
272  fTrackQApT[fb][i] = NULL;
273  fTrackQADphi[fb][i] = NULL;
274  fEbEQRe[fb][i] = NULL;
275  fEbEQIm[fb][i] = NULL;
276  fEbEQMu[fb][i] = NULL;
277  }
278  }
279 }
280 
281 //________________________________________________________________________
282 AliAnalysisTaskCRCZDC::AliAnalysisTaskCRCZDC(const char *name, TString RPtype, Bool_t on, UInt_t iseed, Bool_t bCandidates):
283 AliAnalysisTaskSE(name),
284 fAnalysisType(kAUTOMATIC),
285 fRPType(RPtype),
286 fCFManager1(NULL),
287 fCFManager2(NULL),
288 fCutsEvent(NULL),
289 fCutsRP(NULL),
290 fCutsPOI(NULL),
291 fCutContainer(new TList()),
292 fQAList(NULL),
293 fMinMult(0),
294 fMaxMult(10000000),
295 fMinA(-1.0),
296 fMaxA(-0.01),
297 fMinB(0.01),
298 fMaxB(1.0),
299 fQAon(on),
300 fLoadCandidates(bCandidates),
301 fNbinsMult(10000),
302 fNbinsPt(100),
303 fNbinsPhi(100),
304 fNbinsEta(200),
305 fNbinsQ(500),
306 fNbinsMass(1),
307 fMultMin(0.),
308 fMultMax(10000.),
309 fPtMin(0.),
310 fPtMax(10.),
311 fPhiMin(0.),
312 fPhiMax(TMath::TwoPi()),
313 fEtaMin(-5.),
314 fEtaMax(5.),
315 fQMin(0.),
316 fQMax(3.),
317 fMassMin(-1.),
318 fMassMax(0.),
319 fHistWeightvsPhiMin(0.),
320 fHistWeightvsPhiMax(3.),
321 fExcludedEtaMin(0.),
322 fExcludedEtaMax(0.),
323 fExcludedPhiMin(0.),
324 fExcludedPhiMax(0.),
325 fAfterburnerOn(kFALSE),
326 fNonFlowNumberOfTrackClones(0),
327 fV1(0.),
328 fV2(0.),
329 fV3(0.),
330 fV4(0.),
331 fV5(0.),
332 fDifferentialV2(0),
333 fFlowEvent(NULL),
334 fShuffleTracks(kFALSE),
335 fMyTRandom3(NULL),
336 fAnalysisInput(kAOD),
337 fIsMCInput(kFALSE),
338 fUseMCCen(kTRUE),
339 fRejectPileUp(kTRUE),
340 fRejectPileUpTight(kFALSE),
341 fResetNegativeZDC(kFALSE),
342 fCentrLowLim(0.),
343 fCentrUpLim(100.),
344 fCentrEstimator(kV0M),
345 fOutput(0x0),
346 fhZNCvsZNA(0x0),
347 fhZDCCvsZDCCA(0x0),
348 fhZNCvsZPC(0x0),
349 fhZNAvsZPA(0x0),
350 fhZNvsZP(0x0),
351 fhZNvsVZERO(0x0),
352 fhZDCvsVZERO(0x0),
353 fhAsymm(0x0),
354 fhZNAvsAsymm(0x0),
355 fhZNCvsAsymm(0x0),
356 fhZNCvscentrality(0x0),
357 fhZNAvscentrality(0x0),
358 fhZPCvscentrality(0x0),
359 fhZPAvscentrality(0x0),
360 fDataSet(kAny),
361 fCRCnRun(0),
362 fZDCGainAlpha(0.395),
363 fGenHeader(NULL),
364 fPythiaGenHeader(NULL),
365 fHijingGenHeader(NULL),
366 fFlowTrack(NULL),
367 fAnalysisUtil(NULL),
368 fStack(0x0),
369 fCutTPC(kFALSE),
370 fCenDis(0x0),
371 fVZEROMult(0x0),
372 fMultSelection(0x0),
373 fPileUpCount(0x0),
374 fPileUpMultSelCount(0x0),
375 fMultTOFLowCut(0x0),
376 fMultTOFHighCut(0x0),
377 fUseTowerEq(kFALSE),
378 fTowerEqList(NULL),
379 fUseBadTowerCalib(kFALSE),
380 fBadTowerCalibList(NULL),
381 fVZEROGainEqList(NULL),
382 fVZEROQVecRecList(NULL),
383 fUseZDCSpectraCorr(kFALSE),
384 fCorrectPhiTracklets(kFALSE),
385 fZDCSpectraCorrList(NULL),
386 fSpectraMCList(NULL),
387 fTrackQAList(NULL),
388 fBadTowerStuffList(NULL),
389 fVZEROStuffList(NULL),
390 fVZEROGainEqHist(NULL),
391 fMinRingVZC(1),
392 fMaxRingVZC(4),
393 fMinRingVZA(5),
394 fMaxRingVZA(8),
395 fCachedRunNum(0),
396 fhZNSpectra(0x0),
397 fhZNSpectraCor(0x0),
398 fhZNSpectraPow(0x0),
399 fhZNBCCorr(0x0)
400 {
401 
402  for(int i=0; i<5; i++){
403  fhZNCPM[i] = 0x0;
404  fhZNAPM[i] = 0x0;
405  }
406  for(int i=0; i<4; i++){
407  fhZNCPMQiPMC[i] = 0x0;
408  fhZNAPMQiPMC[i] = 0x0;
409  }
410  for(Int_t r=0; r<fCRCMaxnRun; r++) {
411  fRunList[r] = 0;
412  }
413  for(Int_t c=0; c<2; c++) {
414  for(Int_t i=0; i<5; i++) {
415  fTowerGainEq[c][i] = NULL;
416  }
417  }
418  for(Int_t c=0; c<100; c++) {
419  fBadTowerCalibHist[c] = NULL;
420  }
421  for (Int_t k=0; k<fkVZEROnHar; k++) {
422  // fVZEROQVectorRecQx[k] = NULL;
423  // fVZEROQVectorRecQy[k] = NULL;
424  fVZEROQVectorRecQxStored[k] = NULL;
425  fVZEROQVectorRecQyStored[k] = NULL;
426  for (Int_t t=0; t<fkVZEROnQAplots; t++) {
427  fVZEROQVectorRecFinal[k][t] = NULL;
428  }
429  }
430  for(Int_t i=0; i<8; i++) {
431  SpecCorMu1[i] = NULL;
432  SpecCorMu2[i] = NULL;
433  SpecCorSi[i] = NULL;
434  SpecCorAv[i] = NULL;
435  }
436  this->InitializeRunArrays();
437  fMyTRandom3 = new TRandom3(iseed);
438  gRandom->SetSeed(fMyTRandom3->Integer(65539));
439 
440  DefineInput(0, TChain::Class());
441  // Define output slots here
442  // Define here the flow event output
443  DefineOutput(1, AliFlowEventSimple::Class());
444  DefineOutput(2, TList::Class());
445 
446  for(Int_t j=0; j<2; j++) {
447  for(Int_t c=0; c<10; c++) {
448  fPtSpecGen[j][c] = NULL;
449  fPtSpecFB32[j][c] = NULL;
450  fPtSpecFB96[j][c] = NULL;
451  fPtSpecFB128[j][c] = NULL;
452  fPtSpecFB768[j][c] = NULL;
453  }
454  }
455  for (Int_t c=0; c<2; c++) {
456  fhZNCenDis[c] = NULL;
457  }
458  for(Int_t fb=0; fb<fKNFBs; fb++){
459  for (Int_t i=0; i<4; i++) {
460  fTrackQADCAxy[fb][i] = NULL;
461  fTrackQADCAz[fb][i] = NULL;
462  fTrackQApT[fb][i] = NULL;
463  fTrackQADphi[fb][i] = NULL;
464  fEbEQRe[fb][i] = NULL;
465  fEbEQIm[fb][i] = NULL;
466  fEbEQMu[fb][i] = NULL;
467  }
468  }
469 }
470 
471 //________________________________________________________________________
473 {
474  // Destructor
475  if(fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()){
476  delete fOutput; fOutput=0;
477  }
478  delete fMyTRandom3;
479  delete fFlowEvent;
480  delete fFlowTrack;
481  delete fCutsEvent;
482  if (fTowerEqList) delete fTowerEqList;
487  if (fAnalysisUtil) delete fAnalysisUtil;
488  if (fQAList) delete fQAList;
489  if (fCutContainer) fCutContainer->Delete(); delete fCutContainer;
490 }
491 
492 //________________________________________________________________________
494 {
495  for(Int_t r=0;r<fCRCMaxnRun;r++) {
496  fCRCQVecListRun[r] = NULL;
497  for(Int_t k=0;k<fCRCnTow;k++) {
498  fZNCTower[r][k] = NULL;
499  fZNATower[r][k] = NULL;
500  }
501 // fhZNSpectraRbR[r] = NULL;
502  }
503  // for(Int_t i=0;i<fnCen;i++) {
504  // fPtPhiEtaRbRFB128[r][i] = NULL;
505  // fPtPhiEtaRbRFB768[r][i] = NULL;
506  // }
507  // }
508 }
509 
510 //________________________________________________________________________
512 {
513  // Create the output containers
514 
515  //set the common constants
518  cc->SetNbinsPt(fNbinsPt);
519  cc->SetNbinsPhi(fNbinsPhi);
520  cc->SetNbinsEta(fNbinsEta);
521  cc->SetNbinsQ(fNbinsQ);
523  cc->SetMultMin(fMultMin);
524  cc->SetMultMax(fMultMax);
525  cc->SetPtMin(fPtMin);
526  cc->SetPtMax(fPtMax);
527  cc->SetPhiMin(fPhiMin);
528  cc->SetPhiMax(fPhiMax);
529  cc->SetEtaMin(fEtaMin);
530  cc->SetEtaMax(fEtaMax);
531  cc->SetQMin(fQMin);
532  cc->SetQMax(fQMax);
533  cc->SetMassMin(fMassMin);
534  cc->SetMassMax(fMassMax);
537 
538  fFlowEvent = new AliFlowEvent(20000);
539  fFlowTrack = new AliFlowTrack();
540 
541  //printf(" AliAnalysisTaskCRCZDC::UserCreateOutputObjects()\n\n");
542  fOutput = new TList();
543  fOutput->SetOwner(kTRUE);
544  //fOutput->SetName("output");
545 
546  if (fQAon) {
547  fQAList = new TList();
548  fQAList->SetOwner(kTRUE);
549  fQAList->SetName("AliFlowEventCuts QA");
550  if (fCutsEvent->GetQA()) fQAList->Add(fCutsEvent->GetQA()); //0
551  if (fCutsRP->GetQA()) fQAList->Add(fCutsRP->GetQA()); //1
552  if (fCutsPOI->GetQA())fQAList->Add(fCutsPOI->GetQA()); //2
553  fOutput->Add(fQAList);
554  }
555 
556  fVZEROStuffList = new TList();
557  fVZEROStuffList->SetOwner(kTRUE);
558  fVZEROStuffList->SetName("VZERO stuff");
559  fOutput->Add(fVZEROStuffList);
560 
561  fBadTowerStuffList = new TList();
562  fBadTowerStuffList->SetOwner(kTRUE);
563  fBadTowerStuffList->SetName("BadTowerCalib");
565 
566  fCenDis = new TH1F("fCenDis", "fCenDis", 100, 0., 100.);
567  fOutput->Add(fCenDis);
568  fPileUpCount = new TH1F("fPileUpCount", "fPileUpCount", 9, 0., 9.);
569  fPileUpCount->GetXaxis()->SetBinLabel(1,"plpMV");
570  fPileUpCount->GetXaxis()->SetBinLabel(2,"fromSPD");
571  fPileUpCount->GetXaxis()->SetBinLabel(3,"RefMultiplicityComb08");
572  fPileUpCount->GetXaxis()->SetBinLabel(4,"IncompleteDAQ");
573  fPileUpCount->GetXaxis()->SetBinLabel(5,"abs(V0M-CL1)>7.5");
574  fPileUpCount->GetXaxis()->SetBinLabel(6,"missingVtx");
575  fPileUpCount->GetXaxis()->SetBinLabel(7,"inconsistentVtx");
576  fPileUpCount->GetXaxis()->SetBinLabel(8,"multESDTPCDif");
577  fPileUpCount->GetXaxis()->SetBinLabel(9,"extraPileUpMultSel");
578  fOutput->Add(fPileUpCount);
579  fPileUpMultSelCount = new TH1F("fPileUpMultSelCount", "fPileUpMultSelCount", 8, 0., 8.);
580  fPileUpMultSelCount->GetXaxis()->SetBinLabel(1,"IsNotPileup");
581  fPileUpMultSelCount->GetXaxis()->SetBinLabel(2,"IsNotPileupMV");
582  fPileUpMultSelCount->GetXaxis()->SetBinLabel(3,"IsNotPileupInMultBins");
583  fPileUpMultSelCount->GetXaxis()->SetBinLabel(4,"InconsistentVertices");
584  fPileUpMultSelCount->GetXaxis()->SetBinLabel(5,"TrackletVsCluster");
585  fPileUpMultSelCount->GetXaxis()->SetBinLabel(6,"AsymmetricInVZERO");
586  fPileUpMultSelCount->GetXaxis()->SetBinLabel(7,"IncompleteDAQ");
587  fPileUpMultSelCount->GetXaxis()->SetBinLabel(8,"GoodVertex2016");
589 
590  fMultTOFLowCut = new TF1("fMultTOFLowCut", "[0]+[1]*x+[2]*x*x+[3]*x*x*x - 4.*([4]+[5]*x+[6]*x*x+[7]*x*x*x+[8]*x*x*x*x+[9]*x*x*x*x*x)", 0, 10000);
591  fMultTOFLowCut->SetParameters(-1.0178, 0.333132, 9.10282e-05, -1.61861e-08, 1.47848, 0.0385923, -5.06153e-05, 4.37641e-08, -1.69082e-11, 2.35085e-15);
592  fOutput->Add(fMultTOFLowCut);
593  fMultTOFHighCut = new TF1("fMultTOFHighCut", "[0]+[1]*x+[2]*x*x+[3]*x*x*x + 4.*([4]+[5]*x+[6]*x*x+[7]*x*x*x+[8]*x*x*x*x+[9]*x*x*x*x*x)", 0, 10000);
594  fMultTOFHighCut->SetParameters(-1.0178, 0.333132, 9.10282e-05, -1.61861e-08, 1.47848, 0.0385923, -5.06153e-05, 4.37641e-08, -1.69082e-11, 2.35085e-15);
595  fOutput->Add(fMultTOFHighCut);
596 
597  for(Int_t c=0; c<2; c++) {
598  for(Int_t i=0; i<5; i++) {
599  fTowerGainEq[c][i] = new TH1D();
600  fOutput->Add(fTowerGainEq[c][i]);
601  }
602  }
603 
604  for (Int_t c=0; c<2; c++) {
605  fhZNCenDis[c] = new TH3D(Form("fhZNCenDis[%d]",c), Form("fhZNCenDis[%d]",c), 100, 0., 100., 100, -2., 2. , 100., -2., 2.);
606  fOutput->Add(fhZNCenDis[c]);
607  }
608 
609  if(fBadTowerCalibList) {
610  for(Int_t c=0; c<100; c++) {
611  fBadTowerCalibHist[c] = (TH2D*)fBadTowerCalibList->FindObject(Form("TH2Resp[%d]",c));
613  }
614  }
615  if(fZDCSpectraCorrList) {
616  for(Int_t i=0; i<8; i++) {
617  SpecCorMu1[i] = (TH1D*)fZDCSpectraCorrList->FindObject(Form("SpecCorMu1[%d]",i));
618  fOutput->Add(SpecCorMu1[i]);
619  SpecCorMu2[i] = (TH1D*)fZDCSpectraCorrList->FindObject(Form("SpecCorMu2[%d]",i));
620  fOutput->Add(SpecCorMu2[i]);
621  SpecCorAv[i] = (TH1D*)fZDCSpectraCorrList->FindObject(Form("SpecCorAv[%d]",i));
622  fOutput->Add(SpecCorAv[i]);
623  SpecCorSi[i] = (TH1D*)fZDCSpectraCorrList->FindObject(Form("SpecCorSi[%d]",i));
624  fOutput->Add(SpecCorSi[i]);
625  }
626  }
627 
628  fhZNSpectra = new TH3D("fhZNSpectra","fhZNSpectra",100,0.,100.,8,0.,8.,1000,0.,1.E5);
629  fOutput->Add(fhZNSpectra);
630  fhZNSpectraCor = new TH3D("fhZNSpectraCor","fhZNSpectraCor",100,0.,100.,8,0.,8.,1000,0.,1.E5);
631  fOutput->Add(fhZNSpectraCor);
632  fhZNSpectraPow = new TH3D("fhZNSpectraPow","fhZNSpectraPow",100,0.,100.,8,0.,8.,1000,0.,TMath::Power(1.E5,fZDCGainAlpha));
633  fOutput->Add(fhZNSpectraPow);
634  fhZNBCCorr = new TH3D("fhZNBCCorr","fhZNBCCorr",100,0.,100.,500,0.,1.E5,500,0.,1.E5);
635  fOutput->Add(fhZNBCCorr);
636 
637  if(fAnalysisType == kMCAOD) {
638 
639  fSpectraMCList = new TList();
640  fSpectraMCList->SetOwner(kTRUE);
641  fSpectraMCList->SetName("Spectra");
642  fOutput->Add(fSpectraMCList);
643 
644  Double_t xmin[] = {0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.,1.2,1.4,1.6,1.8,2.,2.33,2.66,3.,3.5,4.,5.,6.,9.,20.};
645  for(Int_t j=0; j<2; j++) {
646  for(Int_t c=0; c<10; c++) {
647  fPtSpecGen[j][c] = new TH1F(Form("fPtSpecGen[%d][%d]",j,c), Form("fPtSpecGen[%d][%d]",j,c), 23, xmin);
648  fSpectraMCList->Add(fPtSpecGen[j][c]);
649  fPtSpecFB32[j][c] = new TH1F(Form("fPtSpecFB32[%d][%d]",j,c), Form("fPtSpecFB32[%d][%d]",j,c), 23, xmin);
650  fSpectraMCList->Add(fPtSpecFB32[j][c]);
651  fPtSpecFB96[j][c] = new TH1F(Form("fPtSpecFB96[%d][%d]",j,c), Form("fPtSpecFB96[%d][%d]",j,c), 23, xmin);
652  fSpectraMCList->Add(fPtSpecFB96[j][c]);
653  fPtSpecFB128[j][c] = new TH1F(Form("fPtSpecFB128[%d][%d]",j,c), Form("fPtSpecFB128[%d][%d]",j,c), 23, xmin);
654  fSpectraMCList->Add(fPtSpecFB128[j][c]);
655  fPtSpecFB768[j][c] = new TH1F(Form("fPtSpecFB768[%d][%d]",j,c), Form("fPtSpecFB768[%d][%d]",j,c), 23, xmin);
656  fSpectraMCList->Add(fPtSpecFB768[j][c]);
657  }
658  }
659  }
660 
661  fAnalysisUtil = new AliAnalysisUtils();
662  fAnalysisUtil->SetUseMVPlpSelection(kTRUE);
663  fAnalysisUtil->SetUseOutOfBunchPileUp(kTRUE);
664 
665  for(int i=0; i<5; i++){
666  char hname[20];
667  sprintf(hname,"hZNCPM%d",i);
668  fhZNCPM[i] = new TH1F(hname, hname, 200, -50., 140000);
669  fOutput->Add(fhZNCPM[i]);
670  //
671  sprintf(hname,"hZNAPM%d",i);
672  fhZNAPM[i] = new TH1F(hname, hname, 200, -50., 140000);
673  fOutput->Add(fhZNAPM[i]);
674  //
675  if(i<4){
676  //
677  char hnamenc[20];
678  sprintf(hnamenc, "hZNCPMQ%dPMC",i+1);
679  fhZNCPMQiPMC[i] = new TH1F(hnamenc, hnamenc, 100, 0., 1.);
680  fOutput->Add(fhZNCPMQiPMC[i]);
681  //
682  char hnamena[20];
683  sprintf(hnamena, "hZNAPMQ%dPMC",i+1);
684  fhZNAPMQiPMC[i] = new TH1F(hnamena, hnamena, 100, 0., 1.);
685  fOutput->Add(fhZNAPMQiPMC[i]);
686  }
687  }
688 
689  fhZNCvsZNA = new TH2F("hZNCvsZNA","hZNCvsZNA",200,-50.,2.E5,200,-50.,2.E5);
690  fOutput->Add(fhZNCvsZNA);
691  fhZDCCvsZDCCA = new TH2F("hZDCCvsZDCCA","hZDCCvsZDCCA",200,0.,2.8E5,200,0.,2.8E5);
692  fOutput->Add(fhZDCCvsZDCCA);
693  fhZNCvsZPC = new TH2F("hZNCvsZPC","hZNCvsZPC",200,-50.,0.8E5,200,-50.,2.E5);
694  fOutput->Add(fhZNCvsZPC);
695  fhZNAvsZPA = new TH2F("hZNAvsZPA","hZNAvsZPA",200,-50.,0.8E5,200,-50.,2.E5);
696  fOutput->Add(fhZNAvsZPA);
697  fhZNvsZP = new TH2F("hZNvsZP","hZNvsZP",200,-50.,1.6E5,200,-50.,4.E5);
698  fOutput->Add(fhZNvsZP);
699  fhZNvsVZERO = new TH2F("hZNvsVZERO","hZNvsVZERO",250,0.,35000.,200,0.,4.E5);
700  fOutput->Add(fhZNvsVZERO);
701  fhZDCvsVZERO = new TH2F("hZDCvsVZERO","hZDCvsVZERO",250,0.,35000.,250,0.,5.6E5);
702  fOutput->Add(fhZDCvsVZERO);
703 
704  fhAsymm = new TH1F("hAsymm" , "Asimmetry ",200,-1.,1.);
705  fOutput->Add(fhAsymm);
706  fhZNAvsAsymm = new TH2F("hZNAvsAsymm","ZNA vs. asymm.",200,-1.,1.,200,-50.,2.E5);
707  fOutput->Add(fhZNAvsAsymm);
708  fhZNCvsAsymm = new TH2F("hZNCvsAsymm","ZNC vs. asymm.",200,-1.,1.,200,-50.,2.E5);
709  fOutput->Add(fhZNCvsAsymm);
710 
711  fhZNCvscentrality = new TH2F("hZNCvscentrality","hZNCvscentrality",100,0.,100.,200,-50.,2.E5);
713  fhZNAvscentrality = new TH2F("hZNAvscentrality","hZNAvscentrality",100,0.,100.,200,-50.,2.E5);
715  fhZPCvscentrality = new TH2F("hZPCvscentrality","hZPCvscentrality",100,0.,100.,200,-50.,0.8E5);
717  fhZPAvscentrality = new TH2F("hZPAvscentrality","hZPAvscentrality",100,0.,100.,200,-50.,0.8E5);
719 
720  //********************************************************************
721 
722  Int_t dRun10h[] = {139510, 139507, 139505, 139503, 139465, 139438, 139437, 139360, 139329, 139328, 139314, 139310, 139309, 139173, 139107, 139105, 139038, 139037, 139036, 139029, 139028, 138872, 138871, 138870, 138837, 138732, 138730, 138666, 138662, 138653, 138652, 138638, 138624, 138621, 138583, 138582, 138579, 138578, 138534, 138469, 138442, 138439, 138438, 138396, 138364, 138275, 138225, 138201, 138197, 138192, 138190, 137848, 137844, 137752, 137751, 137724, 137722, 137718, 137704, 137693, 137692, 137691, 137686, 137685, 137639, 137638, 137608, 137595, 137549, 137546, 137544, 137541, 137539, 137531, 137530, 137443, 137441, 137440, 137439, 137434, 137432, 137431, 137430, 137366, 137243, 137236, 137235, 137232, 137231, 137230, 137162, 137161};
723 
724  Int_t dRun11h[] = {167902, 167903, 167915, 167920, 167985, 167987, 167988, 168066, 168068, 168069, 168076, 168104, 168105, 168107, 168108, 168115, 168212, 168310, 168311, 168322, 168325, 168341, 168342, 168361, 168362, 168458, 168460, 168461, 168464, 168467, 168511, 168512, 168514, 168777, 168826, 168984, 168988, 168992, 169035, 169040, 169044, 169045, 169091, 169094, 169099, 169138, 169143, 169144, 169145, 169148, 169156, 169160, 169167, 169238, 169411, 169415, 169417, 169418, 169419, 169420, 169475, 169498, 169504, 169506, 169512, 169515, 169550, 169553, 169554, 169555, 169557, 169586, 169587, 169588, 169590, 169591, 169835, 169837, 169838, 169846, 169855, 169858, 169859, 169923, 169956, 169965, 170027, 170036,170040, 170081, 170083, 170084, 170085, 170088, 170089, 170091, 170155, 170159, 170163, 170193, 170203, 170204, 170207, 170228, 170230, 170268, 170269, 170270, 170306, 170308, 170309, 170311, 170312, 170313, 170315, 170387, 170388, 170572, 170593};
725 
726  Int_t dRun15o[] = {244917, 244918, 244975, 244980, 244982, 244983, 245064, 245066, 245068, 246390, 246391, 246392, 246994, 246991, 246989, 246984, 246982, 246980, 246948, 246945, 246928, 246851, 246847, 246846, 246845, 246844, 246810, 246809, 246808, 246807, 246805, 246804, 246766, 246765, 246763, 246760, 246759, 246758, 246757, 246751, 246750, 246495, 246493, 246488, 246487, 246434, 246431, 246428, 246424, 246276, 246275, 246272, 246271, 246225, 246222, 246217, 246185, 246182, 246181, 246180, 246178, 246153, 246152, 246151, 246115, 246113, 246089, 246087, 246053, 246052, 246049, 246048, 246042, 246037, 246036, 246012, 246003, 246001, 245954, 245952, 245949, 245923, 245833, 245831, 245829, 245705, 245702, 245700, 245692, 245683};
727 
728  Int_t dRun15ov6[] = {244918, 244975, 244980, 244982, 244983, 245064, 245066, 245068, 246390, 246391, 246392, 246994, 246991, 246989, 246984, 246982, 246980, 246948, 246945, 246928, 246851, 246847, 246846, 246845, 246844, 246810, 246809, 246808, 246807, 246805, 246804, 246766, 246765, 246763, 246760, 246759, 246758, 246757, 246751, 246750, 246495, 246493, 246488, 246487, 246434, 246431, 246428, 246424, 246276, 246275, 246272, 246271, 246225, 246222, 246217, 246185, 246182, 246181, 246180, 246178, 246153, 246152, 246151, 246148, 246115, 246113, 246089, 246087, 246053, 246052, 246049, 246048, 246042, 246037, 246036, 246012, 246003, 246001, 245963, 245954, 245952, 245949, 245923, 245833, 245831, 245829, 245705, 245702, 245700, 245692, 245683};
729 
730  Int_t dRun15opidfix[] = {245145, 245146, 245151, 245152, 245231, 245232, 245259, 245343, 245345, 245346, 245347, 245349, 245353, 245396, 245397, 245401, 245407, 245409, 245441, 245446, 245450, 245454, 245496, 245497, 245501, 245504, 245505, 245507, 245535, 245540, 245542, 245543, 245544, 245545, 245554};
731 
732  if(fDataSet==k2010) {fCRCnRun=92;}
733  if(fDataSet==k2011) {fCRCnRun=119;}
734  if(fDataSet==k2015) {fCRCnRun=90;}
735  if(fDataSet==k2015v6) {fCRCnRun=91;}
736  if(fDataSet==k2015pidfix) {fCRCnRun=35;}
737  if(fDataSet==kAny) {fCRCnRun=1;}
738 
739  Int_t d=0;
740  for(Int_t r=0; r<fCRCnRun; r++) {
741  if(fDataSet==k2010) fRunList[d] = dRun10h[r];
742  if(fDataSet==k2011) fRunList[d] = dRun11h[r];
743  if(fDataSet==k2015) fRunList[d] = dRun15o[r];
744  if(fDataSet==k2015v6) fRunList[d] = dRun15ov6[r];
745  if(fDataSet==k2015pidfix) fRunList[d] = dRun15opidfix[r];
746  if(fDataSet==kAny) fRunList[d] = 1;
747  d++;
748  }
749 
750  fVZEROMult = new TProfile2D("fVZEROMult","fVZEROMult",fCRCnRun,0.,1.*fCRCnRun,64,0.,64.);
751  for (Int_t i=0; i<fCRCnRun; i++) {
752  fVZEROMult->GetXaxis()->SetBinLabel(i+1,Form("%d",fRunList[i]));
753  }
755 
756  if(fVZEROGainEqList) {
757  fVZEROGainEqHist = (TH2D*)fVZEROGainEqList->FindObject("VZEROEqGain");
759  }
760  if(fVZEROQVecRecList) {
761  for (Int_t k=0; k<fkVZEROnHar; k++) {
762  fVZEROQVectorRecQxStored[k] = (TProfile3D*)fVZEROQVecRecList->FindObject(Form("fVZEROQVectorRecQx[%d]",k));
763  fVZEROQVectorRecQxStored[k]->SetTitle(Form("fVZEROQVectorRecQxStored[%d]",k));
764  fVZEROQVectorRecQxStored[k]->SetName(Form("fVZEROQVectorRecQxStored[%d]",k));
766  fVZEROQVectorRecQyStored[k] = (TProfile3D*)fVZEROQVecRecList->FindObject(Form("fVZEROQVectorRecQy[%d]",k));
767  fVZEROQVectorRecQyStored[k]->SetTitle(Form("fVZEROQVectorRecQyStored[%d]",k));
768  fVZEROQVectorRecQyStored[k]->SetName(Form("fVZEROQVectorRecQyStored[%d]",k));
770  for (Int_t t=0; t<fkVZEROnQAplots; t++) {
771  fVZEROQVectorRecFinal[k][t] = new TProfile2D(Form("fVZEROQVectorRecFinal[%d][%d]",k,t),Form("fVZEROQVectorRecFinal[%d][%d]",k,t),fCRCnRun,0.,1.*fCRCnRun,100,0.,100.,"s");
772  fVZEROQVectorRecFinal[k][t]->Sumw2();
774  }
775  }
776  }
777 
778 // for (Int_t k=0; k<fkVZEROnHar; k++) {
779 // fVZEROQVectorRecQx[k] = new TProfile3D(Form("fVZEROQVectorRecQx[%d]",k),Form("fVZEROQVectorRecQx[%d]",k),fCRCnRun,0.,1.*fCRCnRun,100,0.,100.,8,0.,8.,"s");
780 // fVZEROQVectorRecQx[k]->Sumw2();
781 // fVZEROStuffList->Add(fVZEROQVectorRecQx[k]);
782 // fVZEROQVectorRecQy[k] = new TProfile3D(Form("fVZEROQVectorRecQy[%d]",k),Form("fVZEROQVectorRecQy[%d]",k),fCRCnRun,0.,1.*fCRCnRun,100,0.,100.,8,0.,8.,"s");
783 // fVZEROQVectorRecQy[k]->Sumw2();
784 // fVZEROStuffList->Add(fVZEROQVectorRecQy[k]);
785 // }
786 
787  // track QA
788 
789  if(fAnalysisType == kTrackQA) {
790 
791  fTrackQAList = new TList();
792  fTrackQAList->SetOwner(kTRUE);
793  fTrackQAList->SetName("TrackQA");
794  fOutput->Add(fTrackQAList);
795 
796  Int_t nBinsEta=10;
797  Double_t dBinsEta[]={-0.8,-0.64,-0.48,-0.32,-0.16,0.,0.16,0.32,0.48,0.64,0.8};
798  Int_t nBinsPt = 26;
799  Double_t dBinsPt[] = {0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.,1.2,1.4,1.6,1.8,2.,2.33,2.66,3.,3.5,4.,5.,6.,8.,10.,14.,20.,30.,50.};
800 
801  for(Int_t fb=0; fb<fKNFBs; fb++){
802  for (Int_t i=0; i<4; i++) {
803  // DCA
804  fTrackQADCAxy[fb][i] = new TH3D(Form("fTrackQADCAxy[%d][%d]",fb,i),";#eta;p_{T} [GeV/c];DCAxy [cm]",10,-0.8,0.8,10,0.,5.,480,-2.4,2.4);
805  fTrackQAList->Add(fTrackQADCAxy[fb][i]);
806  fTrackQADCAz[fb][i] = new TH3D(Form("fTrackQADCAz[%d][%d]",fb,i),";#eta;p_{T} [GeV/c];DCAz [cm]",10,-0.8,0.8,10,0.,5.,640,-3.2,3.2);
807  fTrackQAList->Add(fTrackQADCAz[fb][i]);
808  // pT
809  fTrackQApT[fb][i] = new TH2D(Form("fTrackQApT[%d][%d]",fb,i),";#eta;p_{T} [GeV/c]",nBinsEta,dBinsEta,nBinsPt,dBinsPt);
810  fTrackQAList->Add(fTrackQApT[fb][i]);
811  // Dphi
812  fTrackQADphi[fb][i] = new TProfile2D(Form("fTrackQADphi[%d][%d]",fb,i),";#eta;p_{T} [GeV/c];cos(#Delta#phi)",nBinsEta,dBinsEta,nBinsPt,dBinsPt);
813  fTrackQAList->Add(fTrackQADphi[fb][i]);
814  // EbE
815  fEbEQRe[fb][i] = new TH2D(Form("fEbEQRe[%d][%d]",fb,i),";#eta;p_{T} [GeV/c];QRe(EbE)",nBinsEta,dBinsEta,nBinsPt,dBinsPt);
816  fTrackQAList->Add(fEbEQRe[fb][i]);
817  fEbEQIm[fb][i] = new TH2D(Form("fEbEQIm[%d][%d]",fb,i),";#eta;p_{T} [GeV/c];QRe(EbE)",nBinsEta,dBinsEta,nBinsPt,dBinsPt);
818  fTrackQAList->Add(fEbEQIm[fb][i]);
819  fEbEQMu[fb][i] = new TH2D(Form("fEbEQMu[%d][%d]",fb,i),";#eta;p_{T} [GeV/c];QRe(EbE)",nBinsEta,dBinsEta,nBinsPt,dBinsPt);
820  fTrackQAList->Add(fEbEQMu[fb][i]);
821  }
822  }
823 
824  }
825 
826  // run-by-run stuff
827 
828  if(!fUseTowerEq) {
829  Double_t ptmin[] = {0.2,0.4,0.6,0.8,1.,1.2,1.4,1.8,2.2,3.,4.,6.,8.,12.,20.};
830  Double_t phimin[] = {0.,TMath::Pi()/8.,2*TMath::Pi()/8.,3*TMath::Pi()/8.,4*TMath::Pi()/8.,5*TMath::Pi()/8.,6*TMath::Pi()/8.,7*TMath::Pi()/8.,8*TMath::Pi()/8.,9*TMath::Pi()/8.,10*TMath::Pi()/8.,11*TMath::Pi()/8.,12*TMath::Pi()/8.,13*TMath::Pi()/8.,14*TMath::Pi()/8.,15*TMath::Pi()/8.,16*TMath::Pi()/8.};
831  Double_t etamin[] = {-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8};
832 
833  for(Int_t r=0;r<fCRCnRun;r++) {
834  fCRCQVecListRun[r] = new TList();
835  fCRCQVecListRun[r]->SetName(Form("Run %d",fRunList[r]));
836  fCRCQVecListRun[r]->SetOwner(kTRUE);
837  fOutput->Add(fCRCQVecListRun[r]);
838 
839  for(Int_t k=0;k<fCRCnTow;k++) {
840  fZNCTower[r][k] = new TProfile(Form("fZNCTower[%d][%d]",fRunList[r],k),Form("fZNCTower[%d][%d]",fRunList[r],k),100,0.,100.,"s");
841  fZNCTower[r][k]->Sumw2();
842  fCRCQVecListRun[r]->Add(fZNCTower[r][k]);
843  fZNATower[r][k] = new TProfile(Form("fZNATower[%d][%d]",fRunList[r],k),Form("fZNATower[%d][%d]",fRunList[r],k),100,0.,100.,"s");
844  fZNATower[r][k]->Sumw2();
845  fCRCQVecListRun[r]->Add(fZNATower[r][k]);
846  }
847 
848  // fhZNSpectraRbR[r] = new TH3D(Form("fhZNSpectraRbR[%d]",fRunList[r]),Form("fhZNSpectraRbR[%d]",fRunList[r]),50,0.,100.,8,0.,8.,100,0.,1.E5);
849  // fCRCQVecListRun[r]->Add(fhZNSpectraRbR[r]);
850 
851  // for(Int_t i=0;i<fnCen;i++) {
852  // fPtPhiEtaRbRFB128[r][i] = new TH3F(Form("fPtPhiEtaRbRFB128[%d][%d]",r,i),Form("fPtPhiEtaRbRFB128[%d][%d]",r,i),14, ptmin, 16, phimin, 16, etamin);
853  // fCRCQVecListRun[r]->Add(fPtPhiEtaRbRFB128[r][i]);
854  // fPtPhiEtaRbRFB768[r][i] = new TH3F(Form("fPtPhiEtaRbRFB768[%d][%d]",r,i),Form("fPtPhiEtaRbRFB768[%d][%d]",r,i),14, ptmin, 16, phimin, 16, etamin);
855  // fCRCQVecListRun[r]->Add(fPtPhiEtaRbRFB768[r][i]);
856  // }
857  }
858  }
859 
860  PostData(2, fOutput);
861 }
862 
863 //________________________________________________________________________
865 {
866  // Execute analysis for current event:
867  AliMCEvent* McEvent = MCEvent();
868  AliAODEvent *aod = dynamic_cast<AliAODEvent*>(InputEvent());
869  AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
870  // AliMultiplicity* myTracklets = NULL;
871  // AliESDPmdTrack* pmdtracks = NULL;
872  // int availableINslot=1;
873 
874  if (!(fCutsRP&&fCutsPOI&&fCutsEvent)) {
875  AliError("cuts not set");
876  return;
877  }
878 
879  Int_t RunBin=-1, bin=0;
880  Int_t RunNum = aod->GetRunNumber();
881  for(Int_t c=0;c<fCRCnRun;c++) {
882  if(fRunList[c]==RunNum) RunBin=bin;
883  else bin++;
884  }
885  if(RunBin==-1) return;
886  if(fDataSet==kAny) RunBin=0;
887 
888  //DEFAULT - automatically takes care of everything
890 
891  // get centrality
892  Double_t centrV0M=300, centrCL1=300, centrCL0=300, centrTRK=300;
894  centrV0M = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("V0M");
895  centrCL1 = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("CL1");
896  centrCL0 = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("CL0");
897  centrTRK = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("TRK");
898  } else {
899  fMultSelection = (AliMultSelection*) InputEvent()->FindListObject("MultSelection");
900  if(!fMultSelection) {
901  //If you get this warning (and lPercentiles 300) please check that the AliMultSelectionTask actually ran (before your task)
902  AliWarning("AliMultSelection object not found!");
903  } else {
904  centrV0M = fMultSelection->GetMultiplicityPercentile("V0M");
905  centrCL1 = fMultSelection->GetMultiplicityPercentile("CL1");
906  centrCL0 = fMultSelection->GetMultiplicityPercentile("CL0");
907  centrTRK = fMultSelection->GetMultiplicityPercentile("TRK");
908  }
909  }
910 
911  //check event cuts
912  if (InputEvent()) {
913  if(!fCutsEvent->IsSelected(InputEvent(),MCEvent())) return;
914  if(fRejectPileUp) {
915  Bool_t IsPileUp = SelectPileup(aod);
916  if(IsPileUp) return;
917  }
918  }
919 
920  //first attach all possible information to the cuts
921  fCutsRP->SetEvent( InputEvent(), MCEvent() ); //attach event
922  fCutsPOI->SetEvent( InputEvent(), MCEvent() );
923 
924  //then make the event
926 
927  if(fAnalysisType == kTracklets) {
928  // fill with tracklets
929  AliAODTracklets* anInputTracklets = (AliAODTracklets*)aod->GetTracklets();
930  Int_t multSPD = anInputTracklets->GetNumberOfTracklets();
931  Double_t BField = aod->GetMagneticField();
932  //loop over tracklets
933  for (Int_t itracklet=0; itracklet<multSPD; ++itracklet) {
934  Float_t thetaTr= anInputTracklets->GetTheta(itracklet);
935  Float_t phiTr= anInputTracklets->GetPhi(itracklet);
936  // calculate eta
937  Float_t etaTr = -TMath::Log(TMath::Tan(thetaTr/2.));
938  //make new AliFLowTrackSimple
939  AliFlowTrack* pTrack = new AliFlowTrack();
940  pTrack->SetPt(0.5);
941  pTrack->SetEta(etaTr);
942  // set charge: "according to Ruben, with positive magnetic field, the positive tracks rotate clockwise. Since the angle stored in AOD is the difference between the hit on the first and second layers of SPD, in positive mag field, positive delta_phi -> positive track charge."
943  Double_t DeltaPhi = anInputTracklets->GetDeltaPhi(itracklet);
944  if(BField>0. && DeltaPhi>0.) pTrack->SetCharge(1);
945  if(BField>0. && DeltaPhi<0.) pTrack->SetCharge(-1);
946  if(BField<0. && DeltaPhi>0.) pTrack->SetCharge(-1);
947  if(BField<0. && DeltaPhi<0.) pTrack->SetCharge(1);
948  // correction of phi
950  phiTr += 39./34.*DeltaPhi;
951  if (phiTr < 0.) phiTr += 2.*TMath::Pi();
952  if (phiTr > 2.*TMath::Pi()) phiTr -= 2.*TMath::Pi();
953  }
954  pTrack->SetPhi(phiTr);
955  //marking the particles as POI type 2:
957  pTrack->SetPOItype(2,kTRUE);
959  //Add the track to the flowevent
960  fFlowEvent->AddTrack(pTrack);
961  }
962  }
963 
965 
970  fFlowEvent->SetCentralityCL1(centrCL1);
971  fFlowEvent->SetCentralityTRK(centrTRK);
972  // fFlowEvent->SetNITSCL1(((AliVAODHeader*)aod->GetHeader())->GetNumberOfITSClusters(1));
973  Double_t SumV0=0.;
974  for(Int_t i=0; i<64; i++) {
975  if(std::isfinite(aod->GetVZEROEqMultiplicity(i))) SumV0 += aod->GetVZEROEqMultiplicity(i);
976  }
977  fFlowEvent->SetNITSCL1(SumV0);
978 
979  // set absolute orbit number
980  UInt_t period = aod->GetPeriodNumber();
981  UInt_t orbit24 = aod->GetOrbitNumber(); // wrapped down to 24 bits
982  if (period > 255) { // 8 bits
983  period = 255;
984  orbit24 = (1<<24) - 1;
985  }
986  if (orbit24 >= (1<<24)) { // 24 bits
987  period = 255;
988  orbit24 = (1<<24) - 1;
989  }
990  UInt_t orbit = period * (1<<24) + orbit24;
991  fFlowEvent->SetAbsOrbit(orbit);
992 
993  Double_t vtxpos[3]={0.,0.,0.};
994  vtxpos[0] = ((AliAODVertex*)aod->GetPrimaryVertex())->GetX();
995  vtxpos[1] = ((AliAODVertex*)aod->GetPrimaryVertex())->GetY();
996  vtxpos[2] = ((AliAODVertex*)aod->GetPrimaryVertex())->GetZ();
997  fFlowEvent->SetVertexPosition(vtxpos);
998 
999  if (McEvent && McEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(McEvent);
1000 
1001  // run-by-run QA
1002  // for(Int_t jTracks = 0; jTracks<aod->GetNumberOfTracks(); jTracks++){
1003  // AliAODTrack* track = (AliAODTrack*)aod->GetTrack(jTracks);
1004  // if(!track) continue;
1005  // // general kinematic & quality cuts
1006  // if (track->Pt() < .2 || track->Pt() > 20. || TMath::Abs(track->Eta()) > .8 || track->GetTPCNcls() < 70) continue;
1007  // if (track->TestFilterBit(128)) fPtPhiEtaRbRFB128[RunBin][CenBin]->Fill(track->Pt(),track->Phi(),track->Eta());
1008  // if (track->TestFilterBit(768)) fPtPhiEtaRbRFB768[RunBin][CenBin]->Fill(track->Pt(),track->Phi(),track->Eta());
1009  // }
1010  fCenDis->Fill(centrV0M);
1011 
1012  }
1013 
1014  if (fAnalysisType == kTrackQA) {
1015 
1016  // empty flow event
1017  fFlowEvent->ClearFast();
1018 
1019  // get centrality
1020  Double_t centrV0M=300;
1022  centrV0M = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("V0M");
1023  } else {
1024  fMultSelection = (AliMultSelection*) InputEvent()->FindListObject("MultSelection");
1025  if(!fMultSelection) {
1026  //If you get this warning (and lPercentiles 300) please check that the AliMultSelectionTask actually ran (before your task)
1027  AliWarning("AliMultSelection object not found!");
1028  } else {
1029  centrV0M = fMultSelection->GetMultiplicityPercentile("V0M");
1030  }
1031  }
1032  if(centrV0M<10. || centrV0M>40.) return;
1033 
1034  //check event cuts
1035  if (InputEvent()) {
1036  if(!fCutsEvent->IsSelected(InputEvent(),MCEvent())) return;
1037  if(fRejectPileUp) {
1038  Bool_t IsPileUp = SelectPileup(aod);
1039  if(IsPileUp) return;
1040  }
1041  }
1042 
1043  Double_t BField = aod->GetMagneticField();
1044 
1045  for(Int_t jTracks = 0; jTracks<aod->GetNumberOfTracks(); jTracks++){
1046 
1047  AliAODTrack* track = (AliAODTrack*)aod->GetTrack(jTracks);
1048  if(!track) continue;
1049 
1050  Double_t dPt = track->Pt();
1051  Double_t dEta = track->Eta();
1052  Double_t dPhi = track->Phi();
1053  Double_t dCharge = track->Charge();
1054  Int_t cw = 0;
1055  if(BField>0. && dCharge>0.) cw=0;
1056  if(BField>0. && dCharge<0.) cw=1;
1057  if(BField<0. && dCharge>0.) cw=2;
1058  if(BField<0. && dCharge<0.) cw=3;
1059 
1060  // general kinematic cuts
1061  if (dPt < .2 || dPt > 50. || TMath::Abs(dEta) > 0.8) continue;
1062 
1063  // cut on DCA
1064  Double_t DCAxy = track->DCA();
1065  Double_t DCAz = track->ZAtDCA();
1066  if (std::abs((Int_t)DCAxy)==999 || std::abs((Int_t)DCAz)==999) {
1067  // re-evaluate the dca as it seems to not be natively present
1068  // allowed only for tracks inside the beam pipe
1069  Double_t pos[3] = {-99., -99., -99.};
1070  track->GetPosition(pos);
1071  if(pos[0]*pos[0]+pos[1]*pos[1] <= 3.*3.) {
1072  AliAODTrack copy(*track); // stack copy
1073  Double_t b[2] = {-99., -99.};
1074  Double_t bCov[3] = {-99., -99., -99.};
1075  if(copy.PropagateToDCA(aod->GetPrimaryVertex(), aod->GetMagneticField(), 100., b, bCov)) {
1076  DCAxy = b[0];
1077  DCAz = b[1];
1078  }
1079  }
1080  }
1081  if(fabs(DCAxy)>2.4 || fabs(DCAz)>3.2) continue;
1082 
1083  // various cuts on TPC clusters
1084  if (track->GetTPCNcls() < 70) continue;
1085  Double_t chi2_per_tpc = track->Chi2perNDF();
1086  if (chi2_per_tpc < 0.1 || chi2_per_tpc > 4.) continue;
1087  Double_t fraction_shared_tpccls = 1.*track->GetTPCnclsS()/track->GetTPCncls();
1088  if (fraction_shared_tpccls > 0.4) continue;
1089 
1090  Int_t FBarray[4] = {32,96,128,768};
1091 
1092  // test filter bits
1093  for(Int_t fb=0; fb<fKNFBs; fb++){
1094  if (track->TestFilterBit(FBarray[fb])) {
1095  fTrackQADCAxy[fb][cw]->Fill(dEta,dPt,DCAxy);
1096  fTrackQADCAz[fb][cw]->Fill(dEta,dPt,DCAz);
1097  fTrackQApT[fb][cw]->Fill(dEta,dPt);
1098  fEbEQRe[fb][cw]->Fill(dEta,dPt,TMath::Cos(dPhi));
1099  fEbEQIm[fb][cw]->Fill(dEta,dPt,TMath::Sin(dPhi));
1100  fEbEQMu[fb][cw]->Fill(dEta,dPt);
1101  }
1102  }
1103  }
1104 
1105  // compute cos(#Delta#phi)
1106  for(Int_t bx=1; bx<=fEbEQRe[0][0]->GetXaxis()->GetNbins(); bx++) {
1107  for(Int_t by=1; by<=fEbEQRe[0][0]->GetYaxis()->GetNbins(); by++) {
1108 
1109  Double_t dEta = fEbEQRe[0][0]->GetXaxis()->GetBinCenter(bx);
1110  Double_t dPt = fEbEQRe[0][0]->GetYaxis()->GetBinCenter(by);
1111 
1112  for(Int_t fb=0; fb<fKNFBs; fb++){
1113  for(Int_t c=0; c<4; c++){
1114 
1115  Double_t QRe = fEbEQRe[fb][c]->GetBinContent(bx,by);
1116  Double_t QIm = fEbEQIm[fb][c]->GetBinContent(bx,by);
1117  Double_t M = 1.*fEbEQMu[fb][c]->GetBinContent(bx,by);
1118 
1119  if(M>1.) {
1120  Double_t c1 = (QRe*QRe+QIm*QIm-M)/(M*(M-1.));
1121  fTrackQADphi[fb][c]->Fill(dEta,dPt,c1);
1122  }
1123  }
1124  }
1125 
1126  }
1127  }
1128 
1129  // reset EbE histograms
1130  for(Int_t fb=0; fb<fKNFBs; fb++){
1131  for(Int_t c=0; c<4; c++){
1132  fEbEQRe[fb][c]->Reset();
1133  fEbEQIm[fb][c]->Reset();
1134  fEbEQMu[fb][c]->Reset();
1135  }
1136  }
1137 
1138  }
1139 
1140  if (fAnalysisType == kMCAOD) {
1141 
1142  //check event cuts
1143  if (InputEvent()) {
1144  if(!fCutsEvent->IsSelected(InputEvent(),MCEvent())) return;
1145  if(fRejectPileUp && fAnalysisUtil->IsPileUpEvent(InputEvent())) return;
1146  }
1147 
1148  fFlowEvent->ClearFast();
1149 
1150  if(!McEvent) {
1151  AliError("ERROR: Could not retrieve MCEvent");
1152  return;
1153  }
1154  fStack = (TClonesArray*)aod->FindListObject(AliAODMCParticle::StdBranchName());
1155  if(!fStack){
1156  AliError("ERROR: Could not retrieve MCStack");
1157  return;
1158  }
1159 
1160  // get centrality (from AliMultSelection or AliCentrality)
1161  Double_t centr = 300;
1163  fMultSelection = (AliMultSelection*)aod->FindListObject("MultSelection");
1164  if(!fMultSelection) {
1165  //If you get this warning (and lPercentiles 300) please check that the AliMultSelectionTask actually ran (before your task)
1166  AliWarning("AliMultSelection object not found!");
1167  } else {
1168  centr = fMultSelection->GetMultiplicityPercentile("V0M");
1169  }
1170  } else {
1171  centr = (((AliVAODHeader*)aod->GetHeader())->GetCentralityP())->GetCentralityPercentile("V0M");
1172  }
1173  // centrality bin
1174  if (centr<fCentrLowLim || centr>=fCentrUpLim ) return;
1175  Int_t CenBin = -1;
1176  CenBin = GetCenBin(centr);
1177  if(CenBin==-1) return;
1178  fCenDis->Fill(centr);
1179 
1180  // reconstructed
1181  for(Int_t jTracks = 0; jTracks<aod->GetNumberOfTracks(); jTracks++){
1182 
1183  AliAODTrack* track = (AliAODTrack*)aod->GetTrack(jTracks);
1184  if(!track) continue;
1185 
1186  // to select primaries
1187  Int_t lp = TMath::Abs(track->GetLabel());
1188 
1189  // general kinematic cuts
1190  if (track->Pt() < .2 || track->Pt() > 20. || TMath::Abs(track->Eta()) > 0.8) continue;
1191 
1192  // cut on DCA
1193  Double_t DCAxy = track->DCA();
1194  Double_t DCAz = track->ZAtDCA();
1195  if(fabs(DCAxy)>2.4 || fabs(DCAz)>3.2) continue;
1196 
1197  // various cuts on TPC clusters
1198  if (track->GetTPCNcls() < 70) continue;
1199  Double_t chi2_per_tpc = track->Chi2perNDF();
1200  if (chi2_per_tpc < 0.1 || chi2_per_tpc > 4.) continue;
1201  Double_t fraction_shared_tpccls = 1.*track->GetTPCnclsS()/track->GetTPCncls();
1202  if (fraction_shared_tpccls > 0.4) continue;
1203 
1204  // test filter bits
1205  if (((AliAODMCParticle*)fStack->At(lp))->IsPhysicalPrimary()) {
1206  if (track->TestFilterBit(32)) fPtSpecFB32[0][CenBin]->Fill(track->Pt());
1207  if (track->TestFilterBit(96)) fPtSpecFB96[0][CenBin]->Fill(track->Pt());
1208  if (track->TestFilterBit(128)) fPtSpecFB128[0][CenBin]->Fill(track->Pt());
1209  if (track->TestFilterBit(768)) fPtSpecFB768[0][CenBin]->Fill(track->Pt());
1210  } else {
1211  if (track->TestFilterBit(32)) fPtSpecFB32[1][CenBin]->Fill(track->Pt());
1212  if (track->TestFilterBit(96)) fPtSpecFB96[1][CenBin]->Fill(track->Pt());
1213  if (track->TestFilterBit(128)) fPtSpecFB128[1][CenBin]->Fill(track->Pt());
1214  if (track->TestFilterBit(768)) fPtSpecFB768[1][CenBin]->Fill(track->Pt());
1215  }
1216 
1217  }
1218 
1219  // generated (physical primaries)
1220 
1221  for(Int_t jTracks = 0; jTracks<fStack->GetEntriesFast(); jTracks++) {
1222  AliAODMCParticle *MCpart = (AliAODMCParticle*)fStack->At(jTracks);
1223  if (!MCpart) {
1224  printf("ERROR: Could not receive MC track %d\n", jTracks);
1225  continue;
1226  }
1227 
1228  // kinematic cuts
1229  if ( MCpart->Pt() < .2 || MCpart->Pt() > 20. || TMath::Abs(MCpart->Eta()) > .8 ) continue;
1230  // select charged primaries
1231  if ( MCpart->Charge() == 0. || !MCpart->IsPhysicalPrimary()) continue;
1232 
1233  fPtSpecGen[0][CenBin]->Fill(MCpart->Pt());
1234  }
1235 
1236  // fGenHeader = McEvent->GenEventHeader();
1237  // if(fGenHeader) fPythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(fGenHeader);
1238  // printf("#reconstructed : %d (rejected from cuts %d), #MC primaries : %d (rejected from cuts %d) \n",AODPOIs,AODbads,MCPrims,MCSecos);
1239  fFlowEvent->SetReferenceMultiplicity(aod->GetNumberOfTracks());
1240  fFlowEvent->SetCentrality(centr);
1241  // if (McEvent && McEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(McEvent);
1242  fFlowEvent->SetRun(RunNum);
1243  // printf("Run : %d, RefMult : %d, Cent : %f \n",fFlowEvent->GetRun(),fFlowEvent->GetReferenceMultiplicity(),fFlowEvent->GetCentrality());
1244  }
1245 
1246  if(fAnalysisType == kMCESD) {
1247 
1248  fFlowEvent->ClearFast();
1249 
1250  if(!esd) {
1251  AliError("ERROR: Could not retrieve ESDEvent");
1252  return;
1253  }
1254  if(!McEvent) {
1255  AliError("ERROR: Could not retrieve MCEvent");
1256  return;
1257  }
1258  AliStack* fStack = fMCEvent->Stack();
1259  if(!fStack) {
1260  AliError("ERROR: Could not retrieve MCStack");
1261  return;
1262  }
1263 
1264  AliESDVertex *vertex = (AliESDVertex*) esd->GetPrimaryVertex();
1265  if (!vertex) return;
1266  if (TMath::Abs(vertex->GetZ()) > 10. ) return;
1267  if (vertex->GetNContributors() < 1 ) return;
1268  AliCentrality *centrality = esd->GetCentrality();
1269  if (!centrality) return;
1270  Double_t centr = centrality->GetCentralityPercentile("V0M");
1271  if (centr<fCentrLowLim || centr>=fCentrUpLim ) return;
1272  Int_t CenBin = -1;
1273  if (centr>0. && centr<5.) CenBin=0;
1274  if (centr>5. && centr<10.) CenBin=1;
1275  if (centr>10. && centr<20.) CenBin=2;
1276  if (centr>20. && centr<30.) CenBin=3;
1277  if (centr>30. && centr<40.) CenBin=4;
1278  if (centr>40. && centr<50.) CenBin=5;
1279  if (centr>50. && centr<60.) CenBin=6;
1280  if (centr>60. && centr<70.) CenBin=7;
1281  if (centr>70. && centr<80.) CenBin=8;
1282  if (centr>80. && centr<90.) CenBin=9;
1283  if(CenBin==-1) return;
1284 
1285  //Generated
1286  Int_t MCPrims = 0;
1287  for ( Int_t i=0 ; i<fStack->GetNtrack() ; i++ ) {
1288 
1289  //Primaries Selection
1290  TParticle *particle = (TParticle*)fStack->Particle(i);
1291  if (!particle) continue;
1292  if (!fStack->IsPhysicalPrimary(i)) continue;
1293  if ( particle->GetPDG()->Charge() == 0.) continue;
1294 
1295  //Kinematic Cuts
1296  if ( particle->Pt()<0.2 || particle->Pt()>10. ) continue;
1297  if ( TMath::Abs(particle->Eta())>0.8 ) continue;
1298 
1299  fFlowTrack->SetPhi(particle->Phi());
1300  fFlowTrack->SetEta(particle->Eta());
1301  fFlowTrack->SetPt(particle->Pt());
1303  fFlowTrack->SetForRPSelection(kTRUE);
1305  fFlowTrack->SetForPOISelection(kFALSE);
1307  MCPrims++;
1308 
1309  fPtSpecGen[0][CenBin]->Fill(particle->Pt());
1310 
1311  }
1312 
1313  //Reconstructed
1314  Int_t ESDPrims = 0;
1315  for (Int_t i=0 ; i<esd->GetNumberOfTracks() ; i++) {
1316 
1317  //Get reconstructed track
1318  AliVTrack *vtrack = static_cast<AliVTrack*>(esd->GetTrack(i));
1319  AliESDtrack *track = dynamic_cast<AliESDtrack*>(vtrack);
1320  if (!track) continue;
1321 
1322  //Primaries selection
1323  Int_t lp = TMath::Abs(track->GetLabel());
1324  if (!fStack->IsPhysicalPrimary(lp)) continue;
1325  TParticle *particle = (TParticle*)fStack->Particle(lp);
1326  if (!particle) continue;
1327  if (particle->GetPDG()->Charge() == 0.) continue;
1328 
1329  // if(!fCutsPOI->PassesESDcuts(track)) continue;
1330 
1331  Bool_t pass = kTRUE;
1332 
1333  if(fCutTPC) {
1334  // printf("******* cutting TPC ******** \n");
1335  UShort_t ntpccls = track->GetTPCNcls();
1336  Double_t tpcchi2 = track->GetTPCchi2();
1337  if (tpcchi2<0.2 || tpcchi2 >=4.) {
1338  // printf("TPCchi2 : %e %e ",tpcchi2,track->GetTPCchi2Iter1());
1339  pass=kFALSE;
1340  }
1341  if (ntpccls < 70) {
1342  // printf("#TPCcluster : %u %u %u %u ",ntpccls,track->GetTPCNclsF(),track->GetTPCNclsFIter1(),track->GetTPCNclsIter1());
1343  pass=kFALSE;
1344  }
1345  }
1346 
1347  Float_t dcaxy=0.0;
1348  Float_t dcaz=0.0;
1349  track->GetImpactParameters(dcaxy,dcaz);
1350  if (dcaxy > 0.3 || dcaz > 0.3) {
1351  // printf("DCA : %e %e ",dcaxy,dcaz);
1352  pass=kFALSE;
1353  }
1354  if(!pass) continue;
1355 
1356  //Kinematic Cuts
1357  if ( track->Pt()<0.2 || track->Pt()>10. ) continue;
1358  if ( TMath::Abs(track->Eta())>0.8 ) continue;
1359 
1360  fFlowTrack->SetPhi(track->Phi());
1361  fFlowTrack->SetEta(track->Eta());
1362  fFlowTrack->SetPt(track->Pt());
1364  fFlowTrack->SetForRPSelection(kFALSE);
1368  ESDPrims++;
1369 
1370  }
1371 
1372  // printf("#reconstructed : %d , #MC primaries : %d \n",ESDPrims,MCPrims);
1373  fFlowEvent->SetReferenceMultiplicity(esd->GetNumberOfTracks());
1374  fFlowEvent->SetCentrality(centr);
1375  if (McEvent && McEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(McEvent);
1376  fFlowEvent->SetRun(esd->GetRunNumber());
1377  // printf("Run : %d, RefMult : %d, Cent : %f \n",fFlowEvent->GetRun(),fFlowEvent->GetReferenceMultiplicity(),fFlowEvent->GetCentrality());
1378 
1379  } // end of if(fAnalysisType == kMCESD)
1380 
1381  if(fAnalysisType == kMCkine) {
1382 
1383  fFlowEvent->ClearFast();
1384 
1385  AliInputEventHandler* McHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
1386  if(!McHandler) {
1387  AliError("ERROR: Could not retrieve MCtruthEventHandler");
1388  return;
1389  }
1390  McEvent = McHandler->MCEvent();
1391  if(!McEvent) {
1392  AliError("ERROR: Could not retrieve MC event");
1393  return;
1394  }
1395 
1396  Int_t nTracks = McEvent->GetNumberOfTracks();
1397  // Int_t nPrimTr = McEvent->GetNumberOfPrimaries();
1398 
1399  //loop over tracks
1400  for (Int_t itrkN=0; itrkN<nTracks; itrkN++) {
1401  //get input particle
1402  AliMCParticle* pParticle = dynamic_cast<AliMCParticle*>(McEvent->GetTrack(itrkN));
1403  if (!pParticle) continue;
1404 
1405  //check if track passes the cuts
1406  if (McEvent->IsPhysicalPrimary(itrkN) && pParticle->Charge()!=0) {
1407  fFlowTrack->Set(pParticle);
1409  fFlowTrack->SetForRPSelection(kTRUE);
1414  }
1415  }// for all tracks
1416 
1417  // if monte carlo event get reaction plane from monte carlo (depends on generator)
1418  if (McEvent && McEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(McEvent);
1419  // set reference multiplicity
1420  fFlowEvent->SetReferenceMultiplicity(McEvent->GetNumberOfTracks());
1421  // tag subevents
1423  // set centrality from impact parameter
1424  Double_t ImpPar=0., CenPer=0.;
1425  fGenHeader = McEvent->GenEventHeader();
1426  if(fGenHeader){
1427  fPythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(fGenHeader);
1428  if(fPythiaGenHeader) ImpPar = fPythiaGenHeader->GetImpactParameter();
1429  fHijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(fGenHeader);
1430  if(fHijingGenHeader) ImpPar = fHijingGenHeader->ImpactParameter();
1431  if(ImpPar) CenPer = 0.3859796743103508*pow(ImpPar,2.);
1432  if(CenPer>0. && CenPer<100.) fFlowEvent->SetCentrality(CenPer);
1433  else return;
1434  fFlowEvent->SetRun(1);
1435  }
1436 
1437  } // end of if(fAnalysisType == kMCkine)
1438 
1439  if (!fFlowEvent) return; //shuts up coverity
1440 
1441  //check final event cuts
1442  Int_t mult = fFlowEvent->NumberOfTracks();
1443  // AliInfo(Form("FlowEvent has %i tracks",mult));
1444  if (mult<fMinMult || mult>fMaxMult) {
1445  AliWarning("FlowEvent cut on multiplicity"); return;
1446  }
1447 
1448  //define dead zone
1450 
1453  if (fAfterburnerOn)
1454  {
1455  //if reaction plane not set from elsewhere randomize it before adding flow
1457  fFlowEvent->SetMCReactionPlaneAngle(gRandom->Uniform(0.0,TMath::TwoPi()));
1458 
1459  if(fDifferentialV2)
1461  else
1462  fFlowEvent->AddFlow(fV1,fV2,fV3,fV4,fV5); //add flow
1463  fFlowEvent->CloneTracks(fNonFlowNumberOfTrackClones); //add nonflow by cloning tracks
1464  }
1466 
1467  //tag subEvents
1469 
1470  //do we want to serve shullfed tracks to everybody?
1472 
1473  // associate the mother particles to their daughters in the flow event (if any)
1475 
1476  //fListHistos->Print();
1477  //fOutputFile->WriteObject(fFlowEvent,"myFlowEventSimple");
1478 
1479  //********************************************************************************************************************************
1480 
1482 
1483  // PHYSICS SELECTION
1484  AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1485  AliInputEventHandler *hdr = (AliInputEventHandler*)am->GetInputEventHandler();
1486 
1487  if(hdr->IsEventSelected() && AliVEvent::kAny) {
1488 
1489  Double_t centrperc = fFlowEvent->GetCentrality();
1490  Int_t cenb = (Int_t)centrperc;
1491 
1492  AliAODTracklets *trackl = aod->GetTracklets();
1493  Int_t nTracklets = trackl->GetNumberOfTracklets();
1494 
1495  // get VZERO data
1496  AliAODVZERO *vzeroAOD = aod->GetVZEROData();
1497  Double_t multV0A = vzeroAOD->GetMTotV0A();
1498  Double_t multV0C = vzeroAOD->GetMTotV0C();
1499 
1500  // set VZERO Q-vectors
1501  if(fDataSet==k2015 || fDataSet==k2015v6) {
1502  Int_t CachednRing = 1;
1503  Double_t QxTot[fkVZEROnHar] = {0.}, QyTot[fkVZEROnHar] = {0.};
1504  Double_t denom = 0.;
1505  Double_t V0TotQC[fkVZEROnHar][2] = {{0.}}, V0TotQA[fkVZEROnHar][2] = {{0.}};
1506  Double_t MultC[fkVZEROnHar] = {0.}, MultA[fkVZEROnHar] = {0.};
1507 
1508  for(Int_t i=0; i<64; i++) {
1509 
1510  // correct multiplicity per channel
1511  Double_t mult = vzeroAOD->GetMultiplicity(i);
1512  if(fVZEROGainEqHist) {
1513  Double_t EqFactor = fVZEROGainEqHist->GetBinContent(RunBin+1,i+1);
1514  if(EqFactor>0.) mult *= EqFactor;
1515  }
1516  fVZEROMult->Fill(RunBin+0.5,i+0.5,mult);
1517 
1518  // build Q-vector per ring
1519  Int_t nRing = (Int_t)i/8 + 1;
1520  Double_t ChPhi = TMath::PiOver4()*(0.5+i%8);
1521 
1522  if(i == 63) {
1523  for (Int_t k=0; k<fkVZEROnHar; k++) {
1524  QxTot[k] += mult*TMath::Cos((k+1.)*ChPhi);
1525  QyTot[k] += mult*TMath::Sin((k+1.)*ChPhi);
1526  }
1527  denom += mult;
1528  nRing++;
1529  }
1530 
1531  if(nRing!=CachednRing) {
1532  for (Int_t k=0; k<fkVZEROnHar; k++) {
1533  Double_t QxRec = QxTot[k]/denom;
1534  Double_t QyRec = QyTot[k]/denom;
1535  // store values for re-centering
1536  // fVZEROQVectorRecQx[k]->Fill(RunBin+0.5,centrperc,CachednRing-0.5,QxRec);
1537  // fVZEROQVectorRecQy[k]->Fill(RunBin+0.5,centrperc,CachednRing-0.5,QyRec);
1538  // do re-centering
1539  if(fVZEROQVectorRecQxStored[k]) {
1540  if(!std::isnan(fVZEROQVectorRecQxStored[k]->GetBinContent(fVZEROQVectorRecQxStored[k]->FindBin(RunBin+0.5,centrperc,CachednRing-0.5)))) QxRec -= fVZEROQVectorRecQxStored[k]->GetBinContent(fVZEROQVectorRecQxStored[k]->FindBin(RunBin+0.5,centrperc,CachednRing-0.5));
1541  if(!std::isnan(fVZEROQVectorRecQyStored[k]->GetBinContent(fVZEROQVectorRecQyStored[k]->FindBin(RunBin+0.5,centrperc,CachednRing-0.5)))) QyRec -= fVZEROQVectorRecQyStored[k]->GetBinContent(fVZEROQVectorRecQyStored[k]->FindBin(RunBin+0.5,centrperc,CachednRing-0.5));
1542  }
1543  // sum of Q-vectors over all rings (total V0 Q-vector)
1544  if (CachednRing >= fMinRingVZC && CachednRing <= fMaxRingVZC) {
1545  V0TotQC[k][0] += QxRec*denom;
1546  V0TotQC[k][1] += QyRec*denom;
1547  MultC[k] += denom;
1548  }
1549  if (CachednRing >= fMinRingVZA && CachednRing <= fMaxRingVZA) {
1550  V0TotQA[k][0] += QxRec*denom;
1551  V0TotQA[k][1] += QyRec*denom;
1552  MultA[k] += denom;
1553  }
1554  QxTot[k] = 0.;
1555  QyTot[k] = 0.;
1556  }
1557  denom = 0.;
1558  CachednRing = nRing;
1559  }
1560  for (Int_t k=0; k<fkVZEROnHar; k++) {
1561  QxTot[k] += mult*TMath::Cos((k+1.)*ChPhi);
1562  QyTot[k] += mult*TMath::Sin((k+1.)*ChPhi);
1563  }
1564  denom += mult;
1565  }
1566 
1567  for (Int_t k=0; k<fkVZEROnHar; k++) {
1568  if(MultC[k]>0. && MultA[k]>0.) {
1569  Double_t QCx = V0TotQC[k][0]/MultC[k], QCy = V0TotQC[k][1]/MultC[k], QAx = V0TotQA[k][0]/MultA[k], QAy = V0TotQA[k][1]/MultA[k];
1570  if(!std::isnan(QCx) && !std::isnan(QCy) && !std::isnan(QAx) && !std::isnan(QAy)) {
1571  fFlowEvent->SetV02Qsub(QCx,QCy,MultC[k],QAx,QAy,MultA[k],k+1);
1572  fVZEROQVectorRecFinal[k][0]->Fill(RunBin+0.5,centrperc,QCx);
1573  fVZEROQVectorRecFinal[k][1]->Fill(RunBin+0.5,centrperc,QCy);
1574  fVZEROQVectorRecFinal[k][2]->Fill(RunBin+0.5,centrperc,QAx);
1575  fVZEROQVectorRecFinal[k][3]->Fill(RunBin+0.5,centrperc,QAy);
1576  fVZEROQVectorRecFinal[k][4]->Fill(RunBin+0.5,centrperc,QCx*QAx);
1577  fVZEROQVectorRecFinal[k][5]->Fill(RunBin+0.5,centrperc,QCy*QAy);
1578  fVZEROQVectorRecFinal[k][6]->Fill(RunBin+0.5,centrperc,QCx*QAy);
1579  fVZEROQVectorRecFinal[k][7]->Fill(RunBin+0.5,centrperc,QCy*QAx);
1580  } else {
1581  fFlowEvent->SetV02Qsub(0.,0.,0.,0.,0.,0.,k+1);
1582  }
1583  } else {
1584  fFlowEvent->SetV02Qsub(0.,0.,0.,0.,0.,0.,k+1);
1585  }
1586  }
1587  }
1588 
1589 // AliAODForwardMult* aodForward = static_cast<AliAODForwardMult*>(aodEvent->FindListObject("Forward"));
1590 // const TH2D& d2Ndetadphi = aodForward->GetHistogram();
1591 // Int_t nEta = d2Ndetadphi.GetXaxis()->GetNbins();
1592 // Int_t nPhi = d2Ndetadphi.GetYaxis()->GetNbins();
1593 // Double_t ret = 0.;
1594 // // Loop over eta
1595 // for (Int_t iEta = 1; iEta <= nEta; iEta++) {
1596 // Int_t valid = d2Ndetadphi.GetBinContent(iEta, 0);
1597 // if (!valid) continue; // No data expected for this eta
1598 // // Loop over phi
1599 // for (Int_t iPhi = 1; i <= nPhi; i++) {
1600 // ret = d2Ndetadphi.GetBinContent(iEta, iPhi);
1601 // printf("eta %e phi %e : %e \n",d2Ndetadphi.GetXaxis()->GetBinCenter(iEta),d2Ndetadphi.GetYaxis()->GetBinCenter(iPhi),ret);
1602 // }
1603 // }
1604 
1605  AliAODZDC *aodZDC = aod->GetZDCData();
1606 
1607  const Double_t * towZNCraw = aodZDC->GetZNCTowerEnergy();
1608  const Double_t * towZNAraw = aodZDC->GetZNATowerEnergy();
1609 
1610  // Get centroid from ZDCs *******************************************************
1611 
1612  Double_t Enucl = (RunNum < 209122 ? 1380. : 2511.);
1613  Double_t xyZNC[2]={0.,0.}, xyZNA[2]={0.,0.};
1614  Double_t towZNC[5]={0.}, towZNA[5]={0.};
1615 
1616  Double_t ZNCcalib=1., ZNAcalib=1.;
1617  if(fUseTowerEq) {
1618  if(RunNum!=fCachedRunNum) {
1619  for(Int_t i=0; i<5; i++) {
1620  fTowerGainEq[0][i] = (TH1D*)(fTowerEqList->FindObject(Form("fZNCTower[%d][%d]",RunNum,i)));
1621  fTowerGainEq[1][i] = (TH1D*)(fTowerEqList->FindObject(Form("fZNATower[%d][%d]",RunNum,i)));
1622  }
1623  }
1624  for(Int_t i=0; i<5; i++) {
1625  if(fTowerGainEq[0][i]) towZNC[i] = towZNCraw[i]*fTowerGainEq[0][i]->GetBinContent(fTowerGainEq[0][i]->FindBin(centrperc));
1626  if(fTowerGainEq[1][i]) towZNA[i] = towZNAraw[i]*fTowerGainEq[1][i]->GetBinContent(fTowerGainEq[1][i]->FindBin(centrperc));
1627  if(fResetNegativeZDC) {
1628  if(towZNC[i]<0.) towZNC[i] = 0.;
1629  if(towZNA[i]<0.) towZNA[i] = 0.;
1630  }
1631  }
1632  } else {
1633  for(Int_t i=0; i<5; i++) {
1634  towZNC[i] = towZNCraw[i];
1635  towZNA[i] = towZNAraw[i];
1636  if(fResetNegativeZDC) {
1637  if(towZNC[i]<0.) towZNC[i] = 0.;
1638  if(towZNA[i]<0.) towZNA[i] = 0.;
1639  }
1640  fZNCTower[RunBin][i]->Fill(centrperc,towZNC[i]);
1641  fZNATower[RunBin][i]->Fill(centrperc,towZNA[i]);
1642  }
1643  }
1644 
1645  if(RunNum>=245829) towZNA[2] = 0.;
1646  Double_t zncEnergy=0., znaEnergy=0.;
1647  for(Int_t i=0; i<5; i++){
1648  zncEnergy += towZNC[i];
1649  znaEnergy += towZNA[i];
1650  }
1651  if(RunNum>=245829) znaEnergy *= 8./7.;
1652  fFlowEvent->SetZNCQ0(towZNC[0]);
1653  fFlowEvent->SetZNAQ0(towZNA[0]);
1654 
1655  Double_t energyZNC = ((AliVAODHeader*)aod->GetHeader())->GetZDCN1Energy();
1656  Double_t energyZNA = ((AliVAODHeader*)aod->GetHeader())->GetZDCN2Energy();
1657  fFlowEvent->SetZNCEnergy(energyZNC);
1658  fFlowEvent->SetZNAEnergy(energyZNA);
1659 
1660  Double_t energyZPC = ((AliVAODHeader*)aod->GetHeader())->GetZDCP1Energy();
1661  Double_t energyZPA = ((AliVAODHeader*)aod->GetHeader())->GetZDCP2Energy();
1662  fFlowEvent->SetZPCEnergy(energyZPC);
1663  fFlowEvent->SetZPAEnergy(energyZPA);
1664 
1665  const Double_t x[4] = {-1.75, 1.75, -1.75, 1.75};
1666  const Double_t y[4] = {-1.75, -1.75, 1.75, 1.75};
1667  Double_t numXZNC=0., numYZNC=0., denZNC=0., cZNC, wZNC, EZNC, SumEZNC=0.;
1668  Double_t numXZNA=0., numYZNA=0., denZNA=0., cZNA, wZNA, EZNA, SumEZNA=0., BadChOr;
1669  Bool_t fAllChONZNC=kTRUE, fAllChONZNA=kTRUE;
1670 
1671  if (fUseMCCen) {
1672  for(Int_t i=0; i<4; i++){
1673 
1674  // get energy
1675  EZNC = towZNC[i+1];
1676  fhZNSpectra->Fill(centrperc,i+0.5,EZNC);
1677 // fhZNSpectraRbR[RunBin]->Fill(centrperc,i+0.5,EZNC);
1678  if(fUseZDCSpectraCorr && EZNC>0.) {
1679  Double_t mu1 = SpecCorMu1[i]->Interpolate(centrperc);
1680  Double_t mu2 = SpecCorMu2[i]->Interpolate(centrperc);
1681  Double_t av = SpecCorAv[i]->Interpolate(centrperc);
1682  Double_t cor1 = SpecCorSi[i]->Interpolate(centrperc);
1683  EZNC = exp( (log(EZNC) - mu1 + mu2*cor1)/cor1 ) + av;
1684  fhZNSpectraCor->Fill(centrperc,i+0.5,EZNC);
1685  }
1686  if(fUseZDCSpectraCorr && EZNC<=0.) fAllChONZNC=kFALSE;
1687 
1688  SumEZNC += EZNC;
1689 
1690  // build centroid
1691  wZNC = TMath::Power(EZNC, fZDCGainAlpha);
1692  numXZNC += x[i]*wZNC;
1693  numYZNC += y[i]*wZNC;
1694  denZNC += wZNC;
1695  fhZNSpectraPow->Fill(centrperc,i+0.5,wZNC);
1696 
1697  // get energy
1698  if(fDataSet==k2015 || fDataSet==k2015v6) {
1699  if(i==1) {
1700  EZNA = towZNA[0]-towZNA[1]-towZNA[3]-towZNA[4];
1701  if(fUseBadTowerCalib && fBadTowerCalibHist[cenb]) {
1702  EZNA = GetBadTowerResp(EZNA, fBadTowerCalibHist[cenb]);
1703  }
1704  } else {
1705  EZNA = towZNA[i+1];
1706  }
1707  } else {
1708  EZNA = towZNA[i+1];
1709  }
1710  fhZNSpectra->Fill(centrperc,i+4.5,EZNA);
1711 // fhZNSpectraRbR[RunBin]->Fill(centrperc,i+4.5,EZNA);
1712  if(fUseZDCSpectraCorr && EZNA>0.) {
1713  Double_t mu1 = SpecCorMu1[i+4]->Interpolate(centrperc);
1714  Double_t mu2 = SpecCorMu2[i+4]->Interpolate(centrperc);
1715  Double_t av = SpecCorAv[i+4]->Interpolate(centrperc);
1716  Double_t cor1 = SpecCorSi[i+4]->Interpolate(centrperc);
1717  EZNA = exp( (log(EZNA) - mu1 + mu2*cor1)/cor1 ) + av;
1718  fhZNSpectraCor->Fill(centrperc,i+4.5,EZNA);
1719  }
1720  if(fUseZDCSpectraCorr && EZNA<=0.) fAllChONZNA=kFALSE;
1721  SumEZNA += EZNA;
1722 
1723  // build centroid
1724  wZNA = TMath::Power(EZNA, fZDCGainAlpha);
1725  numXZNA += x[i]*wZNA;
1726  numYZNA += y[i]*wZNA;
1727  denZNA += wZNA;
1728  fhZNSpectraPow->Fill(centrperc,i+4.5,wZNA);
1729  }
1730  // store distribution for unfolding
1731  if(RunNum<245829) {
1732  Double_t recoE = towZNA[0]-towZNA[1]-towZNA[3]-towZNA[4];
1733  Double_t trueE = towZNA[2];
1734  fhZNBCCorr->Fill(centrperc,trueE,recoE);
1735  }
1736  if(denZNC>0.){
1737  Double_t nSpecnC = SumEZNC/Enucl;
1738  cZNC = 1.89358-0.71262/(nSpecnC+0.71789);
1739  xyZNC[0] = cZNC*numXZNC/denZNC;
1740  xyZNC[1] = cZNC*numYZNC/denZNC;
1741  denZNC *= cZNC;
1742  }
1743  else{
1744  xyZNC[0] = xyZNC[1] = 0.;
1745  }
1746  if(denZNA>0.){
1747  Double_t nSpecnA = SumEZNA/Enucl;
1748  cZNA = 1.89358-0.71262/(nSpecnA+0.71789);
1749  xyZNA[0] = cZNA*numXZNA/denZNA;
1750  xyZNA[1] = cZNA*numYZNA/denZNA;
1751  denZNA *= cZNA;
1752  }
1753  else{
1754  xyZNA[0] = xyZNA[1] = 0.;
1755  }
1756  } else {
1757  for(Int_t i=0; i<4; i++) {
1758  if(towZNC[i+1]>0.) {
1759  wZNC = TMath::Power(towZNC[i+1], fZDCGainAlpha);
1760  numXZNC += x[i]*wZNC;
1761  numYZNC += y[i]*wZNC;
1762  denZNC += wZNC;
1763  }
1764  if(towZNA[i+1]>0.) {
1765  wZNA = TMath::Power(towZNA[i+1], fZDCGainAlpha);
1766  numXZNA += x[i]*wZNA;
1767  numYZNA += y[i]*wZNA;
1768  denZNA += wZNA;
1769  }
1770  }
1771  if(denZNC!=0) {
1772  xyZNC[0] = numXZNC/denZNC;
1773  xyZNC[1] = numYZNC/denZNC;
1774  }
1775  else{
1776  xyZNC[0] = xyZNC[1] = 999.;
1777  zncEnergy = 0.;
1778  }
1779  if(denZNA!=0) {
1780  xyZNA[0] = numXZNA/denZNA;
1781  xyZNA[1] = numYZNA/denZNA;
1782  }
1783  else{
1784  xyZNA[0] = xyZNA[1] = 999.;
1785  znaEnergy = 0.;
1786  }
1787  }
1788 
1789  if(!fAllChONZNC) denZNC=-1.;
1790  if(!fAllChONZNA) denZNA=-1.;
1791 
1792  if(denZNC>0. && pow(xyZNC[0]*xyZNC[0]+xyZNC[1]*xyZNC[1],0.5)>1.E-6) fhZNCenDis[0]->Fill(centrperc,xyZNC[0],xyZNC[1]);
1793  if(denZNA>0. && pow(xyZNA[0]*xyZNA[0]+xyZNA[1]*xyZNA[1],0.5)>1.E-6) fhZNCenDis[1]->Fill(centrperc,-xyZNA[0], xyZNA[1]);
1794 
1795  fFlowEvent->SetZDC2Qsub(xyZNC,denZNC,xyZNA,denZNA);
1796 
1797  // ******************************************************************************
1798 
1799  for(int i=0; i<5; i++){
1800  fhZNCPM[i]->Fill(towZNC[i]);
1801  if((i<4) && (towZNC[0]>0.)) fhZNCPMQiPMC[i]->Fill(towZNC[i+1]/towZNC[0]);
1802  }
1803  for(int i=0; i<5; i++){
1804  fhZNAPM[i]->Fill(towZNA[i]);
1805  if(((i<4) && towZNA[0]>0.)) fhZNAPMQiPMC[i]->Fill(towZNA[i+1]/towZNA[0]);
1806  }
1807 
1808  fhZNCvsZNA->Fill(energyZNA, energyZNC);
1809  fhZDCCvsZDCCA->Fill(energyZNA+energyZPA, energyZNC+energyZPC);
1810  fhZNCvsZPC->Fill(energyZPC, energyZNC);
1811  fhZNAvsZPA->Fill(energyZPA, energyZNA);
1812  fhZNvsZP->Fill(energyZPA+energyZPC, energyZNA+energyZNC);
1813  fhZNvsVZERO->Fill(multV0A+multV0C, energyZNC+energyZNA);
1814  fhZDCvsVZERO->Fill(multV0A+multV0C, energyZNA+energyZPA+energyZNC+energyZPC);
1815 
1816  Double_t asymmetry = -999.;
1817  if((energyZNC+energyZNA)>0.) asymmetry = (energyZNC-energyZNA)/(energyZNC+energyZNA);
1818  fhAsymm->Fill(asymmetry);
1819  fhZNAvsAsymm->Fill(asymmetry, energyZNA);
1820  fhZNCvsAsymm->Fill(asymmetry, energyZNC);
1821 
1822  fhZNCvscentrality->Fill(centrperc, energyZNC);
1823  fhZNAvscentrality->Fill(centrperc, energyZNA);
1824  fhZPCvscentrality->Fill(centrperc, energyZPC);
1825  fhZPAvscentrality->Fill(centrperc, energyZPA);
1826 
1827  } // PHYSICS SELECTION
1828 
1829  }
1830 
1831  // p) cache run number
1833 
1834  // printf("debug: NoRPs %e, NoPOIs %e, RunNum %d, Cen %e \n",fFlowEvent->GetNumberOfRPs(),fFlowEvent->GetNumberOfPOIs(),fCachedRunNum,fFlowEvent->GetCentrality());
1835 
1836  PostData(1, fFlowEvent);
1837 
1838  PostData(2, fOutput);
1839 }
1840 //________________________________________________________________________
1841 
1843 {
1844  Bool_t BisPileup=kFALSE;
1845  Double_t centrV0M=300., centrCL1=300.;
1846 
1848 
1849  // pileup for LHC10h and LHC11h
1850 
1851  centrV0M = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("V0M");
1852  centrCL1 = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("CL1");
1853 
1854  // check anyway pileup
1855  if (plpMV(aod)) {
1856  fPileUpCount->Fill(0.5);
1857  BisPileup=kTRUE;
1858  }
1859 
1860  Short_t isPileup = aod->IsPileupFromSPD(3);
1861  if (isPileup != 0) {
1862  fPileUpCount->Fill(1.5);
1863  BisPileup=kTRUE;
1864  }
1865 
1866  if (((AliAODHeader*)aod->GetHeader())->GetRefMultiplicityComb08() < 0) {
1867  fPileUpCount->Fill(2.5);
1868  BisPileup=kTRUE;
1869  }
1870 
1871  if (aod->IsIncompleteDAQ()) {
1872  fPileUpCount->Fill(3.5);
1873  BisPileup=kTRUE;
1874  }
1875 
1876  // check vertex consistency
1877  const AliAODVertex* vtTrc = aod->GetPrimaryVertex();
1878  const AliAODVertex* vtSPD = aod->GetPrimaryVertexSPD();
1879 
1880  if (vtTrc->GetNContributors() < 2 || vtSPD->GetNContributors()<1) {
1881  fPileUpCount->Fill(5.5);
1882  BisPileup=kTRUE;
1883  }
1884 
1885  double covTrc[6], covSPD[6];
1886  vtTrc->GetCovarianceMatrix(covTrc);
1887  vtSPD->GetCovarianceMatrix(covSPD);
1888 
1889  double dz = vtTrc->GetZ() - vtSPD->GetZ();
1890 
1891  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
1892  double errTrc = TMath::Sqrt(covTrc[5]);
1893  double nsigTot = dz/errTot;
1894  double nsigTrc = dz/errTrc;
1895 
1896  if (TMath::Abs(dz)>0.2 || TMath::Abs(nsigTot)>10 || TMath::Abs(nsigTrc)>20) {
1897  fPileUpCount->Fill(6.5);
1898  BisPileup=kTRUE;
1899  }
1900 
1901  if (fAnalysisUtil->IsPileUpEvent(InputEvent())) {
1902  fPileUpCount->Fill(7.5);
1903  BisPileup=kTRUE;
1904  }
1905 
1906  }
1907  else {
1908 
1909  // pileup for LHC15o, using AliMultSelection
1910 
1911  if(fMultSelection) {
1912  centrV0M = fMultSelection->GetMultiplicityPercentile("V0M");
1913  centrCL1 = fMultSelection->GetMultiplicityPercentile("CL1");
1914  } else {
1915  BisPileup=kTRUE;
1916  }
1917 
1918  // pileup from AliMultSelection
1919  if(!fMultSelection->GetThisEventIsNotPileup()) fPileUpMultSelCount->Fill(0.5);
1920  if(!fMultSelection->GetThisEventIsNotPileupMV()) fPileUpMultSelCount->Fill(1.5);
1921  if(!fMultSelection->GetThisEventIsNotPileupInMultBins()) fPileUpMultSelCount->Fill(2.5);
1922  if(!fMultSelection->GetThisEventHasNoInconsistentVertices()) fPileUpMultSelCount->Fill(3.5);
1923  if(!fMultSelection->GetThisEventPassesTrackletVsCluster()) fPileUpMultSelCount->Fill(4.5);
1924  if(!fMultSelection->GetThisEventIsNotAsymmetricInVZERO()) fPileUpMultSelCount->Fill(5.5);
1925  if(!fMultSelection->GetThisEventIsNotIncompleteDAQ()) fPileUpMultSelCount->Fill(6.5);
1926  if(!fMultSelection->GetThisEventHasGoodVertex2016()) fPileUpMultSelCount->Fill(7.5);
1927 
1928  // pile-up a la Dobrin for LHC15o
1929  if (plpMV(aod)) {
1930  fPileUpCount->Fill(0.5);
1931  BisPileup=kTRUE;
1932  }
1933 
1934  Short_t isPileup = aod->IsPileupFromSPD(3);
1935  if (isPileup != 0) {
1936  fPileUpCount->Fill(1.5);
1937  BisPileup=kTRUE;
1938  }
1939 
1940  if (((AliAODHeader*)aod->GetHeader())->GetRefMultiplicityComb08() < 0) {
1941  fPileUpCount->Fill(2.5);
1942  BisPileup=kTRUE;
1943  }
1944 
1945  if (aod->IsIncompleteDAQ()) {
1946  fPileUpCount->Fill(3.5);
1947  BisPileup=kTRUE;
1948  }
1949 
1950  if(fabs(centrV0M-centrCL1)>7.5) {
1951  fPileUpCount->Fill(4.5);
1952  BisPileup=kTRUE;
1953  }
1954 
1955  // check vertex consistency
1956  const AliAODVertex* vtTrc = aod->GetPrimaryVertex();
1957  const AliAODVertex* vtSPD = aod->GetPrimaryVertexSPD();
1958 
1959  if (vtTrc->GetNContributors() < 2 || vtSPD->GetNContributors()<1) {
1960  fPileUpCount->Fill(5.5);
1961  BisPileup=kTRUE;
1962  }
1963 
1964  double covTrc[6], covSPD[6];
1965  vtTrc->GetCovarianceMatrix(covTrc);
1966  vtSPD->GetCovarianceMatrix(covSPD);
1967 
1968  double dz = vtTrc->GetZ() - vtSPD->GetZ();
1969 
1970  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
1971  double errTrc = TMath::Sqrt(covTrc[5]);
1972  double nsigTot = dz/errTot;
1973  double nsigTrc = dz/errTrc;
1974 
1975  if (TMath::Abs(dz)>0.2 || TMath::Abs(nsigTot)>10 || TMath::Abs(nsigTrc)>20) {
1976  fPileUpCount->Fill(6.5);
1977  BisPileup=kTRUE;
1978  }
1979 
1980  // cuts on tracks
1981  const Int_t nTracks = aod->GetNumberOfTracks();
1982  Int_t multEsd = ((AliAODHeader*)aod->GetHeader())->GetNumberOfESDTracks();
1983 
1984  Int_t multTrk = 0;
1985  Int_t multTrkBefC = 0;
1986  Int_t multTrkTOFBefC = 0;
1987  Int_t multTPC = 0;
1988 
1989  for (Int_t it = 0; it < nTracks; it++) {
1990 
1991  AliAODTrack* aodTrk = (AliAODTrack*)aod->GetTrack(it);
1992  if (!aodTrk){
1993  delete aodTrk;
1994  continue;
1995  }
1996 
1997 // if (aodTrk->TestFilterBit(32)){
1998 // multTrkBefC++;
1999 //
2000 // if ( TMath::Abs(aodTrk->GetTOFsignalDz()) <= 10. && aodTrk->GetTOFsignal() >= 12000. && aodTrk->GetTOFsignal() <= 25000.)
2001 // multTrkTOFBefC++;
2002 //
2003 // if ((TMath::Abs(aodTrk->Eta()) < 0.8) && (aodTrk->GetTPCNcls() >= 70) && (aodTrk->Pt() >= 0.2) && (aodTrk->Pt() < 20.))
2004 // multTrk++;
2005 // }
2006 
2007  if (aodTrk->TestFilterBit(128))
2008  multTPC++;
2009 
2010  } // end of for (Int_t it = 0; it < nTracks; it++)
2011 
2012  Double_t multTPCn = multTPC;
2013  Double_t multEsdn = multEsd;
2014  Double_t multESDTPCDif = multEsdn - multTPCn*3.38;
2015 
2016  if (multESDTPCDif > (fRejectPileUpTight?700.:15000.)) {
2017  fPileUpCount->Fill(7.5);
2018  BisPileup=kTRUE;
2019  }
2020 
2021  if(fRejectPileUpTight) {
2022  if(BisPileup==kFALSE) {
2023  if(!fMultSelection->GetThisEventIsNotPileup()) BisPileup=kTRUE;
2024  if(!fMultSelection->GetThisEventIsNotPileupMV()) BisPileup=kTRUE;
2025  if(!fMultSelection->GetThisEventIsNotPileupInMultBins()) BisPileup=kTRUE;
2026  if(!fMultSelection->GetThisEventHasNoInconsistentVertices()) BisPileup=kTRUE;
2027  if(!fMultSelection->GetThisEventPassesTrackletVsCluster()) BisPileup=kTRUE;
2028  if(!fMultSelection->GetThisEventIsNotIncompleteDAQ()) BisPileup=kTRUE;
2029  if(!fMultSelection->GetThisEventHasGoodVertex2016()) BisPileup=kTRUE;
2030  if(BisPileup) fPileUpCount->Fill(8.5);
2031  }
2032  }
2033  }
2034 
2035  return BisPileup;
2036 }
2037 
2038 //________________________________________________________________________
2039 
2041 {
2042  Double_t EtC = BadTowerCalibHist->ProjectionY("",BadTowerCalibHist->GetXaxis()->FindBin(Et),BadTowerCalibHist->GetXaxis()->FindBin(Et))->GetRandom();
2043  return EtC;
2044 }
2045 
2046 //________________________________________________________________________
2047 
2049 {
2050  Int_t CenBin=-1;
2051  if (Centrality>0. && Centrality<5.) CenBin=0;
2052  if (Centrality>5. && Centrality<10.) CenBin=1;
2053  if (Centrality>10. && Centrality<20.) CenBin=2;
2054  if (Centrality>20. && Centrality<30.) CenBin=3;
2055  if (Centrality>30. && Centrality<40.) CenBin=4;
2056  if (Centrality>40. && Centrality<50.) CenBin=5;
2057  if (Centrality>50. && Centrality<60.) CenBin=6;
2058  if (Centrality>60. && Centrality<70.) CenBin=7;
2059  if (Centrality>70. && Centrality<80.) CenBin=8;
2060  if (Centrality>80. && Centrality<90.) CenBin=9;
2061  if (CenBin>=fnCen) CenBin=-1;
2062  if (fnCen==1) CenBin=0;
2063  return CenBin;
2064 } // end of AliFlowAnalysisCRC::GetCRCCenBin(Double_t Centrality)
2065 //_____________________________________________________________________________
2066 
2067 Double_t AliAnalysisTaskCRCZDC::GetWDist(const AliVVertex* v0, const AliVVertex* v1)
2068 {
2069  // calculate sqrt of weighted distance to other vertex
2070  if (!v0 || !v1) {
2071  printf("One of vertices is not valid\n");
2072  return 0;
2073  }
2074  static TMatrixDSym vVb(3);
2075  double dist = -1;
2076  double dx = v0->GetX()-v1->GetX();
2077  double dy = v0->GetY()-v1->GetY();
2078  double dz = v0->GetZ()-v1->GetZ();
2079  double cov0[6],cov1[6];
2080  v0->GetCovarianceMatrix(cov0);
2081  v1->GetCovarianceMatrix(cov1);
2082  vVb(0,0) = cov0[0]+cov1[0];
2083  vVb(1,1) = cov0[2]+cov1[2];
2084  vVb(2,2) = cov0[5]+cov1[5];
2085  vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
2086  vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
2087  vVb.InvertFast();
2088  if (!vVb.IsValid()) {printf("Singular Matrix\n"); return dist;}
2089  dist = vVb(0,0)*dx*dx + vVb(1,1)*dy*dy + vVb(2,2)*dz*dz
2090  + 2*vVb(0,1)*dx*dy + 2*vVb(0,2)*dx*dz + 2*vVb(1,2)*dy*dz;
2091  return dist>0 ? TMath::Sqrt(dist) : -1;
2092 }
2093 //________________________________________________________________________
2094 
2096 {
2097  // check for multi-vertexer pile-up
2098 
2099  const int kMinPlpContrib = 5;
2100  const double kMaxPlpChi2 = 5.0;
2101  const double kMinWDist = 15;
2102 
2103  const AliVVertex* vtPrm = 0;
2104  const AliVVertex* vtPlp = 0;
2105  int nPlp = 0;
2106 
2107  if ( !(nPlp=aod->GetNumberOfPileupVerticesTracks()) ) return kFALSE;
2108  vtPrm = aod->GetPrimaryVertex();
2109  if (vtPrm == aod->GetPrimaryVertexSPD()) return kTRUE; // there are pile-up vertices but no primary
2110 
2111  //int bcPrim = vtPrm->GetBC();
2112 
2113  for (int ipl=0;ipl<nPlp;ipl++) {
2114  vtPlp = (const AliVVertex*)aod->GetPileupVertexTracks(ipl);
2115  //
2116  if (vtPlp->GetNContributors() < kMinPlpContrib) continue;
2117  if (vtPlp->GetChi2perNDF() > kMaxPlpChi2) continue;
2118  // int bcPlp = vtPlp->GetBC();
2119  // if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2) return kTRUE; // pile-up from other BC
2120  //
2121  double wDst = GetWDist(vtPrm,vtPlp);
2122  if (wDst<kMinWDist) continue;
2123  //
2124  return kTRUE; // pile-up: well separated vertices
2125  }
2126 
2127  return kFALSE;
2128 }
2129 
2130 //________________________________________________________________________
2132  fCutContainer->AddAt(cutsRP,0); fCutsRP=cutsRP; cutsRP->SetPOItype(0);
2133 }
2134 
2135 //________________________________________________________________________
2137  fCutContainer->AddAt(cutsPOI,1); fCutsPOI=cutsPOI; cutsPOI->SetPOItype(1);
2138 }
2139 
2140 //________________________________________________________________________
2142 {
2143  // Terminate analysis
2144  //
2145  /* if(fDebug > 1) printf(" **** AliAnalysisTaskCRCZDC::Terminate() \n");
2146 
2147  //fOutput = dynamic_cast<TList*> (GetOutputData(1));
2148  //if(!fOutput) printf("ERROR: fOutput not available\n");
2149  */
2150 }
void SetCharge(Int_t charge)
TH2F * fhZDCCvsZDCCA
ZNC vs ZNA;.
static const Int_t fkVZEROnQAplots
TH2F * fhZNvsVZERO
ZNC+ZNA vs ZPC+ZPA;.
TH1F * fPtSpecFB128[2][10]
PtSpecRec FB96.
void SetPOItype(Int_t poiType, Bool_t b=kTRUE)
void FindDaughters(Bool_t keepDaughtersInRPselection=kFALSE)
TH3D * fhZNSpectraCor
ZNA vs. centrality.
const Color_t cc[]
Definition: DrawKs.C:1
virtual void SetV02Qsub(Double_t QVCx, Double_t QVCy, Double_t MC, Double_t QVAx, Double_t QVAy, Double_t MA, Int_t har)
double Double_t
Definition: External.C:58
void SetEta(Double_t eta)
TH1F * fPileUpCount
centrality distribution
Int_t fRunList[fCRCMaxnRun]
Definition: External.C:236
static const Int_t fCRCnTow
TH1F * fhZNCPMQiPMC[4]
ZNA PM high res.
void SetZNCQ0(Double_t const en)
TH2F * fhZNCvsZPC
ZDCC vs ZDCCA.
static const Int_t fkVZEROnHar
virtual void SetZDC2Qsub(Double_t *QVC, Double_t MC, Double_t *QVA, Double_t MA)
void AddFlow(Double_t v1, Double_t v2, Double_t v3, Double_t v4, Double_t v5)
Int_t GetReferenceMultiplicity(AliVEvent *event, AliMCEvent *mcEvent)
centrality
void Set(const AliVParticle *p)
void SetZPAEnergy(Double_t const en)
TH2F * fhZNCvsAsymm
ZNA vs asymmetry.
void SetMCReactionPlaneAngle(const AliMCEvent *mcEvent)
TH3D * fhZNCenDis[2]
ZDC vs VZERO;.
void SetCutsPOI(AliFlowTrackCuts *cutsPOI)
void SetZNAQ0(Double_t const en)
TH1F * fPtSpecFB768[2][10]
PtSpecRec FB128.
TH1F * fPtSpecGen[2][10]
list with pt spectra
virtual void UserExec(Option_t *option)
void SetHistWeightvsPhiMax(Double_t d)
TCanvas * c
Definition: TestFitELoss.C:172
TH2F * fhZNCvscentrality
ZNC vs asymmetry.
void SetZNAEnergy(Double_t const en)
void AddTrack(AliFlowTrackSimple *track)
void SetReferenceMultiplicity(Int_t m)
virtual Int_t GetCenBin(Double_t Centrality)
void SetCutsRP(AliFlowTrackCuts *cutsRP)
TRandom * gRandom
TH3D * fTrackQADCAxy[fKNFBs][4]
TH2F * fhZNCvsZNA
PMQi/PMC for ZNA.
TH1F * fhAsymm
ZN centroid vs centrality.
void SetCentralityCL1(Double_t c)
virtual void SetVertexPosition(Double_t *pos)
void IncrementNumberOfPOIs(Int_t poiType=1)
TList * GetQA() const
TH3D * fhZNBCCorr
ZNA vs. centrality.
void SetForRPSelection(Bool_t b=kTRUE)
const Double_t etamin
void Fill(AliFlowTrackCuts *rpCuts, AliFlowTrackCuts *poiCuts)
TH2F * fhZNAvsZPA
ZNC vs ZPC;.
int Int_t
Definition: External.C:63
TH1F * fhZNAPM[5]
ZNC PM high res.
TProfile3D * fVZEROQVectorRecQxStored[fkVZEROnHar]
TF1 * fMultTOFLowCut
centrality distribution
unsigned int UInt_t
Definition: External.C:33
void SetShuffleTracks(Bool_t b)
float Float_t
Definition: External.C:68
AliMultSelection * fMultSelection
void SetNITSCL1(Double_t c)
Definition: External.C:252
Definition: External.C:228
virtual void UserCreateOutputObjects()
Definition: External.C:212
static AliFlowCommonConstants * GetMaster()
TH3D * fhZNSpectra
ZNA vs. centrality.
AliFlowEventCuts * fCutsEvent
const Double_t ptmin
void SetAbsOrbit(UInt_t const en)
TH1F * fPileUpMultSelCount
centrality distribution
Double_t GetCentrality() const
Bool_t fUseTowerEq
MultSelection (RUN2 centrality estimator)
Bool_t plpMV(const AliAODEvent *aod)
void SetCentrality(Double_t c)
AliFlowTrackCuts * fCutsPOI
void SetRun(Int_t const run)
void SetSource(trackSource s)
Definition: AliFlowTrack.h:37
TH2F * fhZPAvscentrality
ZNC vs. centrality.
void SetPhi(Double_t phi)
TProfile3D * fVZEROQVectorRecQyStored[fkVZEROnHar]
short Short_t
Definition: External.C:23
void DefineDeadZone(Double_t etaMin, Double_t etaMax, Double_t phiMin, Double_t phiMax)
AliGenEventHeader * fGenHeader
void SetPt(Double_t pt)
TProfile2D * fTrackQADphi[fKNFBs][4]
TProfile * fZNCTower[fCRCMaxnRun][fCRCnTow]
Q Vectors list per run.
TList * fCRCQVecListRun[fCRCMaxnRun]
Run list.
void SetZNCEnergy(Double_t const en)
virtual void Terminate(Option_t *option)
TProfile * fZNATower[fCRCMaxnRun][fCRCnTow]
ZNC tower spectra.
void SetEvent(AliVEvent *event, AliMCEvent *mcEvent=NULL)
TH2F * fhZPCvscentrality
ZNA vs. centrality.
Bool_t IsSetMCReactionPlaneAngle() const
TClonesArray * fStack
ZNA tower spectra.
TList * GetQA() const
TH2D * fTrackQApT[fKNFBs][4]
TProfile2D * fVZEROQVectorRecFinal[fkVZEROnHar][fkVZEROnQAplots]
static const Int_t fCRCMaxnRun
ZNA vs. centrality.
Double_t GetBadTowerResp(Double_t Et, TH2D *BadTowerCalibHist)
void SetForPOISelection(Bool_t b=kTRUE)
void SetHistWeightvsPhiMin(Double_t d)
TH1F * fhZNAPMQiPMC[4]
PMQi/PMC for ZNC.
void AddV2(Double_t v2)
TH3D * fTrackQADCAz[fKNFBs][4]
TH2F * fhZDCvsVZERO
ZN vs VZERO;.
void TagSubeventsInEta(Double_t etaMinA, Double_t etaMaxA, Double_t etaMinB, Double_t etaMaxB)
void SetZPCEnergy(Double_t const en)
unsigned short UShort_t
Definition: External.C:28
TH1F * fPtSpecFB96[2][10]
PtSpecRec FB32.
TH3D * fhZNSpectraPow
ZNA vs. centrality.
const char Option_t
Definition: External.C:48
TH2F * fhZNAvsAsymm
ZN asymmetry.
TH2F * fhZNAvscentrality
ZNC vs. centrality.
bool Bool_t
Definition: External.C:53
virtual Bool_t IsSelected(TObject *obj, TObject *objmc)
TH1F * fPtSpecFB32[2][10]
PtSpecGen.
TH1F * fhZNCPM[5]
list send on output slot 0
void SetCentralityTRK(Double_t c)
Double_t GetWDist(const AliVVertex *v0, const AliVVertex *v1)
Bool_t SelectPileup(AliAODEvent *aod)
Bool_t fCutTPC
PtSpecRec FB768.
AliFlowTrackCuts * fCutsRP
void InsertTrack(AliFlowTrack *)
virtual void ClearFast()
Bool_t fUseBadTowerCalib
list for storing calib files
AliAnalysisUtils * fAnalysisUtil
Event selection.
TH2F * fhZNvsZP
ZNA vs ZPA;.
AliGenHijingEventHeader * fHijingGenHeader
const Double_t phimin
AliGenPythiaEventHeader * fPythiaGenHeader
Int_t NumberOfTracks() const