AliPhysics  48852ec (48852ec)
 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 
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 fhZDCvsTracklets(0x0),
177 fhZDCvsNclu1(0x0),
178 fhDebunch(0x0),
179 fhAsymm(0x0),
180 fhZNAvsAsymm(0x0),
181 fhZNCvsAsymm(0x0),
182 fhZNCvscentrality(0x0),
183 fhZNAvscentrality(0x0),
184 fCRCnRun(0),
185 fZDCGainAlpha(0.395),
186 fDataSet(kAny),
187 fStack(0x0),
188 fCutTPC(kFALSE),
189 fCenDis(0x0),
190 fVZEROMult(0x0),
191 fMultSelection(0x0),
192 fPileUpCount(0x0),
193 fPileUpMultSelCount(0x0),
194 fMultTOFLowCut(0x0),
195 fMultTOFHighCut(0x0),
196 fUseTowerEq(kFALSE),
197 fTowerEqList(NULL),
198 fUseBadTowerCalib(kFALSE),
199 fBadTowerCalibList(NULL),
200 fVZEROGainEqList(NULL),
201 fVZEROQVecRecList(NULL),
202 fUseZDCSpectraCorr(kFALSE),
203 fCorrectPhiTracklets(kFALSE),
204 fZDCSpectraCorrList(NULL),
205 fSpectraMCList(NULL),
206 fTrackQAList(NULL),
207 fBadTowerStuffList(NULL),
208 fVZEROStuffList(NULL),
209 fCachedRunNum(0),
210 fhZNSpectra(0x0),
211 fhZNSpectraCor(0x0),
212 fhZNSpectraPow(0x0),
213 fhZNBCCorr(0x0)
214 {
215  for(int i=0; i<5; i++){
216  fhZNCPM[i] = 0x0;
217  fhZNAPM[i] = 0x0;
218  }
219  for(int i=0; i<4; i++){
220  fhZNCPMQiPMC[i] = 0x0;
221  fhZNAPMQiPMC[i] = 0x0;
222  }
223  for(Int_t r=0; r<fCRCMaxnRun; r++) {
224  fRunList[r] = 0;
225  }
226  for(Int_t c=0; c<2; c++) {
227  for(Int_t i=0; i<5; i++) {
228  fTowerGainEq[c][i] = NULL;
229  }
230  }
231  for(Int_t c=0; c<100; c++) {
232  fBadTowerCalibHist[c] = NULL;
233  }
234  fVZEROGainEqHist = NULL;
235  for (Int_t k=0; k<fkVZEROnHar; k++) {
236 // fVZEROQVectorRecQx[k] = NULL;
237 // fVZEROQVectorRecQy[k] = NULL;
238  fVZEROQVectorRecQxStored[k] = NULL;
239  fVZEROQVectorRecQyStored[k] = NULL;
240  for (Int_t t=0; t<fkVZEROnQAplots; t++) {
241  fVZEROQVectorRecFinal[k][t] = NULL;
242  }
243  }
244  for(Int_t i=0; i<8; i++) {
245  SpecCorMu1[i] = NULL;
246  SpecCorMu2[i] = NULL;
247  SpecCorSi[i] = NULL;
248  SpecCorAv[i] = NULL;
249  }
250  this->InitializeRunArrays();
251  fMyTRandom3 = new TRandom3(1);
252  gRandom->SetSeed(fMyTRandom3->Integer(65539));
253  for(Int_t j=0; j<2; j++) {
254  for(Int_t c=0; c<10; c++) {
255  fPtSpecGen[j][c] = NULL;
256  fPtSpecFB32[j][c] = NULL;
257  fPtSpecFB96[j][c] = NULL;
258  fPtSpecFB128[j][c] = NULL;
259  fPtSpecFB768[j][c] = NULL;
260  }
261  }
262  for (Int_t c=0; c<2; c++) {
263  fhZNCenDis[c] = NULL;
264  }
265  fMinRingVZC=1;
266  fMaxRingVZC=4;
267  fMinRingVZA=5;
268  fMaxRingVZA=8;
269  for(Int_t fb=0; fb<fKNFBs; fb++){
270  for (Int_t i=0; i<4; i++) {
271  fTrackQADCAxy[fb][i] = NULL;
272  fTrackQADCAz[fb][i] = NULL;
273  fTrackQApT[fb][i] = NULL;
274  fTrackQADphi[fb][i] = NULL;
275  fEbEQRe[fb][i] = NULL;
276  fEbEQIm[fb][i] = NULL;
277  fEbEQMu[fb][i] = NULL;
278  }
279  }
280 }
281 
282 //________________________________________________________________________
283 AliAnalysisTaskCRCZDC::AliAnalysisTaskCRCZDC(const char *name, TString RPtype, Bool_t on, UInt_t iseed, Bool_t bCandidates):
284 AliAnalysisTaskSE(name),
285 fAnalysisType(kAUTOMATIC),
286 fRPType(RPtype),
287 fCFManager1(NULL),
288 fCFManager2(NULL),
289 fCutsEvent(NULL),
290 fCutsRP(NULL),
291 fCutsPOI(NULL),
292 fCutContainer(new TList()),
293 fQAList(NULL),
294 fMinMult(0),
295 fMaxMult(10000000),
296 fMinA(-1.0),
297 fMaxA(-0.01),
298 fMinB(0.01),
299 fMaxB(1.0),
300 fQAon(on),
301 fLoadCandidates(bCandidates),
302 fNbinsMult(10000),
303 fNbinsPt(100),
304 fNbinsPhi(100),
305 fNbinsEta(200),
306 fNbinsQ(500),
307 fNbinsMass(1),
308 fMultMin(0.),
309 fMultMax(10000.),
310 fPtMin(0.),
311 fPtMax(10.),
312 fPhiMin(0.),
313 fPhiMax(TMath::TwoPi()),
314 fEtaMin(-5.),
315 fEtaMax(5.),
316 fQMin(0.),
317 fQMax(3.),
318 fMassMin(-1.),
319 fMassMax(0.),
320 fHistWeightvsPhiMin(0.),
321 fHistWeightvsPhiMax(3.),
322 fExcludedEtaMin(0.),
323 fExcludedEtaMax(0.),
324 fExcludedPhiMin(0.),
325 fExcludedPhiMax(0.),
326 fAfterburnerOn(kFALSE),
327 fNonFlowNumberOfTrackClones(0),
328 fV1(0.),
329 fV2(0.),
330 fV3(0.),
331 fV4(0.),
332 fV5(0.),
333 fDifferentialV2(0),
334 fFlowEvent(NULL),
335 fShuffleTracks(kFALSE),
336 fMyTRandom3(NULL),
337 fAnalysisInput(kAOD),
338 fIsMCInput(kFALSE),
339 fUseMCCen(kTRUE),
340 fRejectPileUp(kTRUE),
341 fRejectPileUpTight(kFALSE),
342 fResetNegativeZDC(kFALSE),
343 fCentrLowLim(0.),
344 fCentrUpLim(100.),
345 fCentrEstimator(kV0M),
346 fOutput(0x0),
347 fhZNCvsZNA(0x0),
348 fhZDCCvsZDCCA(0x0),
349 fhZNCvsZPC(0x0),
350 fhZNAvsZPA(0x0),
351 fhZNvsZP(0x0),
352 fhZNvsVZERO(0x0),
353 fhZDCvsVZERO(0x0),
354 fhZDCvsTracklets(0x0),
355 fhZDCvsNclu1(0x0),
356 fhDebunch(0x0),
357 fhAsymm(0x0),
358 fhZNAvsAsymm(0x0),
359 fhZNCvsAsymm(0x0),
360 fhZNCvscentrality(0x0),
361 fhZNAvscentrality(0x0),
362 fDataSet(kAny),
363 fCRCnRun(0),
364 fZDCGainAlpha(0.395),
365 fGenHeader(NULL),
366 fPythiaGenHeader(NULL),
367 fHijingGenHeader(NULL),
368 fFlowTrack(NULL),
369 fAnalysisUtil(NULL),
370 fStack(0x0),
371 fCutTPC(kFALSE),
372 fCenDis(0x0),
373 fVZEROMult(0x0),
374 fMultSelection(0x0),
375 fPileUpCount(0x0),
376 fPileUpMultSelCount(0x0),
377 fMultTOFLowCut(0x0),
378 fMultTOFHighCut(0x0),
379 fUseTowerEq(kFALSE),
380 fTowerEqList(NULL),
381 fUseBadTowerCalib(kFALSE),
382 fBadTowerCalibList(NULL),
383 fVZEROGainEqList(NULL),
384 fVZEROQVecRecList(NULL),
385 fUseZDCSpectraCorr(kFALSE),
386 fCorrectPhiTracklets(kFALSE),
387 fZDCSpectraCorrList(NULL),
388 fSpectraMCList(NULL),
389 fTrackQAList(NULL),
390 fBadTowerStuffList(NULL),
391 fVZEROStuffList(NULL),
392 fCachedRunNum(0),
393 fhZNSpectra(0x0),
394 fhZNSpectraCor(0x0),
395 fhZNSpectraPow(0x0),
396 fhZNBCCorr(0x0)
397 {
398 
399  for(int i=0; i<5; i++){
400  fhZNCPM[i] = 0x0;
401  fhZNAPM[i] = 0x0;
402  }
403  for(int i=0; i<4; i++){
404  fhZNCPMQiPMC[i] = 0x0;
405  fhZNAPMQiPMC[i] = 0x0;
406  }
407  for(Int_t r=0; r<fCRCMaxnRun; r++) {
408  fRunList[r] = 0;
409  }
410  for(Int_t c=0; c<2; c++) {
411  for(Int_t i=0; i<5; i++) {
412  fTowerGainEq[c][i] = NULL;
413  }
414  }
415  for(Int_t c=0; c<100; c++) {
416  fBadTowerCalibHist[c] = NULL;
417  }
418  fVZEROGainEqHist = NULL;
419  for (Int_t k=0; k<fkVZEROnHar; k++) {
420  // fVZEROQVectorRecQx[k] = NULL;
421  // fVZEROQVectorRecQy[k] = NULL;
422  fVZEROQVectorRecQxStored[k] = NULL;
423  fVZEROQVectorRecQyStored[k] = NULL;
424  for (Int_t t=0; t<fkVZEROnQAplots; t++) {
425  fVZEROQVectorRecFinal[k][t] = NULL;
426  }
427  }
428  for(Int_t i=0; i<8; i++) {
429  SpecCorMu1[i] = NULL;
430  SpecCorMu2[i] = NULL;
431  SpecCorSi[i] = NULL;
432  SpecCorAv[i] = NULL;
433  }
434  this->InitializeRunArrays();
435  fMyTRandom3 = new TRandom3(iseed);
436  gRandom->SetSeed(fMyTRandom3->Integer(65539));
437 
438  DefineInput(0, TChain::Class());
439  // Define output slots here
440  // Define here the flow event output
441  DefineOutput(1, AliFlowEventSimple::Class());
442  DefineOutput(2, TList::Class());
443 
444  for(Int_t j=0; j<2; j++) {
445  for(Int_t c=0; c<10; c++) {
446  fPtSpecGen[j][c] = NULL;
447  fPtSpecFB32[j][c] = NULL;
448  fPtSpecFB96[j][c] = NULL;
449  fPtSpecFB128[j][c] = NULL;
450  fPtSpecFB768[j][c] = NULL;
451  }
452  }
453  for (Int_t c=0; c<2; c++) {
454  fhZNCenDis[c] = NULL;
455  }
456  fMinRingVZC=1;
457  fMaxRingVZC=4;
458  fMinRingVZA=5;
459  fMaxRingVZA=8;
460  for(Int_t fb=0; fb<fKNFBs; fb++){
461  for (Int_t i=0; i<4; i++) {
462  fTrackQADCAxy[fb][i] = NULL;
463  fTrackQADCAz[fb][i] = NULL;
464  fTrackQApT[fb][i] = NULL;
465  fTrackQADphi[fb][i] = NULL;
466  fEbEQRe[fb][i] = NULL;
467  fEbEQIm[fb][i] = NULL;
468  fEbEQMu[fb][i] = NULL;
469  }
470  }
471 }
472 
473 //________________________________________________________________________
475 {
476  // Destructor
477  if(fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()){
478  delete fOutput; fOutput=0;
479  }
480  delete fMyTRandom3;
481  delete fFlowEvent;
482  delete fFlowTrack;
483  delete fCutsEvent;
484  if (fTowerEqList) delete fTowerEqList;
489  if (fAnalysisUtil) delete fAnalysisUtil;
490  if (fQAList) delete fQAList;
491  if (fCutContainer) fCutContainer->Delete(); delete fCutContainer;
492 }
493 
494 //________________________________________________________________________
496 {
497  for(Int_t r=0;r<fCRCMaxnRun;r++) {
498  fCRCQVecListRun[r] = NULL;
499  for(Int_t k=0;k<fCRCnTow;k++) {
500  fZNCTower[r][k] = NULL;
501  fZNATower[r][k] = NULL;
502  }
503 // fhZNSpectraRbR[r] = NULL;
504  }
505  // for(Int_t i=0;i<fnCen;i++) {
506  // fPtPhiEtaRbRFB128[r][i] = NULL;
507  // fPtPhiEtaRbRFB768[r][i] = NULL;
508  // }
509  // }
510 }
511 
512 //________________________________________________________________________
514 {
515  // Create the output containers
516 
517  //set the common constants
520  cc->SetNbinsPt(fNbinsPt);
521  cc->SetNbinsPhi(fNbinsPhi);
522  cc->SetNbinsEta(fNbinsEta);
523  cc->SetNbinsQ(fNbinsQ);
525  cc->SetMultMin(fMultMin);
526  cc->SetMultMax(fMultMax);
527  cc->SetPtMin(fPtMin);
528  cc->SetPtMax(fPtMax);
529  cc->SetPhiMin(fPhiMin);
530  cc->SetPhiMax(fPhiMax);
531  cc->SetEtaMin(fEtaMin);
532  cc->SetEtaMax(fEtaMax);
533  cc->SetQMin(fQMin);
534  cc->SetQMax(fQMax);
535  cc->SetMassMin(fMassMin);
536  cc->SetMassMax(fMassMax);
539 
540  fFlowEvent = new AliFlowEvent(20000);
541  fFlowTrack = new AliFlowTrack();
542 
543  //printf(" AliAnalysisTaskCRCZDC::UserCreateOutputObjects()\n\n");
544  fOutput = new TList();
545  fOutput->SetOwner(kTRUE);
546  //fOutput->SetName("output");
547 
548  if (fQAon) {
549  fQAList = new TList();
550  fQAList->SetOwner(kTRUE);
551  fQAList->SetName("AliFlowEventCuts QA");
552  if (fCutsEvent->GetQA()) fQAList->Add(fCutsEvent->GetQA()); //0
553  if (fCutsRP->GetQA()) fQAList->Add(fCutsRP->GetQA()); //1
554  if (fCutsPOI->GetQA())fQAList->Add(fCutsPOI->GetQA()); //2
555  fOutput->Add(fQAList);
556  }
557 
558  fVZEROStuffList = new TList();
559  fVZEROStuffList->SetOwner(kTRUE);
560  fVZEROStuffList->SetName("VZERO stuff");
561  fOutput->Add(fVZEROStuffList);
562 
563  fBadTowerStuffList = new TList();
564  fBadTowerStuffList->SetOwner(kTRUE);
565  fBadTowerStuffList->SetName("BadTowerCalib");
567 
568  fCenDis = new TH1F("fCenDis", "fCenDis", 100, 0., 100.);
569  fOutput->Add(fCenDis);
570  fPileUpCount = new TH1F("fPileUpCount", "fPileUpCount", 9, 0., 9.);
571  fPileUpCount->GetXaxis()->SetBinLabel(1,"plpMV");
572  fPileUpCount->GetXaxis()->SetBinLabel(2,"fromSPD");
573  fPileUpCount->GetXaxis()->SetBinLabel(3,"RefMultiplicityComb08");
574  fPileUpCount->GetXaxis()->SetBinLabel(4,"IncompleteDAQ");
575  fPileUpCount->GetXaxis()->SetBinLabel(5,"abs(V0M-CL1)>7.5");
576  fPileUpCount->GetXaxis()->SetBinLabel(6,"missingVtx");
577  fPileUpCount->GetXaxis()->SetBinLabel(7,"inconsistentVtx");
578  fPileUpCount->GetXaxis()->SetBinLabel(8,"multESDTPCDif");
579  fPileUpCount->GetXaxis()->SetBinLabel(9,"extraPileUpMultSel");
580  fOutput->Add(fPileUpCount);
581  fPileUpMultSelCount = new TH1F("fPileUpMultSelCount", "fPileUpMultSelCount", 8, 0., 8.);
582  fPileUpMultSelCount->GetXaxis()->SetBinLabel(1,"IsNotPileup");
583  fPileUpMultSelCount->GetXaxis()->SetBinLabel(2,"IsNotPileupMV");
584  fPileUpMultSelCount->GetXaxis()->SetBinLabel(3,"IsNotPileupInMultBins");
585  fPileUpMultSelCount->GetXaxis()->SetBinLabel(4,"InconsistentVertices");
586  fPileUpMultSelCount->GetXaxis()->SetBinLabel(5,"TrackletVsCluster");
587  fPileUpMultSelCount->GetXaxis()->SetBinLabel(6,"AsymmetricInVZERO");
588  fPileUpMultSelCount->GetXaxis()->SetBinLabel(7,"IncompleteDAQ");
589  fPileUpMultSelCount->GetXaxis()->SetBinLabel(8,"GoodVertex2016");
591 
592  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);
593  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);
594  fOutput->Add(fMultTOFLowCut);
595  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);
596  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);
597  fOutput->Add(fMultTOFHighCut);
598 
599  for(Int_t c=0; c<2; c++) {
600  for(Int_t i=0; i<5; i++) {
601  fTowerGainEq[c][i] = new TH1D();
602  fOutput->Add(fTowerGainEq[c][i]);
603  }
604  }
605 
606  for (Int_t c=0; c<2; c++) {
607  fhZNCenDis[c] = new TH3D(Form("fhZNCenDis[%d]",c), Form("fhZNCenDis[%d]",c), 100, 0., 100., 100, -2., 2. , 100., -2., 2.);
608  fOutput->Add(fhZNCenDis[c]);
609  }
610 
611  if(fBadTowerCalibList) {
612  for(Int_t c=0; c<100; c++) {
613  fBadTowerCalibHist[c] = (TH2D*)fBadTowerCalibList->FindObject(Form("TH2Resp[%d]",c));
615  }
616  }
617  if(fZDCSpectraCorrList) {
618  for(Int_t i=0; i<8; i++) {
619  SpecCorMu1[i] = (TH1D*)fZDCSpectraCorrList->FindObject(Form("SpecCorMu1[%d]",i));
620  fOutput->Add(SpecCorMu1[i]);
621  SpecCorMu2[i] = (TH1D*)fZDCSpectraCorrList->FindObject(Form("SpecCorMu2[%d]",i));
622  fOutput->Add(SpecCorMu2[i]);
623  SpecCorAv[i] = (TH1D*)fZDCSpectraCorrList->FindObject(Form("SpecCorAv[%d]",i));
624  fOutput->Add(SpecCorAv[i]);
625  SpecCorSi[i] = (TH1D*)fZDCSpectraCorrList->FindObject(Form("SpecCorSi[%d]",i));
626  fOutput->Add(SpecCorSi[i]);
627  }
628  }
629 
630  fhZNSpectra = new TH3D("fhZNSpectra","fhZNSpectra",100,0.,100.,8,0.,8.,1000,0.,1.E5);
631  fOutput->Add(fhZNSpectra);
632  fhZNSpectraCor = new TH3D("fhZNSpectraCor","fhZNSpectraCor",100,0.,100.,8,0.,8.,1000,0.,1.E5);
633  fOutput->Add(fhZNSpectraCor);
634  fhZNSpectraPow = new TH3D("fhZNSpectraPow","fhZNSpectraPow",100,0.,100.,8,0.,8.,1000,0.,TMath::Power(1.E5,fZDCGainAlpha));
635  fOutput->Add(fhZNSpectraPow);
636  fhZNBCCorr = new TH3D("fhZNBCCorr","fhZNBCCorr",100,0.,100.,500,0.,1.E5,500,0.,1.E5);
637  fOutput->Add(fhZNBCCorr);
638 
639  if(fAnalysisType == kMCAOD) {
640 
641  fSpectraMCList = new TList();
642  fSpectraMCList->SetOwner(kTRUE);
643  fSpectraMCList->SetName("Spectra");
644  fOutput->Add(fSpectraMCList);
645 
646  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.};
647  for(Int_t j=0; j<2; j++) {
648  for(Int_t c=0; c<10; c++) {
649  fPtSpecGen[j][c] = new TH1F(Form("fPtSpecGen[%d][%d]",j,c), Form("fPtSpecGen[%d][%d]",j,c), 23, xmin);
650  fSpectraMCList->Add(fPtSpecGen[j][c]);
651  fPtSpecFB32[j][c] = new TH1F(Form("fPtSpecFB32[%d][%d]",j,c), Form("fPtSpecFB32[%d][%d]",j,c), 23, xmin);
652  fSpectraMCList->Add(fPtSpecFB32[j][c]);
653  fPtSpecFB96[j][c] = new TH1F(Form("fPtSpecFB96[%d][%d]",j,c), Form("fPtSpecFB96[%d][%d]",j,c), 23, xmin);
654  fSpectraMCList->Add(fPtSpecFB96[j][c]);
655  fPtSpecFB128[j][c] = new TH1F(Form("fPtSpecFB128[%d][%d]",j,c), Form("fPtSpecFB128[%d][%d]",j,c), 23, xmin);
656  fSpectraMCList->Add(fPtSpecFB128[j][c]);
657  fPtSpecFB768[j][c] = new TH1F(Form("fPtSpecFB768[%d][%d]",j,c), Form("fPtSpecFB768[%d][%d]",j,c), 23, xmin);
658  fSpectraMCList->Add(fPtSpecFB768[j][c]);
659  }
660  }
661  }
662 
663  fAnalysisUtil = new AliAnalysisUtils();
664  fAnalysisUtil->SetUseMVPlpSelection(kTRUE);
665  fAnalysisUtil->SetUseOutOfBunchPileUp(kTRUE);
666 
667  for(int i=0; i<5; i++){
668  char hname[20];
669  sprintf(hname,"hZNCPM%d",i);
670  fhZNCPM[i] = new TH1F(hname, hname, 200, -50., 140000);
671  fOutput->Add(fhZNCPM[i]);
672  //
673  sprintf(hname,"hZNAPM%d",i);
674  fhZNAPM[i] = new TH1F(hname, hname, 200, -50., 140000);
675  fOutput->Add(fhZNAPM[i]);
676  //
677  if(i<4){
678  //
679  char hnamenc[20];
680  sprintf(hnamenc, "hZNCPMQ%dPMC",i+1);
681  fhZNCPMQiPMC[i] = new TH1F(hnamenc, hnamenc, 100, 0., 1.);
682  fOutput->Add(fhZNCPMQiPMC[i]);
683  //
684  char hnamena[20];
685  sprintf(hnamena, "hZNAPMQ%dPMC",i+1);
686  fhZNAPMQiPMC[i] = new TH1F(hnamena, hnamena, 100, 0., 1.);
687  fOutput->Add(fhZNAPMQiPMC[i]);
688  }
689  }
690 
691  fhZNCvsZNA = new TH2F("hZNCvsZNA","hZNCvsZNA",200,-50.,140000,200,-50.,140000);
692  fOutput->Add(fhZNCvsZNA);
693  fhZDCCvsZDCCA = new TH2F("hZDCCvsZDCCA","hZDCCvsZDCCA",200,0.,180000.,200,0.,200000.);
694  fOutput->Add(fhZDCCvsZDCCA);
695  fhZNCvsZPC = new TH2F("hZNCvsZPC","hZNCvsZPC",200,-50.,50000,200,-50.,140000);
696  fOutput->Add(fhZNCvsZPC);
697  fhZNAvsZPA = new TH2F("hZNAvsZPA","hZNAvsZPA",200,-50.,50000,200,-50.,140000);
698  fOutput->Add(fhZNAvsZPA);
699  fhZNvsZP = new TH2F("hZNvsZP","hZNvsZP",200,-50.,80000,200,-50.,200000);
700  fOutput->Add(fhZNvsZP);
701  fhZNvsVZERO = new TH2F("hZNvsVZERO","hZNvsVZERO",250,0.,25000.,200,0.,200000.);
702  fOutput->Add(fhZNvsVZERO);
703  fhZDCvsVZERO = new TH2F("hZDCvsVZERO","hZDCvsVZERO",250,0.,25000.,250,0.,250000.);
704  fOutput->Add(fhZDCvsVZERO);
705  fhZDCvsTracklets = new TH2F("hZDCvsTracklets","hZDCvsTracklets",200,0.,4000.,250,0.,250000.);
707  fhZDCvsNclu1 = new TH2F("hZDCvsNclu1", "hZDCvsNclu1", 200, 0.,8000.,200,0.,250000.);
708  fOutput->Add(fhZDCvsNclu1);
709  fhDebunch = new TH2F("hDebunch","hDebunch",240,-100.,-40.,240,-30.,30.);
710  fOutput->Add(fhDebunch);
711 
712  fhAsymm = new TH1F("hAsymm" , "Asimmetry ",200,-1.,1.);
713  fOutput->Add(fhAsymm);
714  fhZNAvsAsymm = new TH2F("hZNAvsAsymm","ZNA vs. asymm.",200,-1.,1.,200,0.,80.);
715  fOutput->Add(fhZNAvsAsymm);
716  fhZNCvsAsymm = new TH2F("hZNCvsAsymm","ZNC vs. asymm.",200,-1.,1.,200,0.,80.);
717  fOutput->Add(fhZNCvsAsymm);
718 
719  fhZNCvscentrality = new TH2F("hZNCvscentrality","ZNC vs. centrality",100,0.,100.,100,0.,100.);
721  fhZNAvscentrality = new TH2F("hZNAvscentrality","ZNA vs. centrality",100,0.,100.,100,0.,100.);
723 
724  //********************************************************************
725 
726  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};
727 
728  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};
729 
730  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};
731 
732  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};
733 
734  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};
735 
736  if(fDataSet==k2010) {fCRCnRun=92;}
737  if(fDataSet==k2011) {fCRCnRun=119;}
738  if(fDataSet==k2015) {fCRCnRun=90;}
739  if(fDataSet==k2015v6) {fCRCnRun=91;}
740  if(fDataSet==k2015pidfix) {fCRCnRun=35;}
741  if(fDataSet==kAny) {fCRCnRun=1;}
742 
743  Int_t d=0;
744  for(Int_t r=0; r<fCRCnRun; r++) {
745  if(fDataSet==k2010) fRunList[d] = dRun10h[r];
746  if(fDataSet==k2011) fRunList[d] = dRun11h[r];
747  if(fDataSet==k2015) fRunList[d] = dRun15o[r];
748  if(fDataSet==k2015v6) fRunList[d] = dRun15ov6[r];
749  if(fDataSet==k2015pidfix) fRunList[d] = dRun15opidfix[r];
750  if(fDataSet==kAny) fRunList[d] = 1;
751  d++;
752  }
753 
754  fVZEROMult = new TProfile2D("fVZEROMult","fVZEROMult",fCRCnRun,0.,1.*fCRCnRun,64,0.,64.);
755  for (Int_t i=0; i<fCRCnRun; i++) {
756  fVZEROMult->GetXaxis()->SetBinLabel(i+1,Form("%d",fRunList[i]));
757  }
759 
760  if(fVZEROGainEqList) {
761  fVZEROGainEqHist = (TH2D*)fVZEROGainEqList->FindObject("VZEROEqGain");
763  }
764  if(fVZEROQVecRecList) {
765  for (Int_t k=0; k<fkVZEROnHar; k++) {
766  fVZEROQVectorRecQxStored[k] = (TProfile3D*)fVZEROQVecRecList->FindObject(Form("fVZEROQVectorRecQx[%d]",k));
767  fVZEROQVectorRecQxStored[k]->SetTitle(Form("fVZEROQVectorRecQxStored[%d]",k));
768  fVZEROQVectorRecQxStored[k]->SetName(Form("fVZEROQVectorRecQxStored[%d]",k));
770  fVZEROQVectorRecQyStored[k] = (TProfile3D*)fVZEROQVecRecList->FindObject(Form("fVZEROQVectorRecQy[%d]",k));
771  fVZEROQVectorRecQyStored[k]->SetTitle(Form("fVZEROQVectorRecQyStored[%d]",k));
772  fVZEROQVectorRecQyStored[k]->SetName(Form("fVZEROQVectorRecQyStored[%d]",k));
774  for (Int_t t=0; t<fkVZEROnQAplots; t++) {
775  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");
776  fVZEROQVectorRecFinal[k][t]->Sumw2();
778  }
779  }
780  }
781 
782 // for (Int_t k=0; k<fkVZEROnHar; k++) {
783 // fVZEROQVectorRecQx[k] = new TProfile3D(Form("fVZEROQVectorRecQx[%d]",k),Form("fVZEROQVectorRecQx[%d]",k),fCRCnRun,0.,1.*fCRCnRun,100,0.,100.,8,0.,8.,"s");
784 // fVZEROQVectorRecQx[k]->Sumw2();
785 // fVZEROStuffList->Add(fVZEROQVectorRecQx[k]);
786 // fVZEROQVectorRecQy[k] = new TProfile3D(Form("fVZEROQVectorRecQy[%d]",k),Form("fVZEROQVectorRecQy[%d]",k),fCRCnRun,0.,1.*fCRCnRun,100,0.,100.,8,0.,8.,"s");
787 // fVZEROQVectorRecQy[k]->Sumw2();
788 // fVZEROStuffList->Add(fVZEROQVectorRecQy[k]);
789 // }
790 
791  // track QA
792 
793  if(fAnalysisType == kTrackQA) {
794 
795  fTrackQAList = new TList();
796  fTrackQAList->SetOwner(kTRUE);
797  fTrackQAList->SetName("TrackQA");
798  fOutput->Add(fTrackQAList);
799 
800  Int_t nBinsEta=10;
801  Double_t dBinsEta[]={-0.8,-0.64,-0.48,-0.32,-0.16,0.,0.16,0.32,0.48,0.64,0.8};
802  Int_t nBinsPt = 26;
803  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.};
804 
805  for(Int_t fb=0; fb<fKNFBs; fb++){
806  for (Int_t i=0; i<4; i++) {
807  // DCA
808  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);
809  fTrackQAList->Add(fTrackQADCAxy[fb][i]);
810  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);
811  fTrackQAList->Add(fTrackQADCAz[fb][i]);
812  // pT
813  fTrackQApT[fb][i] = new TH2D(Form("fTrackQApT[%d][%d]",fb,i),";#eta;p_{T} [GeV/c]",nBinsEta,dBinsEta,nBinsPt,dBinsPt);
814  fTrackQAList->Add(fTrackQApT[fb][i]);
815  // Dphi
816  fTrackQADphi[fb][i] = new TProfile2D(Form("fTrackQADphi[%d][%d]",fb,i),";#eta;p_{T} [GeV/c];cos(#Delta#phi)",nBinsEta,dBinsEta,nBinsPt,dBinsPt);
817  fTrackQAList->Add(fTrackQADphi[fb][i]);
818  // EbE
819  fEbEQRe[fb][i] = new TH2D(Form("fEbEQRe[%d][%d]",fb,i),";#eta;p_{T} [GeV/c];QRe(EbE)",nBinsEta,dBinsEta,nBinsPt,dBinsPt);
820  fTrackQAList->Add(fEbEQRe[fb][i]);
821  fEbEQIm[fb][i] = new TH2D(Form("fEbEQIm[%d][%d]",fb,i),";#eta;p_{T} [GeV/c];QRe(EbE)",nBinsEta,dBinsEta,nBinsPt,dBinsPt);
822  fTrackQAList->Add(fEbEQIm[fb][i]);
823  fEbEQMu[fb][i] = new TH2D(Form("fEbEQMu[%d][%d]",fb,i),";#eta;p_{T} [GeV/c];QRe(EbE)",nBinsEta,dBinsEta,nBinsPt,dBinsPt);
824  fTrackQAList->Add(fEbEQMu[fb][i]);
825  }
826  }
827 
828  }
829 
830  // run-by-run stuff
831 
832  if(!fUseTowerEq) {
833  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.};
834  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.};
835  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};
836 
837  for(Int_t r=0;r<fCRCnRun;r++) {
838  fCRCQVecListRun[r] = new TList();
839  fCRCQVecListRun[r]->SetName(Form("Run %d",fRunList[r]));
840  fCRCQVecListRun[r]->SetOwner(kTRUE);
841  fOutput->Add(fCRCQVecListRun[r]);
842 
843  for(Int_t k=0;k<fCRCnTow;k++) {
844  fZNCTower[r][k] = new TProfile(Form("fZNCTower[%d][%d]",fRunList[r],k),Form("fZNCTower[%d][%d]",fRunList[r],k),100,0.,100.,"s");
845  fZNCTower[r][k]->Sumw2();
846  fCRCQVecListRun[r]->Add(fZNCTower[r][k]);
847  fZNATower[r][k] = new TProfile(Form("fZNATower[%d][%d]",fRunList[r],k),Form("fZNATower[%d][%d]",fRunList[r],k),100,0.,100.,"s");
848  fZNATower[r][k]->Sumw2();
849  fCRCQVecListRun[r]->Add(fZNATower[r][k]);
850  }
851 
852  // fhZNSpectraRbR[r] = new TH3D(Form("fhZNSpectraRbR[%d]",fRunList[r]),Form("fhZNSpectraRbR[%d]",fRunList[r]),50,0.,100.,8,0.,8.,100,0.,1.E5);
853  // fCRCQVecListRun[r]->Add(fhZNSpectraRbR[r]);
854 
855  // for(Int_t i=0;i<fnCen;i++) {
856  // fPtPhiEtaRbRFB128[r][i] = new TH3F(Form("fPtPhiEtaRbRFB128[%d][%d]",r,i),Form("fPtPhiEtaRbRFB128[%d][%d]",r,i),14, ptmin, 16, phimin, 16, etamin);
857  // fCRCQVecListRun[r]->Add(fPtPhiEtaRbRFB128[r][i]);
858  // fPtPhiEtaRbRFB768[r][i] = new TH3F(Form("fPtPhiEtaRbRFB768[%d][%d]",r,i),Form("fPtPhiEtaRbRFB768[%d][%d]",r,i),14, ptmin, 16, phimin, 16, etamin);
859  // fCRCQVecListRun[r]->Add(fPtPhiEtaRbRFB768[r][i]);
860  // }
861  }
862  }
863 
864  PostData(2, fOutput);
865 }
866 
867 //________________________________________________________________________
869 {
870  // Execute analysis for current event:
871  AliMCEvent* McEvent = MCEvent();
872  AliAODEvent *aod = dynamic_cast<AliAODEvent*>(InputEvent());
873  AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
874  // AliMultiplicity* myTracklets = NULL;
875  // AliESDPmdTrack* pmdtracks = NULL;
876  // int availableINslot=1;
877 
878  if (!(fCutsRP&&fCutsPOI&&fCutsEvent)) {
879  AliError("cuts not set");
880  return;
881  }
882 
883  Int_t RunBin=-1, bin=0;
884  Int_t RunNum = aod->GetRunNumber();
885  for(Int_t c=0;c<fCRCnRun;c++) {
886  if(fRunList[c]==RunNum) RunBin=bin;
887  else bin++;
888  }
889  if(RunBin==-1) return;
890  if(fDataSet==kAny) RunBin=0;
891 
892  //DEFAULT - automatically takes care of everything
894 
895  // get centrality
896  Double_t centrV0M=300, centrCL1=300, centrCL0=300, centrTRK=300;
898  centrV0M = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("V0M");
899  centrCL1 = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("CL1");
900  centrCL0 = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("CL0");
901  centrTRK = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("TRK");
902  } else {
903  fMultSelection = (AliMultSelection*) InputEvent()->FindListObject("MultSelection");
904  if(!fMultSelection) {
905  //If you get this warning (and lPercentiles 300) please check that the AliMultSelectionTask actually ran (before your task)
906  AliWarning("AliMultSelection object not found!");
907  } else {
908  centrV0M = fMultSelection->GetMultiplicityPercentile("V0M");
909  centrCL1 = fMultSelection->GetMultiplicityPercentile("CL1");
910  centrCL0 = fMultSelection->GetMultiplicityPercentile("CL0");
911  centrTRK = fMultSelection->GetMultiplicityPercentile("TRK");
912  }
913  }
914 
915  //check event cuts
916  if (InputEvent()) {
917  if(!fCutsEvent->IsSelected(InputEvent(),MCEvent())) return;
918  if(fRejectPileUp) {
919  Bool_t IsPileUp = SelectPileup(aod);
920  if(IsPileUp) return;
921  }
922  }
923 
924  //first attach all possible information to the cuts
925  fCutsRP->SetEvent( InputEvent(), MCEvent() ); //attach event
926  fCutsPOI->SetEvent( InputEvent(), MCEvent() );
927 
928  //then make the event
930 
931  if(fAnalysisType == kTracklets) {
932  // fill with tracklets
933  AliAODTracklets* anInputTracklets = (AliAODTracklets*)aod->GetTracklets();
934  Int_t multSPD = anInputTracklets->GetNumberOfTracklets();
935  Double_t BField = aod->GetMagneticField();
936  //loop over tracklets
937  for (Int_t itracklet=0; itracklet<multSPD; ++itracklet) {
938  Float_t thetaTr= anInputTracklets->GetTheta(itracklet);
939  Float_t phiTr= anInputTracklets->GetPhi(itracklet);
940  // calculate eta
941  Float_t etaTr = -TMath::Log(TMath::Tan(thetaTr/2.));
942  //make new AliFLowTrackSimple
943  AliFlowTrack* pTrack = new AliFlowTrack();
944  pTrack->SetPt(0.5);
945  pTrack->SetEta(etaTr);
946  // 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."
947  Double_t DeltaPhi = anInputTracklets->GetDeltaPhi(itracklet);
948  if(BField>0. && DeltaPhi>0.) pTrack->SetCharge(1);
949  if(BField>0. && DeltaPhi<0.) pTrack->SetCharge(-1);
950  if(BField<0. && DeltaPhi>0.) pTrack->SetCharge(-1);
951  if(BField<0. && DeltaPhi<0.) pTrack->SetCharge(1);
952  // correction of phi
954  phiTr += 39./34.*DeltaPhi;
955  if (phiTr < 0.) phiTr += 2.*TMath::Pi();
956  if (phiTr > 2.*TMath::Pi()) phiTr -= 2.*TMath::Pi();
957  }
958  pTrack->SetPhi(phiTr);
959  //marking the particles as POI type 2:
961  pTrack->SetPOItype(2,kTRUE);
963  //Add the track to the flowevent
964  fFlowEvent->AddTrack(pTrack);
965  }
966  }
967 
969 
974  fFlowEvent->SetCentralityCL1(centrCL1);
975  fFlowEvent->SetCentralityTRK(centrTRK);
976  // fFlowEvent->SetNITSCL1(((AliVAODHeader*)aod->GetHeader())->GetNumberOfITSClusters(1));
977  Double_t SumV0=0.;
978  for(Int_t i=0; i<64; i++) {
979  if(std::isfinite(aod->GetVZEROEqMultiplicity(i))) SumV0 += aod->GetVZEROEqMultiplicity(i);
980  }
981  fFlowEvent->SetNITSCL1(SumV0);
982 
983  // set absolute orbit number
984  UInt_t period = aod->GetPeriodNumber();
985  UInt_t orbit24 = aod->GetOrbitNumber(); // wrapped down to 24 bits
986  if (period > 255) { // 8 bits
987  period = 255;
988  orbit24 = (1<<24) - 1;
989  }
990  if (orbit24 >= (1<<24)) { // 24 bits
991  period = 255;
992  orbit24 = (1<<24) - 1;
993  }
994  UInt_t orbit = period * (1<<24) + orbit24;
995  fFlowEvent->SetAbsOrbit(orbit);
996 
997  Double_t vtxpos[3]={0.,0.,0.};
998  vtxpos[0] = ((AliAODVertex*)aod->GetPrimaryVertex())->GetX();
999  vtxpos[1] = ((AliAODVertex*)aod->GetPrimaryVertex())->GetY();
1000  vtxpos[2] = ((AliAODVertex*)aod->GetPrimaryVertex())->GetZ();
1001  fFlowEvent->SetVertexPosition(vtxpos);
1002 
1003  if (McEvent && McEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(McEvent);
1004 
1005  // run-by-run QA
1006  // for(Int_t jTracks = 0; jTracks<aod->GetNumberOfTracks(); jTracks++){
1007  // AliAODTrack* track = (AliAODTrack*)aod->GetTrack(jTracks);
1008  // if(!track) continue;
1009  // // general kinematic & quality cuts
1010  // if (track->Pt() < .2 || track->Pt() > 20. || TMath::Abs(track->Eta()) > .8 || track->GetTPCNcls() < 70) continue;
1011  // if (track->TestFilterBit(128)) fPtPhiEtaRbRFB128[RunBin][CenBin]->Fill(track->Pt(),track->Phi(),track->Eta());
1012  // if (track->TestFilterBit(768)) fPtPhiEtaRbRFB768[RunBin][CenBin]->Fill(track->Pt(),track->Phi(),track->Eta());
1013  // }
1014  fCenDis->Fill(centrV0M);
1015 
1016  }
1017 
1018  if (fAnalysisType == kTrackQA) {
1019 
1020  // empty flow event
1021  fFlowEvent->ClearFast();
1022 
1023  // get centrality
1024  Double_t centrV0M=300;
1026  centrV0M = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("V0M");
1027  } else {
1028  fMultSelection = (AliMultSelection*) InputEvent()->FindListObject("MultSelection");
1029  if(!fMultSelection) {
1030  //If you get this warning (and lPercentiles 300) please check that the AliMultSelectionTask actually ran (before your task)
1031  AliWarning("AliMultSelection object not found!");
1032  } else {
1033  centrV0M = fMultSelection->GetMultiplicityPercentile("V0M");
1034  }
1035  }
1036  if(centrV0M<10. || centrV0M>40.) return;
1037 
1038  //check event cuts
1039  if (InputEvent()) {
1040  if(!fCutsEvent->IsSelected(InputEvent(),MCEvent())) return;
1041  if(fRejectPileUp) {
1042  Bool_t IsPileUp = SelectPileup(aod);
1043  if(IsPileUp) return;
1044  }
1045  }
1046 
1047  Double_t BField = aod->GetMagneticField();
1048 
1049  for(Int_t jTracks = 0; jTracks<aod->GetNumberOfTracks(); jTracks++){
1050 
1051  AliAODTrack* track = (AliAODTrack*)aod->GetTrack(jTracks);
1052  if(!track) continue;
1053 
1054  Double_t dPt = track->Pt();
1055  Double_t dEta = track->Eta();
1056  Double_t dPhi = track->Phi();
1057  Double_t dCharge = track->Charge();
1058  Int_t cw = 0;
1059  if(BField>0. && dCharge>0.) cw=0;
1060  if(BField>0. && dCharge<0.) cw=1;
1061  if(BField<0. && dCharge>0.) cw=2;
1062  if(BField<0. && dCharge<0.) cw=3;
1063 
1064  // general kinematic cuts
1065  if (dPt < .2 || dPt > 50. || TMath::Abs(dEta) > 0.8) continue;
1066 
1067  // cut on DCA
1068  Double_t DCAxy = track->DCA();
1069  Double_t DCAz = track->ZAtDCA();
1070  if (std::abs((Int_t)DCAxy)==999 || std::abs((Int_t)DCAz)==999) {
1071  // re-evaluate the dca as it seems to not be natively present
1072  // allowed only for tracks inside the beam pipe
1073  Double_t pos[3] = {-99., -99., -99.};
1074  track->GetPosition(pos);
1075  if(pos[0]*pos[0]+pos[1]*pos[1] <= 3.*3.) {
1076  AliAODTrack copy(*track); // stack copy
1077  Double_t b[2] = {-99., -99.};
1078  Double_t bCov[3] = {-99., -99., -99.};
1079  if(copy.PropagateToDCA(aod->GetPrimaryVertex(), aod->GetMagneticField(), 100., b, bCov)) {
1080  DCAxy = b[0];
1081  DCAz = b[1];
1082  }
1083  }
1084  }
1085  if(fabs(DCAxy)>2.4 || fabs(DCAz)>3.2) continue;
1086 
1087  // various cuts on TPC clusters
1088  if (track->GetTPCNcls() < 70) continue;
1089  Double_t chi2_per_tpc = track->Chi2perNDF();
1090  if (chi2_per_tpc < 0.1 || chi2_per_tpc > 4.) continue;
1091  Double_t fraction_shared_tpccls = 1.*track->GetTPCnclsS()/track->GetTPCncls();
1092  if (fraction_shared_tpccls > 0.4) continue;
1093 
1094  Int_t FBarray[4] = {32,96,128,768};
1095 
1096  // test filter bits
1097  for(Int_t fb=0; fb<fKNFBs; fb++){
1098  if (track->TestFilterBit(FBarray[fb])) {
1099  fTrackQADCAxy[fb][cw]->Fill(dEta,dPt,DCAxy);
1100  fTrackQADCAz[fb][cw]->Fill(dEta,dPt,DCAz);
1101  fTrackQApT[fb][cw]->Fill(dEta,dPt);
1102  fEbEQRe[fb][cw]->Fill(dEta,dPt,TMath::Cos(dPhi));
1103  fEbEQIm[fb][cw]->Fill(dEta,dPt,TMath::Sin(dPhi));
1104  fEbEQMu[fb][cw]->Fill(dEta,dPt);
1105  }
1106  }
1107  }
1108 
1109  // compute cos(#Delta#phi)
1110  for(Int_t bx=1; bx<=fEbEQRe[0][0]->GetXaxis()->GetNbins(); bx++) {
1111  for(Int_t by=1; by<=fEbEQRe[0][0]->GetYaxis()->GetNbins(); by++) {
1112 
1113  Double_t dEta = fEbEQRe[0][0]->GetXaxis()->GetBinCenter(bx);
1114  Double_t dPt = fEbEQRe[0][0]->GetYaxis()->GetBinCenter(by);
1115 
1116  for(Int_t fb=0; fb<fKNFBs; fb++){
1117  for(Int_t c=0; c<4; c++){
1118 
1119  Double_t QRe = fEbEQRe[fb][c]->GetBinContent(bx,by);
1120  Double_t QIm = fEbEQIm[fb][c]->GetBinContent(bx,by);
1121  Double_t M = 1.*fEbEQMu[fb][c]->GetBinContent(bx,by);
1122 
1123  if(M>1.) {
1124  Double_t c1 = (QRe*QRe+QIm*QIm-M)/(M*(M-1.));
1125  fTrackQADphi[fb][c]->Fill(dEta,dPt,c1);
1126  }
1127  }
1128  }
1129 
1130  }
1131  }
1132 
1133  // reset EbE histograms
1134  for(Int_t fb=0; fb<fKNFBs; fb++){
1135  for(Int_t c=0; c<4; c++){
1136  fEbEQRe[fb][c]->Reset();
1137  fEbEQIm[fb][c]->Reset();
1138  fEbEQMu[fb][c]->Reset();
1139  }
1140  }
1141 
1142  }
1143 
1144  if (fAnalysisType == kMCAOD) {
1145 
1146  //check event cuts
1147  if (InputEvent()) {
1148  if(!fCutsEvent->IsSelected(InputEvent(),MCEvent())) return;
1149  if(fRejectPileUp && fAnalysisUtil->IsPileUpEvent(InputEvent())) return;
1150  }
1151 
1152  fFlowEvent->ClearFast();
1153 
1154  if(!McEvent) {
1155  AliError("ERROR: Could not retrieve MCEvent");
1156  return;
1157  }
1158  fStack = (TClonesArray*)aod->FindListObject(AliAODMCParticle::StdBranchName());
1159  if(!fStack){
1160  AliError("ERROR: Could not retrieve MCStack");
1161  return;
1162  }
1163 
1164  // get centrality (from AliMultSelection or AliCentrality)
1165  Double_t centr = 300;
1167  fMultSelection = (AliMultSelection*)aod->FindListObject("MultSelection");
1168  if(!fMultSelection) {
1169  //If you get this warning (and lPercentiles 300) please check that the AliMultSelectionTask actually ran (before your task)
1170  AliWarning("AliMultSelection object not found!");
1171  } else {
1172  centr = fMultSelection->GetMultiplicityPercentile("V0M");
1173  }
1174  } else {
1175  centr = (((AliVAODHeader*)aod->GetHeader())->GetCentralityP())->GetCentralityPercentile("V0M");
1176  }
1177  // centrality bin
1178  if (centr<fCentrLowLim || centr>=fCentrUpLim ) return;
1179  Int_t CenBin = -1;
1180  CenBin = GetCenBin(centr);
1181  if(CenBin==-1) return;
1182  fCenDis->Fill(centr);
1183 
1184  // reconstructed
1185  for(Int_t jTracks = 0; jTracks<aod->GetNumberOfTracks(); jTracks++){
1186 
1187  AliAODTrack* track = (AliAODTrack*)aod->GetTrack(jTracks);
1188  if(!track) continue;
1189 
1190  // to select primaries
1191  Int_t lp = TMath::Abs(track->GetLabel());
1192 
1193  // general kinematic cuts
1194  if (track->Pt() < .2 || track->Pt() > 20. || TMath::Abs(track->Eta()) > 0.8) continue;
1195 
1196  // cut on DCA
1197  Double_t DCAxy = track->DCA();
1198  Double_t DCAz = track->ZAtDCA();
1199  if(fabs(DCAxy)>2.4 || fabs(DCAz)>3.2) continue;
1200 
1201  // various cuts on TPC clusters
1202  if (track->GetTPCNcls() < 70) continue;
1203  Double_t chi2_per_tpc = track->Chi2perNDF();
1204  if (chi2_per_tpc < 0.1 || chi2_per_tpc > 4.) continue;
1205  Double_t fraction_shared_tpccls = 1.*track->GetTPCnclsS()/track->GetTPCncls();
1206  if (fraction_shared_tpccls > 0.4) continue;
1207 
1208  // test filter bits
1209  if (((AliAODMCParticle*)fStack->At(lp))->IsPhysicalPrimary()) {
1210  if (track->TestFilterBit(32)) fPtSpecFB32[0][CenBin]->Fill(track->Pt());
1211  if (track->TestFilterBit(96)) fPtSpecFB96[0][CenBin]->Fill(track->Pt());
1212  if (track->TestFilterBit(128)) fPtSpecFB128[0][CenBin]->Fill(track->Pt());
1213  if (track->TestFilterBit(768)) fPtSpecFB768[0][CenBin]->Fill(track->Pt());
1214  } else {
1215  if (track->TestFilterBit(32)) fPtSpecFB32[1][CenBin]->Fill(track->Pt());
1216  if (track->TestFilterBit(96)) fPtSpecFB96[1][CenBin]->Fill(track->Pt());
1217  if (track->TestFilterBit(128)) fPtSpecFB128[1][CenBin]->Fill(track->Pt());
1218  if (track->TestFilterBit(768)) fPtSpecFB768[1][CenBin]->Fill(track->Pt());
1219  }
1220 
1221  }
1222 
1223  // generated (physical primaries)
1224 
1225  for(Int_t jTracks = 0; jTracks<fStack->GetEntriesFast(); jTracks++) {
1226  AliAODMCParticle *MCpart = (AliAODMCParticle*)fStack->At(jTracks);
1227  if (!MCpart) {
1228  printf("ERROR: Could not receive MC track %d\n", jTracks);
1229  continue;
1230  }
1231  // kinematic cuts
1232  if ( MCpart->Pt() < .2 || MCpart->Pt() > 20. || TMath::Abs(MCpart->Eta()) > .8 ) continue;
1233  // select charged primaries
1234  if ( MCpart->Charge() == 0. || !MCpart->IsPhysicalPrimary()) continue;
1235 
1236  fPtSpecGen[0][CenBin]->Fill(MCpart->Pt());
1237  }
1238 
1239  // fGenHeader = McEvent->GenEventHeader();
1240  // if(fGenHeader) fPythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(fGenHeader);
1241  // printf("#reconstructed : %d (rejected from cuts %d), #MC primaries : %d (rejected from cuts %d) \n",AODPOIs,AODbads,MCPrims,MCSecos);
1242  fFlowEvent->SetReferenceMultiplicity(aod->GetNumberOfTracks());
1243  fFlowEvent->SetCentrality(centr);
1244  // if (McEvent && McEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(McEvent);
1245  fFlowEvent->SetRun(RunNum);
1246  // printf("Run : %d, RefMult : %d, Cent : %f \n",fFlowEvent->GetRun(),fFlowEvent->GetReferenceMultiplicity(),fFlowEvent->GetCentrality());
1247  }
1248 
1249  if(fAnalysisType == kMCESD) {
1250 
1251  fFlowEvent->ClearFast();
1252 
1253  if(!esd) {
1254  AliError("ERROR: Could not retrieve ESDEvent");
1255  return;
1256  }
1257  if(!McEvent) {
1258  AliError("ERROR: Could not retrieve MCEvent");
1259  return;
1260  }
1261  AliStack* fStack = fMCEvent->Stack();
1262  if(!fStack) {
1263  AliError("ERROR: Could not retrieve MCStack");
1264  return;
1265  }
1266 
1267  AliESDVertex *vertex = (AliESDVertex*) esd->GetPrimaryVertex();
1268  if (!vertex) return;
1269  if (TMath::Abs(vertex->GetZ()) > 10. ) return;
1270  if (vertex->GetNContributors() < 1 ) return;
1271  AliCentrality *centrality = esd->GetCentrality();
1272  if (!centrality) return;
1273  Double_t centr = centrality->GetCentralityPercentile("V0M");
1274  if (centr<fCentrLowLim || centr>=fCentrUpLim ) return;
1275  Int_t CenBin = -1;
1276  if (centr>0. && centr<5.) CenBin=0;
1277  if (centr>5. && centr<10.) CenBin=1;
1278  if (centr>10. && centr<20.) CenBin=2;
1279  if (centr>20. && centr<30.) CenBin=3;
1280  if (centr>30. && centr<40.) CenBin=4;
1281  if (centr>40. && centr<50.) CenBin=5;
1282  if (centr>50. && centr<60.) CenBin=6;
1283  if (centr>60. && centr<70.) CenBin=7;
1284  if (centr>70. && centr<80.) CenBin=8;
1285  if (centr>80. && centr<90.) CenBin=9;
1286  if(CenBin==-1) return;
1287 
1288  //Generated
1289  Int_t MCPrims = 0;
1290  for ( Int_t i=0 ; i<fStack->GetNtrack() ; i++ ) {
1291 
1292  //Primaries Selection
1293  TParticle *particle = (TParticle*)fStack->Particle(i);
1294  if (!particle) continue;
1295  if (!fStack->IsPhysicalPrimary(i)) continue;
1296  if ( particle->GetPDG()->Charge() == 0.) continue;
1297 
1298  //Kinematic Cuts
1299  if ( particle->Pt()<0.2 || particle->Pt()>10. ) continue;
1300  if ( TMath::Abs(particle->Eta())>0.8 ) continue;
1301 
1302  fFlowTrack->SetPhi(particle->Phi());
1303  fFlowTrack->SetEta(particle->Eta());
1304  fFlowTrack->SetPt(particle->Pt());
1306  fFlowTrack->SetForRPSelection(kTRUE);
1308  fFlowTrack->SetForPOISelection(kFALSE);
1310  MCPrims++;
1311 
1312  fPtSpecGen[0][CenBin]->Fill(particle->Pt());
1313 
1314  }
1315 
1316  //Reconstructed
1317  Int_t ESDPrims = 0;
1318  for (Int_t i=0 ; i<esd->GetNumberOfTracks() ; i++) {
1319 
1320  //Get reconstructed track
1321  AliVTrack *vtrack = static_cast<AliVTrack*>(esd->GetTrack(i));
1322  AliESDtrack *track = dynamic_cast<AliESDtrack*>(vtrack);
1323  if (!track) continue;
1324 
1325  //Primaries selection
1326  Int_t lp = TMath::Abs(track->GetLabel());
1327  if (!fStack->IsPhysicalPrimary(lp)) continue;
1328  TParticle *particle = (TParticle*)fStack->Particle(lp);
1329  if (!particle) continue;
1330  if (particle->GetPDG()->Charge() == 0.) continue;
1331 
1332  // if(!fCutsPOI->PassesESDcuts(track)) continue;
1333 
1334  Bool_t pass = kTRUE;
1335 
1336  if(fCutTPC) {
1337  // printf("******* cutting TPC ******** \n");
1338  UShort_t ntpccls = track->GetTPCNcls();
1339  Double_t tpcchi2 = track->GetTPCchi2();
1340  if (tpcchi2<0.2 || tpcchi2 >=4.) {
1341  // printf("TPCchi2 : %e %e ",tpcchi2,track->GetTPCchi2Iter1());
1342  pass=kFALSE;
1343  }
1344  if (ntpccls < 70) {
1345  // printf("#TPCcluster : %u %u %u %u ",ntpccls,track->GetTPCNclsF(),track->GetTPCNclsFIter1(),track->GetTPCNclsIter1());
1346  pass=kFALSE;
1347  }
1348  }
1349 
1350  Float_t dcaxy=0.0;
1351  Float_t dcaz=0.0;
1352  track->GetImpactParameters(dcaxy,dcaz);
1353  if (dcaxy > 0.3 || dcaz > 0.3) {
1354  // printf("DCA : %e %e ",dcaxy,dcaz);
1355  pass=kFALSE;
1356  }
1357  if(!pass) continue;
1358 
1359  //Kinematic Cuts
1360  if ( track->Pt()<0.2 || track->Pt()>10. ) continue;
1361  if ( TMath::Abs(track->Eta())>0.8 ) continue;
1362 
1363  fFlowTrack->SetPhi(track->Phi());
1364  fFlowTrack->SetEta(track->Eta());
1365  fFlowTrack->SetPt(track->Pt());
1367  fFlowTrack->SetForRPSelection(kFALSE);
1371  ESDPrims++;
1372 
1373  }
1374 
1375  // printf("#reconstructed : %d , #MC primaries : %d \n",ESDPrims,MCPrims);
1376  fFlowEvent->SetReferenceMultiplicity(esd->GetNumberOfTracks());
1377  fFlowEvent->SetCentrality(centr);
1378  if (McEvent && McEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(McEvent);
1379  fFlowEvent->SetRun(esd->GetRunNumber());
1380  // printf("Run : %d, RefMult : %d, Cent : %f \n",fFlowEvent->GetRun(),fFlowEvent->GetReferenceMultiplicity(),fFlowEvent->GetCentrality());
1381 
1382  } // end of if(fAnalysisType == kMCESD)
1383 
1384  if(fAnalysisType == kMCkine) {
1385 
1386  fFlowEvent->ClearFast();
1387 
1388  AliInputEventHandler* McHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
1389  if(!McHandler) {
1390  AliError("ERROR: Could not retrieve MCtruthEventHandler");
1391  return;
1392  }
1393  McEvent = McHandler->MCEvent();
1394  if(!McEvent) {
1395  AliError("ERROR: Could not retrieve MC event");
1396  return;
1397  }
1398 
1399  Int_t nTracks = McEvent->GetNumberOfTracks();
1400  // Int_t nPrimTr = McEvent->GetNumberOfPrimaries();
1401 
1402  //loop over tracks
1403  for (Int_t itrkN=0; itrkN<nTracks; itrkN++) {
1404  //get input particle
1405  AliMCParticle* pParticle = dynamic_cast<AliMCParticle*>(McEvent->GetTrack(itrkN));
1406  if (!pParticle) continue;
1407 
1408  //check if track passes the cuts
1409  if (McEvent->IsPhysicalPrimary(itrkN) && pParticle->Charge()!=0) {
1410  fFlowTrack->Set(pParticle);
1412  fFlowTrack->SetForRPSelection(kTRUE);
1417  }
1418  }// for all tracks
1419 
1420  // if monte carlo event get reaction plane from monte carlo (depends on generator)
1421  if (McEvent && McEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(McEvent);
1422  // set reference multiplicity
1423  fFlowEvent->SetReferenceMultiplicity(McEvent->GetNumberOfTracks());
1424  // tag subevents
1426  // set centrality from impact parameter
1427  Double_t ImpPar=0., CenPer=0.;
1428  fGenHeader = McEvent->GenEventHeader();
1429  if(fGenHeader){
1430  fPythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(fGenHeader);
1431  if(fPythiaGenHeader) ImpPar = fPythiaGenHeader->GetImpactParameter();
1432  fHijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(fGenHeader);
1433  if(fHijingGenHeader) ImpPar = fHijingGenHeader->ImpactParameter();
1434  if(ImpPar) CenPer = 0.3859796743103508*pow(ImpPar,2.);
1435  if(CenPer>0. && CenPer<100.) fFlowEvent->SetCentrality(CenPer);
1436  else return;
1437  fFlowEvent->SetRun(1);
1438  }
1439 
1440  } // end of if(fAnalysisType == kMCkine)
1441 
1442  if (!fFlowEvent) return; //shuts up coverity
1443 
1444  //check final event cuts
1445  Int_t mult = fFlowEvent->NumberOfTracks();
1446  // AliInfo(Form("FlowEvent has %i tracks",mult));
1447  if (mult<fMinMult || mult>fMaxMult) {
1448  AliWarning("FlowEvent cut on multiplicity"); return;
1449  }
1450 
1451  //define dead zone
1453 
1456  if (fAfterburnerOn)
1457  {
1458  //if reaction plane not set from elsewhere randomize it before adding flow
1460  fFlowEvent->SetMCReactionPlaneAngle(gRandom->Uniform(0.0,TMath::TwoPi()));
1461 
1462  if(fDifferentialV2)
1464  else
1465  fFlowEvent->AddFlow(fV1,fV2,fV3,fV4,fV5); //add flow
1466  fFlowEvent->CloneTracks(fNonFlowNumberOfTrackClones); //add nonflow by cloning tracks
1467  }
1469 
1470  //tag subEvents
1472 
1473  //do we want to serve shullfed tracks to everybody?
1475 
1476  // associate the mother particles to their daughters in the flow event (if any)
1478 
1479  //fListHistos->Print();
1480  //fOutputFile->WriteObject(fFlowEvent,"myFlowEventSimple");
1481 
1482  //********************************************************************************************************************************
1483 
1485 
1486  // PHYSICS SELECTION
1487  AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1488  AliInputEventHandler *hdr = (AliInputEventHandler*)am->GetInputEventHandler();
1489 
1490  if(hdr->IsEventSelected() && AliVEvent::kAny) {
1491 
1492  Double_t centrperc = fFlowEvent->GetCentrality();
1493  Int_t cenb = (Int_t)centrperc;
1494 
1495  AliAODTracklets *trackl = aod->GetTracklets();
1496  Int_t nTracklets = trackl->GetNumberOfTracklets();
1497 
1498  // get VZERO data
1499  AliAODVZERO *vzeroAOD = aod->GetVZEROData();
1500  Double_t multV0A = vzeroAOD->GetMTotV0A();
1501  Double_t multV0C = vzeroAOD->GetMTotV0C();
1502 
1503  // set VZERO Q-vectors
1504  if(fDataSet==k2015 || fDataSet==k2015v6) {
1505  Int_t CachednRing = 1;
1506  Double_t QxTot[fkVZEROnHar] = {0.}, QyTot[fkVZEROnHar] = {0.};
1507  Double_t denom = 0.;
1508  Double_t V0TotQC[fkVZEROnHar][2] = {{0.}}, V0TotQA[fkVZEROnHar][2] = {{0.}};
1509  Double_t MultC[fkVZEROnHar] = {0.}, MultA[fkVZEROnHar] = {0.};
1510 
1511  for(Int_t i=0; i<64; i++) {
1512 
1513  // correct multiplicity per channel
1514  Double_t mult = vzeroAOD->GetMultiplicity(i);
1515  if(fVZEROGainEqHist) {
1516  Double_t EqFactor = fVZEROGainEqHist->GetBinContent(RunBin+1,i+1);
1517  if(EqFactor>0.) mult *= EqFactor;
1518  }
1519  fVZEROMult->Fill(RunBin+0.5,i+0.5,mult);
1520 
1521  // build Q-vector per ring
1522  Int_t nRing = (Int_t)i/8 + 1;
1523  Double_t ChPhi = TMath::PiOver4()*(0.5+i%8);
1524 
1525  if(i == 63) {
1526  for (Int_t k=0; k<fkVZEROnHar; k++) {
1527  QxTot[k] += mult*TMath::Cos((k+1.)*ChPhi);
1528  QyTot[k] += mult*TMath::Sin((k+1.)*ChPhi);
1529  }
1530  denom += mult;
1531  nRing++;
1532  }
1533 
1534  if(nRing!=CachednRing) {
1535  for (Int_t k=0; k<fkVZEROnHar; k++) {
1536  Double_t QxRec = QxTot[k]/denom;
1537  Double_t QyRec = QyTot[k]/denom;
1538  // store values for re-centering
1539  // fVZEROQVectorRecQx[k]->Fill(RunBin+0.5,centrperc,CachednRing-0.5,QxRec);
1540  // fVZEROQVectorRecQy[k]->Fill(RunBin+0.5,centrperc,CachednRing-0.5,QyRec);
1541  // do re-centering
1542  if(fVZEROQVectorRecQxStored[k]) {
1543  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));
1544  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));
1545  }
1546  // sum of Q-vectors over all rings (total V0 Q-vector)
1547  if (CachednRing >= fMinRingVZC && CachednRing <= fMaxRingVZC) {
1548  V0TotQC[k][0] += QxRec*denom;
1549  V0TotQC[k][1] += QyRec*denom;
1550  MultC[k] += denom;
1551  }
1552  if (CachednRing >= fMinRingVZA && CachednRing <= fMaxRingVZA) {
1553  V0TotQA[k][0] += QxRec*denom;
1554  V0TotQA[k][1] += QyRec*denom;
1555  MultA[k] += denom;
1556  }
1557  QxTot[k] = 0.;
1558  QyTot[k] = 0.;
1559  }
1560  denom = 0.;
1561  CachednRing = nRing;
1562  }
1563  for (Int_t k=0; k<fkVZEROnHar; k++) {
1564  QxTot[k] += mult*TMath::Cos((k+1.)*ChPhi);
1565  QyTot[k] += mult*TMath::Sin((k+1.)*ChPhi);
1566  }
1567  denom += mult;
1568  }
1569 
1570  for (Int_t k=0; k<fkVZEROnHar; k++) {
1571  if(MultC[k]>0. && MultA[k]>0.) {
1572  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];
1573  if(!std::isnan(QCx) && !std::isnan(QCy) && !std::isnan(QAx) && !std::isnan(QAy)) {
1574  fFlowEvent->SetV02Qsub(QCx,QCy,MultC[k],QAx,QAy,MultA[k],k+1);
1575  fVZEROQVectorRecFinal[k][0]->Fill(RunBin+0.5,centrperc,QCx);
1576  fVZEROQVectorRecFinal[k][1]->Fill(RunBin+0.5,centrperc,QCy);
1577  fVZEROQVectorRecFinal[k][2]->Fill(RunBin+0.5,centrperc,QAx);
1578  fVZEROQVectorRecFinal[k][3]->Fill(RunBin+0.5,centrperc,QAy);
1579  fVZEROQVectorRecFinal[k][4]->Fill(RunBin+0.5,centrperc,QCx*QAx);
1580  fVZEROQVectorRecFinal[k][5]->Fill(RunBin+0.5,centrperc,QCy*QAy);
1581  fVZEROQVectorRecFinal[k][6]->Fill(RunBin+0.5,centrperc,QCx*QAy);
1582  fVZEROQVectorRecFinal[k][7]->Fill(RunBin+0.5,centrperc,QCy*QAx);
1583  } else {
1584  fFlowEvent->SetV02Qsub(0.,0.,0.,0.,0.,0.,k+1);
1585  }
1586  } else {
1587  fFlowEvent->SetV02Qsub(0.,0.,0.,0.,0.,0.,k+1);
1588  }
1589  }
1590  }
1591 
1592 // AliAODForwardMult* aodForward = static_cast<AliAODForwardMult*>(aodEvent->FindListObject("Forward"));
1593 // const TH2D& d2Ndetadphi = aodForward->GetHistogram();
1594 // Int_t nEta = d2Ndetadphi.GetXaxis()->GetNbins();
1595 // Int_t nPhi = d2Ndetadphi.GetYaxis()->GetNbins();
1596 // Double_t ret = 0.;
1597 // // Loop over eta
1598 // for (Int_t iEta = 1; iEta <= nEta; iEta++) {
1599 // Int_t valid = d2Ndetadphi.GetBinContent(iEta, 0);
1600 // if (!valid) continue; // No data expected for this eta
1601 // // Loop over phi
1602 // for (Int_t iPhi = 1; i <= nPhi; i++) {
1603 // ret = d2Ndetadphi.GetBinContent(iEta, iPhi);
1604 // printf("eta %e phi %e : %e \n",d2Ndetadphi.GetXaxis()->GetBinCenter(iEta),d2Ndetadphi.GetYaxis()->GetBinCenter(iPhi),ret);
1605 // }
1606 // }
1607 
1608  AliAODZDC *aodZDC = aod->GetZDCData();
1609 
1610  Double_t energyZNC = (Double_t) (aodZDC->GetZNCEnergy());
1611  Double_t energyZPC = (Double_t) (aodZDC->GetZPCEnergy());
1612  Double_t energyZNA = (Double_t) (aodZDC->GetZNAEnergy());
1613  Double_t energyZPA = (Double_t) (aodZDC->GetZPAEnergy());
1614  Double_t energyZEM1 = (Double_t) (aodZDC->GetZEM1Energy());
1615  Double_t energyZEM2 = (Double_t) (aodZDC->GetZEM2Energy());
1616 
1617  const Double_t * towZNCraw = aodZDC->GetZNCTowerEnergy();
1618  const Double_t * towZNAraw = aodZDC->GetZNATowerEnergy();
1619 
1620  // Get centroid from ZDCs *******************************************************
1621 
1622  Double_t Enucl = (RunNum < 209122 ? 1380. : 2511.);
1623  Double_t xyZNC[2]={0.,0.}, xyZNA[2]={0.,0.};
1624  Double_t towZNC[5]={0.}, towZNA[5]={0.};
1625 
1626  Double_t ZNCcalib=1., ZNAcalib=1.;
1627  if(fUseTowerEq) {
1628  if(RunNum!=fCachedRunNum) {
1629  for(Int_t i=0; i<5; i++) {
1630  fTowerGainEq[0][i] = (TH1D*)(fTowerEqList->FindObject(Form("fZNCTower[%d][%d]",RunNum,i)));
1631  fTowerGainEq[1][i] = (TH1D*)(fTowerEqList->FindObject(Form("fZNATower[%d][%d]",RunNum,i)));
1632  }
1633  }
1634  for(Int_t i=0; i<5; i++) {
1635  if(fTowerGainEq[0][i]) towZNC[i] = towZNCraw[i]*fTowerGainEq[0][i]->GetBinContent(fTowerGainEq[0][i]->FindBin(centrperc));
1636  if(fTowerGainEq[1][i]) towZNA[i] = towZNAraw[i]*fTowerGainEq[1][i]->GetBinContent(fTowerGainEq[1][i]->FindBin(centrperc));
1637  if(fResetNegativeZDC) {
1638  if(towZNC[i]<0.) towZNC[i] = 0.;
1639  if(towZNA[i]<0.) towZNA[i] = 0.;
1640  }
1641  }
1642  } else {
1643  for(Int_t i=0; i<5; i++) {
1644  towZNC[i] = towZNCraw[i];
1645  towZNA[i] = towZNAraw[i];
1646  if(fResetNegativeZDC) {
1647  if(towZNC[i]<0.) towZNC[i] = 0.;
1648  if(towZNA[i]<0.) towZNA[i] = 0.;
1649  }
1650  fZNCTower[RunBin][i]->Fill(centrperc,towZNC[i]);
1651  fZNATower[RunBin][i]->Fill(centrperc,towZNA[i]);
1652  }
1653  }
1654 
1655  if(RunNum>=245829) towZNA[2] = 0.;
1656  Double_t zncEnergy=0., znaEnergy=0.;
1657  for(Int_t i=0; i<5; i++){
1658  zncEnergy += towZNC[i];
1659  znaEnergy += towZNA[i];
1660  }
1661  if(RunNum>=245829) znaEnergy *= 8./7.;
1662  fFlowEvent->SetZNCEnergy(towZNC[0]);
1663  fFlowEvent->SetZNAEnergy(towZNA[0]);
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  Double_t tdcSum = aodZDC->GetZDCTimeSum();
1800  Double_t tdcDiff = aodZDC->GetZDCTimeDiff();
1801  fhDebunch->Fill(tdcDiff, tdcSum);
1802 
1803  for(int i=0; i<5; i++){
1804  fhZNCPM[i]->Fill(towZNC[i]);
1805  if((i<4) && (towZNC[0]>0.)) fhZNCPMQiPMC[i]->Fill(towZNC[i+1]/towZNC[0]);
1806  }
1807  for(int i=0; i<5; i++){
1808  fhZNAPM[i]->Fill(towZNA[i]);
1809  if(((i<4) && towZNA[0]>0.)) fhZNAPMQiPMC[i]->Fill(towZNA[i+1]/towZNA[0]);
1810  }
1811 
1812  fhZNCvsZNA->Fill(energyZNA, energyZNC);
1813  fhZDCCvsZDCCA->Fill(energyZNA+energyZPA, energyZNC+energyZPC);
1814  fhZNCvsZPC->Fill(energyZPC, energyZNC);
1815  fhZNAvsZPA->Fill(energyZPA, energyZNA);
1816  fhZNvsZP->Fill(energyZPA+energyZPC, energyZNA+energyZNC);
1817  fhZNvsVZERO->Fill(multV0A+multV0C, energyZNC+energyZNA);
1818  fhZDCvsVZERO->Fill(multV0A+multV0C, energyZNA+energyZPA+energyZNC+energyZPC);
1819  fhZDCvsTracklets->Fill((Double_t) (nTracklets), energyZNA+energyZPA+energyZNC+energyZPC);
1820 
1821  Double_t asymmetry = -999.;
1822  if((energyZNC+energyZNA)>0.) asymmetry = (energyZNC-energyZNA)/(energyZNC+energyZNA);
1823  fhAsymm->Fill(asymmetry);
1824  fhZNAvsAsymm->Fill(asymmetry, energyZNA/1000.);
1825  fhZNCvsAsymm->Fill(asymmetry, energyZNC/1000.);
1826 
1827  fhZNCvscentrality->Fill(centrperc, energyZNC/1000.);
1828  fhZNAvscentrality->Fill(centrperc, energyZNA/1000.);
1829 
1830  } // PHYSICS SELECTION
1831 
1832  }
1833 
1834  // p) cache run number
1836 
1837  // printf("debug: NoRPs %e, NoPOIs %e, RunNum %d, Cen %e \n",fFlowEvent->GetNumberOfRPs(),fFlowEvent->GetNumberOfPOIs(),fCachedRunNum,fFlowEvent->GetCentrality());
1838 
1839  PostData(1, fFlowEvent);
1840 
1841  PostData(2, fOutput);
1842 }
1843 //________________________________________________________________________
1844 
1846 {
1847  Bool_t BisPileup=kFALSE;
1848  Double_t centrV0M=300., centrCL1=300.;
1849 
1851 
1852  // pileup for LHC10h and LHC11h
1853 
1854  centrV0M = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("V0M");
1855  centrCL1 = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("CL1");
1856 
1857  // check anyway pileup
1858  if (plpMV(aod)) {
1859  fPileUpCount->Fill(0.5);
1860  BisPileup=kTRUE;
1861  }
1862 
1863  Short_t isPileup = aod->IsPileupFromSPD(3);
1864  if (isPileup != 0) {
1865  fPileUpCount->Fill(1.5);
1866  BisPileup=kTRUE;
1867  }
1868 
1869  if (((AliAODHeader*)aod->GetHeader())->GetRefMultiplicityComb08() < 0) {
1870  fPileUpCount->Fill(2.5);
1871  BisPileup=kTRUE;
1872  }
1873 
1874  if (aod->IsIncompleteDAQ()) {
1875  fPileUpCount->Fill(3.5);
1876  BisPileup=kTRUE;
1877  }
1878 
1879  // check vertex consistency
1880  const AliAODVertex* vtTrc = aod->GetPrimaryVertex();
1881  const AliAODVertex* vtSPD = aod->GetPrimaryVertexSPD();
1882 
1883  if (vtTrc->GetNContributors() < 2 || vtSPD->GetNContributors()<1) {
1884  fPileUpCount->Fill(5.5);
1885  BisPileup=kTRUE;
1886  }
1887 
1888  double covTrc[6], covSPD[6];
1889  vtTrc->GetCovarianceMatrix(covTrc);
1890  vtSPD->GetCovarianceMatrix(covSPD);
1891 
1892  double dz = vtTrc->GetZ() - vtSPD->GetZ();
1893 
1894  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
1895  double errTrc = TMath::Sqrt(covTrc[5]);
1896  double nsigTot = dz/errTot;
1897  double nsigTrc = dz/errTrc;
1898 
1899  if (TMath::Abs(dz)>0.2 || TMath::Abs(nsigTot)>10 || TMath::Abs(nsigTrc)>20) {
1900  fPileUpCount->Fill(6.5);
1901  BisPileup=kTRUE;
1902  }
1903 
1904  if (fAnalysisUtil->IsPileUpEvent(InputEvent())) {
1905  fPileUpCount->Fill(7.5);
1906  BisPileup=kTRUE;
1907  }
1908 
1909  }
1910  else {
1911 
1912  // pileup for LHC15o, using AliMultSelection
1913 
1914  if(fMultSelection) {
1915  centrV0M = fMultSelection->GetMultiplicityPercentile("V0M");
1916  centrCL1 = fMultSelection->GetMultiplicityPercentile("CL1");
1917  } else {
1918  BisPileup=kTRUE;
1919  }
1920 
1921  // pileup from AliMultSelection
1922  if(!fMultSelection->GetThisEventIsNotPileup()) fPileUpMultSelCount->Fill(0.5);
1923  if(!fMultSelection->GetThisEventIsNotPileupMV()) fPileUpMultSelCount->Fill(1.5);
1924  if(!fMultSelection->GetThisEventIsNotPileupInMultBins()) fPileUpMultSelCount->Fill(2.5);
1925  if(!fMultSelection->GetThisEventHasNoInconsistentVertices()) fPileUpMultSelCount->Fill(3.5);
1926  if(!fMultSelection->GetThisEventPassesTrackletVsCluster()) fPileUpMultSelCount->Fill(4.5);
1927  if(!fMultSelection->GetThisEventIsNotAsymmetricInVZERO()) fPileUpMultSelCount->Fill(5.5);
1928  if(!fMultSelection->GetThisEventIsNotIncompleteDAQ()) fPileUpMultSelCount->Fill(6.5);
1929  if(!fMultSelection->GetThisEventHasGoodVertex2016()) fPileUpMultSelCount->Fill(7.5);
1930 
1931  // pile-up a la Dobrin for LHC15o
1932  if (plpMV(aod)) {
1933  fPileUpCount->Fill(0.5);
1934  BisPileup=kTRUE;
1935  }
1936 
1937  Short_t isPileup = aod->IsPileupFromSPD(3);
1938  if (isPileup != 0) {
1939  fPileUpCount->Fill(1.5);
1940  BisPileup=kTRUE;
1941  }
1942 
1943  if (((AliAODHeader*)aod->GetHeader())->GetRefMultiplicityComb08() < 0) {
1944  fPileUpCount->Fill(2.5);
1945  BisPileup=kTRUE;
1946  }
1947 
1948  if (aod->IsIncompleteDAQ()) {
1949  fPileUpCount->Fill(3.5);
1950  BisPileup=kTRUE;
1951  }
1952 
1953  if(fabs(centrV0M-centrCL1)>7.5) {
1954  fPileUpCount->Fill(4.5);
1955  BisPileup=kTRUE;
1956  }
1957 
1958  // check vertex consistency
1959  const AliAODVertex* vtTrc = aod->GetPrimaryVertex();
1960  const AliAODVertex* vtSPD = aod->GetPrimaryVertexSPD();
1961 
1962  if (vtTrc->GetNContributors() < 2 || vtSPD->GetNContributors()<1) {
1963  fPileUpCount->Fill(5.5);
1964  BisPileup=kTRUE;
1965  }
1966 
1967  double covTrc[6], covSPD[6];
1968  vtTrc->GetCovarianceMatrix(covTrc);
1969  vtSPD->GetCovarianceMatrix(covSPD);
1970 
1971  double dz = vtTrc->GetZ() - vtSPD->GetZ();
1972 
1973  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
1974  double errTrc = TMath::Sqrt(covTrc[5]);
1975  double nsigTot = dz/errTot;
1976  double nsigTrc = dz/errTrc;
1977 
1978  if (TMath::Abs(dz)>0.2 || TMath::Abs(nsigTot)>10 || TMath::Abs(nsigTrc)>20) {
1979  fPileUpCount->Fill(6.5);
1980  BisPileup=kTRUE;
1981  }
1982 
1983  // cuts on tracks
1984  const Int_t nTracks = aod->GetNumberOfTracks();
1985  Int_t multEsd = ((AliAODHeader*)aod->GetHeader())->GetNumberOfESDTracks();
1986 
1987  Int_t multTrk = 0;
1988  Int_t multTrkBefC = 0;
1989  Int_t multTrkTOFBefC = 0;
1990  Int_t multTPC = 0;
1991 
1992  for (Int_t it = 0; it < nTracks; it++) {
1993 
1994  AliAODTrack* aodTrk = (AliAODTrack*)aod->GetTrack(it);
1995  if (!aodTrk){
1996  delete aodTrk;
1997  continue;
1998  }
1999 
2000 // if (aodTrk->TestFilterBit(32)){
2001 // multTrkBefC++;
2002 //
2003 // if ( TMath::Abs(aodTrk->GetTOFsignalDz()) <= 10. && aodTrk->GetTOFsignal() >= 12000. && aodTrk->GetTOFsignal() <= 25000.)
2004 // multTrkTOFBefC++;
2005 //
2006 // if ((TMath::Abs(aodTrk->Eta()) < 0.8) && (aodTrk->GetTPCNcls() >= 70) && (aodTrk->Pt() >= 0.2) && (aodTrk->Pt() < 20.))
2007 // multTrk++;
2008 // }
2009 
2010  if (aodTrk->TestFilterBit(128))
2011  multTPC++;
2012 
2013  } // end of for (Int_t it = 0; it < nTracks; it++)
2014 
2015  Double_t multTPCn = multTPC;
2016  Double_t multEsdn = multEsd;
2017  Double_t multESDTPCDif = multEsdn - multTPCn*3.38;
2018 
2019  if (multESDTPCDif > (fRejectPileUpTight?700.:15000.)) {
2020  fPileUpCount->Fill(7.5);
2021  BisPileup=kTRUE;
2022  }
2023 
2024  if(fRejectPileUpTight) {
2025  if(BisPileup==kFALSE) {
2026  if(!fMultSelection->GetThisEventIsNotPileup()) BisPileup=kTRUE;
2027  if(!fMultSelection->GetThisEventIsNotPileupMV()) BisPileup=kTRUE;
2028  if(!fMultSelection->GetThisEventIsNotPileupInMultBins()) BisPileup=kTRUE;
2029  if(!fMultSelection->GetThisEventHasNoInconsistentVertices()) BisPileup=kTRUE;
2030  if(!fMultSelection->GetThisEventPassesTrackletVsCluster()) BisPileup=kTRUE;
2031  if(!fMultSelection->GetThisEventIsNotIncompleteDAQ()) BisPileup=kTRUE;
2032  if(!fMultSelection->GetThisEventHasGoodVertex2016()) BisPileup=kTRUE;
2033  if(BisPileup) fPileUpCount->Fill(8.5);
2034  }
2035  }
2036  }
2037 
2038  return BisPileup;
2039 }
2040 
2041 //________________________________________________________________________
2042 
2044 {
2045  Double_t EtC = BadTowerCalibHist->ProjectionY("",BadTowerCalibHist->GetXaxis()->FindBin(Et),BadTowerCalibHist->GetXaxis()->FindBin(Et))->GetRandom();
2046  return EtC;
2047 }
2048 
2049 //________________________________________________________________________
2050 
2052 {
2053  Int_t CenBin=-1;
2054  if (Centrality>0. && Centrality<5.) CenBin=0;
2055  if (Centrality>5. && Centrality<10.) CenBin=1;
2056  if (Centrality>10. && Centrality<20.) CenBin=2;
2057  if (Centrality>20. && Centrality<30.) CenBin=3;
2058  if (Centrality>30. && Centrality<40.) CenBin=4;
2059  if (Centrality>40. && Centrality<50.) CenBin=5;
2060  if (Centrality>50. && Centrality<60.) CenBin=6;
2061  if (Centrality>60. && Centrality<70.) CenBin=7;
2062  if (Centrality>70. && Centrality<80.) CenBin=8;
2063  if (Centrality>80. && Centrality<90.) CenBin=9;
2064  if (CenBin>=fnCen) CenBin=-1;
2065  if (fnCen==1) CenBin=0;
2066  return CenBin;
2067 } // end of AliFlowAnalysisCRC::GetCRCCenBin(Double_t Centrality)
2068 //_____________________________________________________________________________
2069 
2070 Double_t AliAnalysisTaskCRCZDC::GetWDist(const AliVVertex* v0, const AliVVertex* v1)
2071 {
2072  // calculate sqrt of weighted distance to other vertex
2073  if (!v0 || !v1) {
2074  printf("One of vertices is not valid\n");
2075  return 0;
2076  }
2077  static TMatrixDSym vVb(3);
2078  double dist = -1;
2079  double dx = v0->GetX()-v1->GetX();
2080  double dy = v0->GetY()-v1->GetY();
2081  double dz = v0->GetZ()-v1->GetZ();
2082  double cov0[6],cov1[6];
2083  v0->GetCovarianceMatrix(cov0);
2084  v1->GetCovarianceMatrix(cov1);
2085  vVb(0,0) = cov0[0]+cov1[0];
2086  vVb(1,1) = cov0[2]+cov1[2];
2087  vVb(2,2) = cov0[5]+cov1[5];
2088  vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
2089  vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
2090  vVb.InvertFast();
2091  if (!vVb.IsValid()) {printf("Singular Matrix\n"); return dist;}
2092  dist = vVb(0,0)*dx*dx + vVb(1,1)*dy*dy + vVb(2,2)*dz*dz
2093  + 2*vVb(0,1)*dx*dy + 2*vVb(0,2)*dx*dz + 2*vVb(1,2)*dy*dz;
2094  return dist>0 ? TMath::Sqrt(dist) : -1;
2095 }
2096 //________________________________________________________________________
2097 
2099 {
2100  // check for multi-vertexer pile-up
2101 
2102  const int kMinPlpContrib = 5;
2103  const double kMaxPlpChi2 = 5.0;
2104  const double kMinWDist = 15;
2105 
2106  const AliVVertex* vtPrm = 0;
2107  const AliVVertex* vtPlp = 0;
2108  int nPlp = 0;
2109 
2110  if ( !(nPlp=aod->GetNumberOfPileupVerticesTracks()) ) return kFALSE;
2111  vtPrm = aod->GetPrimaryVertex();
2112  if (vtPrm == aod->GetPrimaryVertexSPD()) return kTRUE; // there are pile-up vertices but no primary
2113 
2114  //int bcPrim = vtPrm->GetBC();
2115 
2116  for (int ipl=0;ipl<nPlp;ipl++) {
2117  vtPlp = (const AliVVertex*)aod->GetPileupVertexTracks(ipl);
2118  //
2119  if (vtPlp->GetNContributors() < kMinPlpContrib) continue;
2120  if (vtPlp->GetChi2perNDF() > kMaxPlpChi2) continue;
2121  // int bcPlp = vtPlp->GetBC();
2122  // if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2) return kTRUE; // pile-up from other BC
2123  //
2124  double wDst = GetWDist(vtPrm,vtPlp);
2125  if (wDst<kMinWDist) continue;
2126  //
2127  return kTRUE; // pile-up: well separated vertices
2128  }
2129 
2130  return kFALSE;
2131 }
2132 
2133 //________________________________________________________________________
2135  fCutContainer->AddAt(cutsRP,0); fCutsRP=cutsRP; cutsRP->SetPOItype(0);
2136 }
2137 
2138 //________________________________________________________________________
2140  fCutContainer->AddAt(cutsPOI,1); fCutsPOI=cutsPOI; cutsPOI->SetPOItype(1);
2141 }
2142 
2143 //________________________________________________________________________
2145 {
2146  // Terminate analysis
2147  //
2148  /* if(fDebug > 1) printf(" **** AliAnalysisTaskCRCZDC::Terminate() \n");
2149 
2150  //fOutput = dynamic_cast<TList*> (GetOutputData(1));
2151  //if(!fOutput) printf("ERROR: fOutput not available\n");
2152  */
2153 }
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
TH2F * fhZDCvsTracklets
ZDC vs VZERO;.
static const Int_t fCRCnTow
TH1F * fhZNCPMQiPMC[4]
ZNA PM high res.
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)
TH2F * fhZNCvsAsymm
ZNA vs asymmetry.
void SetMCReactionPlaneAngle(const AliMCEvent *mcEvent)
TH3D * fhZNCenDis[2]
Debunch;.
void SetCutsPOI(AliFlowTrackCuts *cutsPOI)
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.
TH2F * fhDebunch
ZDC vs N_cluster layer 1;.
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.
TH2F * fhZDCvsNclu1
ZDC vs N_tracklets;.
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)
ClassImp(AliAnalysisTaskCRCZDC) AliAnalysisTaskCRCZDC
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
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)
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)
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