AliPhysics  8195daa (8195daa)
 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 "TProfile3D.h"
40 
41 #include "AliAnalysisManager.h"
42 #include "AliInputEventHandler.h"
43 #include "AliVEvent.h"
44 #include "AliESD.h"
45 #include "AliESDEvent.h"
46 #include "AliESDHeader.h"
47 #include "AliESDInputHandler.h"
48 #include "AliESDZDC.h"
49 #include "AliMultiplicity.h"
50 #include "AliAnalysisUtils.h"
51 #include "AliAODHandler.h"
52 #include "AliAODTrack.h"
53 #include "AliAODEvent.h"
54 #include "AliAODHeader.h"
55 #include "AliAODVertex.h"
56 #include "AliAODVZERO.h"
57 #include "AliAODZDC.h"
58 #include "AliAODMCHeader.h"
59 #include "AliMCEventHandler.h"
60 #include "AliMCEvent.h"
61 #include "AliHeader.h"
62 #include "AliVParticle.h"
63 #include "AliStack.h"
64 #include "AliAODMCParticle.h"
65 #include "AliAnalysisTaskSE.h"
66 #include "AliGenEventHeader.h"
67 #include "AliPhysicsSelectionTask.h"
68 #include "AliPhysicsSelection.h"
69 #include "AliBackgroundSelection.h"
70 #include "AliTriggerAnalysis.h"
71 #include "AliCentrality.h"
72 #include "AliAnalysisTaskCRCZDC.h"
73 #include "AliMultSelection.h"
74 #include "AliLumiTools.h"
75 
76 // ALICE Correction Framework
77 #include "AliCFManager.h"
78 
79 // Interface to Event generators to get Reaction Plane Angle
80 #include "AliGenCocktailEventHeader.h"
81 #include "AliGenPythiaEventHeader.h"
82 #include "AliGenHijingEventHeader.h"
83 #include "AliGenGeVSimEventHeader.h"
84 #include "AliGenEposEventHeader.h"
85 
86 // Interface to Load short life particles
87 #include "TObjArray.h"
88 #include "AliFlowCandidateTrack.h"
89 
90 // Interface to make the Flow Event Simple used in the flow analysis methods
91 #include "AliFlowEvent.h"
92 #include "AliFlowTrackCuts.h"
93 #include "AliFlowEventCuts.h"
94 #include "AliFlowCommonConstants.h"
95 
97 
98 //________________________________________________________________________
101 fAnalysisType("AUTOMATIC"),
102 fRPType(""),
103 fCFManager1(NULL),
104 fCFManager2(NULL),
105 fCutsEvent(NULL),
106 fCutsRP(NULL),
107 fCutsPOI(NULL),
108 fCutContainer(new TList()),
109 fQAList(NULL),
110 fMinMult(0),
111 fMaxMult(10000000),
112 fMinA(-1.0),
113 fMaxA(-0.01),
114 fMinB(0.01),
115 fMaxB(1.0),
116 fGenHeader(NULL),
117 fPythiaGenHeader(NULL),
118 fHijingGenHeader(NULL),
119 fFlowTrack(NULL),
120 fAnalysisUtil(NULL),
121 fQAon(kFALSE),
122 fLoadCandidates(kFALSE),
123 fNbinsMult(10000),
124 fNbinsPt(100),
125 fNbinsPhi(100),
126 fNbinsEta(200),
127 fNbinsQ(500),
128 fNbinsMass(1),
129 fMultMin(0.),
130 fMultMax(10000.),
131 fPtMin(0.),
132 fPtMax(10.),
133 fPhiMin(0.),
134 fPhiMax(TMath::TwoPi()),
135 fEtaMin(-5.),
136 fEtaMax(5.),
137 fQMin(0.),
138 fQMax(3.),
139 fMassMin(-1.),
140 fMassMax(0.),
141 fHistWeightvsPhiMin(0.),
142 fHistWeightvsPhiMax(3.),
143 fExcludedEtaMin(0.),
144 fExcludedEtaMax(0.),
145 fExcludedPhiMin(0.),
146 fExcludedPhiMax(0.),
147 fAfterburnerOn(kFALSE),
148 fNonFlowNumberOfTrackClones(0),
149 fV1(0.),
150 fV2(0.),
151 fV3(0.),
152 fV4(0.),
153 fV5(0.),
154 fDifferentialV2(0),
155 fFlowEvent(NULL),
156 fShuffleTracks(kFALSE),
157 fMyTRandom3(NULL),
158 fAnalysisInput(kAOD),
159 fIsMCInput(kFALSE),
160 fUseMCCen(kTRUE),
161 fRejectPileUp(kTRUE),
162 fRejectPileUpTight(kFALSE),
163 fCentrLowLim(0.),
164 fCentrUpLim(100.),
165 fCentrEstimator(kV0M),
166 fOutput(0x0),
167 fhZNCvsZNA(0x0),
168 fhZDCCvsZDCCA(0x0),
169 fhZNCvsZPC(0x0),
170 fhZNAvsZPA(0x0),
171 fhZNvsZP(0x0),
172 fhZNvsVZERO(0x0),
173 fhZDCvsVZERO(0x0),
174 fhZDCvsTracklets(0x0),
175 fhZDCvsNclu1(0x0),
176 fhDebunch(0x0),
177 fhZNCcentroid(0x0),
178 fhZNAcentroid(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 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 fUseZDCSpectraCorr(kFALSE),
200 fZDCSpectraCorrList(NULL),
201 fSpectraMCList(NULL),
202 fBadTowerStuffList(NULL),
203 fCachedRunNum(0),
204 fhZNSpectra(0x0),
205 fhZNSpectraCor(0x0),
206 fhZNSpectraPow(0x0),
207 fhZNBCCorr(0x0)
208 {
209  for(int i=0; i<5; i++){
210  fhZNCPM[i] = 0x0;
211  fhZNAPM[i] = 0x0;
212  }
213  for(int i=0; i<4; i++){
214  fhZNCPMQiPMC[i] = 0x0;
215  fhZNAPMQiPMC[i] = 0x0;
216  }
217  for(Int_t r=0; r<fCRCMaxnRun; r++) {
218  fRunList[r] = 0;
219  }
220  for(Int_t c=0; c<2; c++) {
221  for(Int_t i=0; i<5; i++) {
222  fTowerGainEq[c][i] = NULL;
223  }
224  }
225  for(Int_t c=0; c<100; c++) {
226  fBadTowerCalibHist[c] = NULL;
227  }
228  for(Int_t i=0; i<8; i++) {
229  SpecCorMu1[i] = NULL;
230  SpecCorMu2[i] = NULL;
231  FitSpecCorSi[i] = NULL;
232  SpecCorAv[i] = NULL;
233  }
234  this->InitializeRunArrays();
235  fMyTRandom3 = new TRandom3(1);
236  gRandom->SetSeed(fMyTRandom3->Integer(65539));
237  for(Int_t j=0; j<2; j++) {
238  for(Int_t c=0; c<10; c++) {
239  fPtSpecGen[j][c] = NULL;
240  fPtSpecFB32[j][c] = NULL;
241  fPtSpecFB96[j][c] = NULL;
242  fPtSpecFB128[j][c] = NULL;
243  fPtSpecFB768[j][c] = NULL;
244  }
245  }
246 }
247 
248 //________________________________________________________________________
249 AliAnalysisTaskCRCZDC::AliAnalysisTaskCRCZDC(const char *name, TString RPtype, Bool_t on, UInt_t iseed, Bool_t bCandidates):
250 AliAnalysisTaskSE(name),
251 fAnalysisType("AUTOMATIC"),
252 fRPType(RPtype),
253 fCFManager1(NULL),
254 fCFManager2(NULL),
255 fCutsEvent(NULL),
256 fCutsRP(NULL),
257 fCutsPOI(NULL),
258 fCutContainer(new TList()),
259 fQAList(NULL),
260 fMinMult(0),
261 fMaxMult(10000000),
262 fMinA(-1.0),
263 fMaxA(-0.01),
264 fMinB(0.01),
265 fMaxB(1.0),
266 fQAon(on),
267 fLoadCandidates(bCandidates),
268 fNbinsMult(10000),
269 fNbinsPt(100),
270 fNbinsPhi(100),
271 fNbinsEta(200),
272 fNbinsQ(500),
273 fNbinsMass(1),
274 fMultMin(0.),
275 fMultMax(10000.),
276 fPtMin(0.),
277 fPtMax(10.),
278 fPhiMin(0.),
279 fPhiMax(TMath::TwoPi()),
280 fEtaMin(-5.),
281 fEtaMax(5.),
282 fQMin(0.),
283 fQMax(3.),
284 fMassMin(-1.),
285 fMassMax(0.),
286 fHistWeightvsPhiMin(0.),
287 fHistWeightvsPhiMax(3.),
288 fExcludedEtaMin(0.),
289 fExcludedEtaMax(0.),
290 fExcludedPhiMin(0.),
291 fExcludedPhiMax(0.),
292 fAfterburnerOn(kFALSE),
293 fNonFlowNumberOfTrackClones(0),
294 fV1(0.),
295 fV2(0.),
296 fV3(0.),
297 fV4(0.),
298 fV5(0.),
299 fDifferentialV2(0),
300 fFlowEvent(NULL),
301 fShuffleTracks(kFALSE),
302 fMyTRandom3(NULL),
303 fAnalysisInput(kAOD),
304 fIsMCInput(kFALSE),
305 fUseMCCen(kTRUE),
306 fRejectPileUp(kTRUE),
307 fRejectPileUpTight(kFALSE),
308 fCentrLowLim(0.),
309 fCentrUpLim(100.),
310 fCentrEstimator(kV0M),
311 fOutput(0x0),
312 fhZNCvsZNA(0x0),
313 fhZDCCvsZDCCA(0x0),
314 fhZNCvsZPC(0x0),
315 fhZNAvsZPA(0x0),
316 fhZNvsZP(0x0),
317 fhZNvsVZERO(0x0),
318 fhZDCvsVZERO(0x0),
319 fhZDCvsTracklets(0x0),
320 fhZDCvsNclu1(0x0),
321 fhDebunch(0x0),
322 fhZNCcentroid(0x0),
323 fhZNAcentroid(0x0),
324 fhAsymm(0x0),
325 fhZNAvsAsymm(0x0),
326 fhZNCvsAsymm(0x0),
327 fhZNCvscentrality(0x0),
328 fhZNAvscentrality(0x0),
329 fDataSet(kAny),
330 fCRCnRun(0),
331 fZDCGainAlpha(0.395),
332 fGenHeader(NULL),
333 fPythiaGenHeader(NULL),
334 fHijingGenHeader(NULL),
335 fFlowTrack(NULL),
336 fAnalysisUtil(NULL),
337 fStack(0x0),
338 fCutTPC(kFALSE),
339 fCenDis(0x0),
340 fMultSelection(0x0),
341 fPileUpCount(0x0),
342 fPileUpMultSelCount(0x0),
343 fMultTOFLowCut(0x0),
344 fMultTOFHighCut(0x0),
345 fUseTowerEq(kFALSE),
346 fTowerEqList(NULL),
347 fUseBadTowerCalib(kFALSE),
348 fBadTowerCalibList(NULL),
349 fUseZDCSpectraCorr(kFALSE),
350 fZDCSpectraCorrList(NULL),
351 fSpectraMCList(NULL),
352 fBadTowerStuffList(NULL),
353 fCachedRunNum(0),
354 fhZNSpectra(0x0),
355 fhZNSpectraCor(0x0),
356 fhZNSpectraPow(0x0),
357 fhZNBCCorr(0x0)
358 {
359 
360  for(int i=0; i<5; i++){
361  fhZNCPM[i] = 0x0;
362  fhZNAPM[i] = 0x0;
363  }
364  for(int i=0; i<4; i++){
365  fhZNCPMQiPMC[i] = 0x0;
366  fhZNAPMQiPMC[i] = 0x0;
367  }
368  for(Int_t r=0; r<fCRCMaxnRun; r++) {
369  fRunList[r] = 0;
370  }
371  for(Int_t c=0; c<2; c++) {
372  for(Int_t i=0; i<5; i++) {
373  fTowerGainEq[c][i] = NULL;
374  }
375  }
376  for(Int_t c=0; c<100; c++) {
377  fBadTowerCalibHist[c] = NULL;
378  }
379  for(Int_t i=0; i<8; i++) {
380  SpecCorMu1[i] = NULL;
381  SpecCorMu2[i] = NULL;
382  FitSpecCorSi[i] = NULL;
383  SpecCorAv[i] = NULL;
384  }
385  this->InitializeRunArrays();
386  fMyTRandom3 = new TRandom3(iseed);
387  gRandom->SetSeed(fMyTRandom3->Integer(65539));
388 
389  DefineInput(0, TChain::Class());
390  // Define output slots here
391  // Define here the flow event output
392  DefineOutput(1, AliFlowEventSimple::Class());
393  DefineOutput(2, TList::Class());
394 
395  for(Int_t j=0; j<2; j++) {
396  for(Int_t c=0; c<10; c++) {
397  fPtSpecGen[j][c] = NULL;
398  fPtSpecFB32[j][c] = NULL;
399  fPtSpecFB96[j][c] = NULL;
400  fPtSpecFB128[j][c] = NULL;
401  fPtSpecFB768[j][c] = NULL;
402  }
403  }
404 
405 }
406 
407 //________________________________________________________________________
409 {
410  // Destructor
411  if(fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()){
412  delete fOutput; fOutput=0;
413  }
414  delete fMyTRandom3;
415  delete fFlowEvent;
416  delete fFlowTrack;
417  delete fCutsEvent;
418  if (fTowerEqList) delete fTowerEqList;
421  if (fAnalysisUtil) delete fAnalysisUtil;
422  if (fQAList) delete fQAList;
423  if (fCutContainer) fCutContainer->Delete(); delete fCutContainer;
424 }
425 
426 //________________________________________________________________________
428 {
429  for(Int_t r=0;r<fCRCMaxnRun;r++) {
430  fCRCQVecListRun[r] = NULL;
431  for(Int_t k=0;k<fCRCnTow;k++) {
432  fZNCTower[r][k] = NULL;
433  fZNATower[r][k] = NULL;
434  }
435  }
436  // for(Int_t i=0;i<fnCen;i++) {
437  // fPtPhiEtaRbRFB128[r][i] = NULL;
438  // fPtPhiEtaRbRFB768[r][i] = NULL;
439  // }
440  // }
441 }
442 
443 //________________________________________________________________________
445 {
446  // Create the output containers
447 
448  if (!(fAnalysisType == "AOD" || fAnalysisType == "MCkine" || fAnalysisType == "MCAOD" || fAnalysisType == "AUTOMATIC" || fAnalysisType == "MCESD"))
449  {
450  AliError("WRONG ANALYSIS TYPE! only MCESD, MCkine, MCAOD, AOD and AUTOMATIC are allowed.");
451  exit(1);
452  }
453 
454  //set the common constants
457  cc->SetNbinsPt(fNbinsPt);
458  cc->SetNbinsPhi(fNbinsPhi);
459  cc->SetNbinsEta(fNbinsEta);
460  cc->SetNbinsQ(fNbinsQ);
462  cc->SetMultMin(fMultMin);
463  cc->SetMultMax(fMultMax);
464  cc->SetPtMin(fPtMin);
465  cc->SetPtMax(fPtMax);
466  cc->SetPhiMin(fPhiMin);
467  cc->SetPhiMax(fPhiMax);
468  cc->SetEtaMin(fEtaMin);
469  cc->SetEtaMax(fEtaMax);
470  cc->SetQMin(fQMin);
471  cc->SetQMax(fQMax);
472  cc->SetMassMin(fMassMin);
473  cc->SetMassMax(fMassMax);
476 
477  fFlowEvent = new AliFlowEvent(20000);
478  fFlowTrack = new AliFlowTrack();
479 
480  //printf(" AliAnalysisTaskCRCZDC::UserCreateOutputObjects()\n\n");
481  fOutput = new TList();
482  fOutput->SetOwner(kTRUE);
483  //fOutput->SetName("output");
484 
485  if (fQAon) {
486  fQAList = new TList();
487  fQAList->SetOwner(kTRUE);
488  fQAList->SetName("AliFlowEventCuts QA");
489  if (fCutsEvent->GetQA()) fQAList->Add(fCutsEvent->GetQA()); //0
490  if (fCutsRP->GetQA()) fQAList->Add(fCutsRP->GetQA()); //1
491  if (fCutsPOI->GetQA())fQAList->Add(fCutsPOI->GetQA()); //2
492  fOutput->Add(fQAList);
493  }
494 
495  fCenDis = new TH1F("fCenDis", "fCenDis", 100, 0., 100.);
496  fOutput->Add(fCenDis);
497  fPileUpCount = new TH1F("fPileUpCount", "fPileUpCount", 9, 0., 9.);
498  fPileUpCount->GetXaxis()->SetBinLabel(1,"plpMV");
499  fPileUpCount->GetXaxis()->SetBinLabel(2,"fromSPD");
500  fPileUpCount->GetXaxis()->SetBinLabel(3,"RefMultiplicityComb08");
501  fPileUpCount->GetXaxis()->SetBinLabel(4,"IncompleteDAQ");
502  fPileUpCount->GetXaxis()->SetBinLabel(5,"abs(V0M-CL1)>7.5");
503  fPileUpCount->GetXaxis()->SetBinLabel(6,"missingVtx");
504  fPileUpCount->GetXaxis()->SetBinLabel(7,"inconsistentVtx");
505  fPileUpCount->GetXaxis()->SetBinLabel(8,"multESDTPCDif");
506  fPileUpCount->GetXaxis()->SetBinLabel(9,"extraPileUpMultSel");
507  fOutput->Add(fPileUpCount);
508  fPileUpMultSelCount = new TH1F("fPileUpMultSelCount", "fPileUpMultSelCount", 8, 0., 8.);
509  fPileUpMultSelCount->GetXaxis()->SetBinLabel(1,"IsNotPileup");
510  fPileUpMultSelCount->GetXaxis()->SetBinLabel(2,"IsNotPileupMV");
511  fPileUpMultSelCount->GetXaxis()->SetBinLabel(3,"IsNotPileupInMultBins");
512  fPileUpMultSelCount->GetXaxis()->SetBinLabel(4,"InconsistentVertices");
513  fPileUpMultSelCount->GetXaxis()->SetBinLabel(5,"TrackletVsCluster");
514  fPileUpMultSelCount->GetXaxis()->SetBinLabel(6,"AsymmetricInVZERO");
515  fPileUpMultSelCount->GetXaxis()->SetBinLabel(7,"IncompleteDAQ");
516  fPileUpMultSelCount->GetXaxis()->SetBinLabel(8,"GoodVertex2016");
518 
519  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);
520  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);
521  fOutput->Add(fMultTOFLowCut);
522  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);
523  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);
524  fOutput->Add(fMultTOFHighCut);
525 
526  for(Int_t c=0; c<2; c++) {
527  for(Int_t i=0; i<5; i++) {
528  fTowerGainEq[c][i] = new TH1D();
529  fOutput->Add(fTowerGainEq[c][i]);
530  }
531  }
532 
533  fBadTowerStuffList = new TList();
534  fBadTowerStuffList->SetOwner(kTRUE);
535  fBadTowerStuffList->SetName("BadTowerCalib");
537 
538  if(fBadTowerCalibList) {
539  for(Int_t c=0; c<100; c++) {
540  fBadTowerCalibHist[c] = (TH2D*)fBadTowerCalibList->FindObject(Form("TH2Resp[%d]",c));
542  }
543  }
544  if(fZDCSpectraCorrList) {
545  for(Int_t i=0; i<8; i++) {
546  SpecCorMu1[i] = (TH1D*)fZDCSpectraCorrList->FindObject(Form("SpecCorMu1[%d]",i));
547  fOutput->Add(SpecCorMu1[i]);
548  SpecCorMu2[i] = (TH1D*)fZDCSpectraCorrList->FindObject(Form("SpecCorMu2[%d]",i));
549  fOutput->Add(SpecCorMu2[i]);
550  SpecCorAv[i] = (TH1D*)fZDCSpectraCorrList->FindObject(Form("SpecCorAv[%d]",i));
551  fOutput->Add(SpecCorAv[i]);
552  }
553  for(Int_t i=0; i<8; i++) {
554  TH1D* SpecCorMu = (TH1D*)fZDCSpectraCorrList->FindObject(Form("SpecCorSi[%d]",i));
555  TF1* FisrtFit = new TF1("FirstFit","pol5",0.,100.);
556  SpecCorMu->Fit(FisrtFit,"QRN","",10.,80.);
557  for (Int_t bx=1; bx<=SpecCorMu->GetNbinsX(); bx++) {
558  if(SpecCorMu->GetXaxis()->GetBinCenter(bx)>5. && SpecCorMu->GetXaxis()->GetBinCenter(bx)<75.) {
559  if(fabs(SpecCorMu->GetBinContent(bx)-FisrtFit->Eval(SpecCorMu->GetXaxis()->GetBinCenter(bx))) > SpecCorMu->GetBinContent(bx)*0.01) SpecCorMu->SetBinContent(bx,0.);
560  }
561  }
562  for (Int_t bx=1; bx<=SpecCorMu->GetNbinsX(); bx++) {
563  if(SpecCorMu->GetBinContent(bx)==0. && SpecCorMu->GetXaxis()->GetBinCenter(bx)<79.) {
564  Double_t xmin=0.,xmax=0.;
565  Int_t cmin=0,cmax=0;
566  while(xmin==0.) {
567  cmin++;
568  xmin = SpecCorMu->GetBinContent(bx-cmin);
569  }
570  while(xmax==0.) {
571  cmax++;
572  xmax = SpecCorMu->GetBinContent(bx+cmax);
573  }
574  SpecCorMu->SetBinContent(bx,(xmin+xmax)/2.);
575  }
576  if(SpecCorMu->GetXaxis()->GetBinCenter(bx)>79.) {
577  SpecCorMu->SetBinContent(bx,SpecCorMu->GetBinContent(SpecCorMu->GetXaxis()->FindBin(78.5)));
578  }
579  }
580  FitSpecCorSi[i] = new TF1(Form("FitSpecCorSi[%d]",i),"pol7",0.,100.);
581  SpecCorMu->Fit(FitSpecCorSi[i],"QRN","",0.,80.);
582  fOutput->Add(FitSpecCorSi[i]);
583  }
584  }
585 
586  fhZNSpectra = new TH3D("fhZNSpectra","fhZNSpectra",100,0.,100.,8,0.,8.,1000,0.,1.E5);
587  fOutput->Add(fhZNSpectra);
588  fhZNSpectraCor = new TH3D("fhZNSpectraCor","fhZNSpectraCor",100,0.,100.,8,0.,8.,1000,0.,1.E5);
589  fOutput->Add(fhZNSpectraCor);
590  fhZNSpectraPow = new TH3D("fhZNSpectraPow","fhZNSpectraPow",100,0.,100.,8,0.,8.,1000,0.,TMath::Power(1.E5,fZDCGainAlpha));
591  fOutput->Add(fhZNSpectraPow);
592  fhZNBCCorr = new TH3D("fhZNBCCorr","fhZNBCCorr",100,0.,100.,500,0.,1.E5,500,0.,1.E5);
593  fOutput->Add(fhZNBCCorr);
594 
595  if(fAnalysisType == "MCAOD") {
596 
597  fSpectraMCList = new TList();
598  fSpectraMCList->SetOwner(kTRUE);
599  fSpectraMCList->SetName("Spectra");
600  fOutput->Add(fSpectraMCList);
601 
602  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.};
603  for(Int_t j=0; j<2; j++) {
604  for(Int_t c=0; c<10; c++) {
605  fPtSpecGen[j][c] = new TH1F(Form("fPtSpecGen[%d][%d]",j,c), Form("fPtSpecGen[%d][%d]",j,c), 23, xmin);
606  fSpectraMCList->Add(fPtSpecGen[j][c]);
607  fPtSpecFB32[j][c] = new TH1F(Form("fPtSpecFB32[%d][%d]",j,c), Form("fPtSpecFB32[%d][%d]",j,c), 23, xmin);
608  fSpectraMCList->Add(fPtSpecFB32[j][c]);
609  fPtSpecFB96[j][c] = new TH1F(Form("fPtSpecFB96[%d][%d]",j,c), Form("fPtSpecFB96[%d][%d]",j,c), 23, xmin);
610  fSpectraMCList->Add(fPtSpecFB96[j][c]);
611  fPtSpecFB128[j][c] = new TH1F(Form("fPtSpecFB128[%d][%d]",j,c), Form("fPtSpecFB128[%d][%d]",j,c), 23, xmin);
612  fSpectraMCList->Add(fPtSpecFB128[j][c]);
613  fPtSpecFB768[j][c] = new TH1F(Form("fPtSpecFB768[%d][%d]",j,c), Form("fPtSpecFB768[%d][%d]",j,c), 23, xmin);
614  fSpectraMCList->Add(fPtSpecFB768[j][c]);
615  }
616  }
617  }
618 
619  fAnalysisUtil = new AliAnalysisUtils();
620  fAnalysisUtil->SetUseMVPlpSelection(kTRUE);
621  fAnalysisUtil->SetUseOutOfBunchPileUp(kTRUE);
622 
623  for(int i=0; i<5; i++){
624  char hname[20];
625  sprintf(hname,"hZNCPM%d",i);
626  fhZNCPM[i] = new TH1F(hname, hname, 200, -50., 140000);
627  fOutput->Add(fhZNCPM[i]);
628  //
629  sprintf(hname,"hZNAPM%d",i);
630  fhZNAPM[i] = new TH1F(hname, hname, 200, -50., 140000);
631  fOutput->Add(fhZNAPM[i]);
632  //
633  if(i<4){
634  //
635  char hnamenc[20];
636  sprintf(hnamenc, "hZNCPMQ%dPMC",i+1);
637  fhZNCPMQiPMC[i] = new TH1F(hnamenc, hnamenc, 100, 0., 1.);
638  fOutput->Add(fhZNCPMQiPMC[i]);
639  //
640  char hnamena[20];
641  sprintf(hnamena, "hZNAPMQ%dPMC",i+1);
642  fhZNAPMQiPMC[i] = new TH1F(hnamena, hnamena, 100, 0., 1.);
643  fOutput->Add(fhZNAPMQiPMC[i]);
644  }
645  }
646 
647  fhZNCvsZNA = new TH2F("hZNCvsZNA","hZNCvsZNA",200,-50.,140000,200,-50.,140000);
648  fOutput->Add(fhZNCvsZNA);
649  fhZDCCvsZDCCA = new TH2F("hZDCCvsZDCCA","hZDCCvsZDCCA",200,0.,180000.,200,0.,200000.);
650  fOutput->Add(fhZDCCvsZDCCA);
651  fhZNCvsZPC = new TH2F("hZNCvsZPC","hZNCvsZPC",200,-50.,50000,200,-50.,140000);
652  fOutput->Add(fhZNCvsZPC);
653  fhZNAvsZPA = new TH2F("hZNAvsZPA","hZNAvsZPA",200,-50.,50000,200,-50.,140000);
654  fOutput->Add(fhZNAvsZPA);
655  fhZNvsZP = new TH2F("hZNvsZP","hZNvsZP",200,-50.,80000,200,-50.,200000);
656  fOutput->Add(fhZNvsZP);
657  fhZNvsVZERO = new TH2F("hZNvsVZERO","hZNvsVZERO",250,0.,25000.,200,0.,200000.);
658  fOutput->Add(fhZNvsVZERO);
659  fhZDCvsVZERO = new TH2F("hZDCvsVZERO","hZDCvsVZERO",250,0.,25000.,250,0.,250000.);
660  fOutput->Add(fhZDCvsVZERO);
661  fhZDCvsTracklets = new TH2F("hZDCvsTracklets","hZDCvsTracklets",200,0.,4000.,250,0.,250000.);
663  fhZDCvsNclu1 = new TH2F("hZDCvsNclu1", "hZDCvsNclu1", 200, 0.,8000.,200,0.,250000.);
664  fOutput->Add(fhZDCvsNclu1);
665  fhDebunch = new TH2F("hDebunch","hDebunch",240,-100.,-40.,240,-30.,30.);
666  fOutput->Add(fhDebunch);
667  fhZNCcentroid = new TH2F("hZNCcentroid","hZNCcentroid",100,-3.5,3.5,100,-3.5,3.5);
668  fOutput->Add(fhZNCcentroid);
669  fhZNAcentroid = new TH2F("hZNAcentroid","hZNAcentroid",100,-3.5,3.5,100,-3.5,3.5);
670  fOutput->Add(fhZNAcentroid);
671 
672  fhAsymm = new TH1F("hAsymm" , "Asimmetry ",200,-1.,1.);
673  fOutput->Add(fhAsymm);
674  fhZNAvsAsymm = new TH2F("hZNAvsAsymm","ZNA vs. asymm.",200,-1.,1.,200,0.,80.);
675  fOutput->Add(fhZNAvsAsymm);
676  fhZNCvsAsymm = new TH2F("hZNCvsAsymm","ZNC vs. asymm.",200,-1.,1.,200,0.,80.);
677  fOutput->Add(fhZNCvsAsymm);
678 
679  fhZNCvscentrality = new TH2F("hZNCvscentrality","ZNC vs. centrality",100,0.,100.,100,0.,100.);
681  fhZNAvscentrality = new TH2F("hZNAvscentrality","ZNA vs. centrality",100,0.,100.,100,0.,100.);
683 
684  //********************************************************************
685 
686  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};
687 
688  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};
689 
690  // 12 low IR: 244917, 244918, 244975, 244980, 244982, 244983, 245064, 245066, 245068, 246390, 246391, 246392
691  // 80 high IR ("CentralBarrelTracking" good runs): 246994, 246991, 246989, 246984, 246982, 246980, 246948, 246945, 246928, 246871, 246870, 246867, 246865, 246864, 246859, 246858, 246851, 246847, 246846, 246845, 246844, 246810, 246809, 246808, 246807, 246805, 246804, 246766, 246765, 246763, 246760, 246759, 246758, 246757, 246751, 246750, 246676, 246675, 246540, 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
692 
693  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};
694 
695  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};
696 
697  if(fDataSet==k2010) {fCRCnRun=92;}
698  if(fDataSet==k2011) {fCRCnRun=119;}
699  if(fDataSet==k2015) {fCRCnRun=90;}
700  if(fDataSet==k2015v6) {fCRCnRun=91;}
701  if(fDataSet==kAny) {fCRCnRun=1;}
702 
703  Int_t d=0;
704  for(Int_t r=0; r<fCRCnRun; r++) {
705  if(fDataSet==k2010) fRunList[d] = dRun10h[r];
706  if(fDataSet==k2011) fRunList[d] = dRun11h[r];
707  if(fDataSet==k2015) fRunList[d] = dRun15o[r];
708  if(fDataSet==k2015v6) fRunList[d] = dRun15ov6[r];
709  if(fDataSet==kAny) fRunList[d] = 1;
710  d++;
711  }
712 
713  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.};
714  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.};
715  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};
716 
717  for(Int_t r=0;r<fCRCnRun;r++) {
718  fCRCQVecListRun[r] = new TList();
719  fCRCQVecListRun[r]->SetName(Form("Run %d",fRunList[r]));
720  fCRCQVecListRun[r]->SetOwner(kTRUE);
721  fOutput->Add(fCRCQVecListRun[r]);
722 
723  if(!fUseTowerEq) {
724  for(Int_t k=0;k<fCRCnTow;k++) {
725  fZNCTower[r][k] = new TProfile(Form("fZNCTower[%d][%d]",fRunList[r],k),Form("fZNCTower[%d][%d]",fRunList[r],k),100,0.,100.,"s");
726  fZNCTower[r][k]->Sumw2();
727  fCRCQVecListRun[r]->Add(fZNCTower[r][k]);
728  fZNATower[r][k] = new TProfile(Form("fZNATower[%d][%d]",fRunList[r],k),Form("fZNATower[%d][%d]",fRunList[r],k),100,0.,100.,"s");
729  fZNATower[r][k]->Sumw2();
730  fCRCQVecListRun[r]->Add(fZNATower[r][k]);
731  }
732  }
733 
734  // for(Int_t i=0;i<fnCen;i++) {
735  // fPtPhiEtaRbRFB128[r][i] = new TH3F(Form("fPtPhiEtaRbRFB128[%d][%d]",r,i),Form("fPtPhiEtaRbRFB128[%d][%d]",r,i),14, ptmin, 16, phimin, 16, etamin);
736  // fCRCQVecListRun[r]->Add(fPtPhiEtaRbRFB128[r][i]);
737  // fPtPhiEtaRbRFB768[r][i] = new TH3F(Form("fPtPhiEtaRbRFB768[%d][%d]",r,i),Form("fPtPhiEtaRbRFB768[%d][%d]",r,i),14, ptmin, 16, phimin, 16, etamin);
738  // fCRCQVecListRun[r]->Add(fPtPhiEtaRbRFB768[r][i]);
739  // }
740  }
741 
742  PostData(2, fOutput);
743 }
744 
745 //________________________________________________________________________
747 {
748  // Execute analysis for current event:
749  AliMCEvent* McEvent = MCEvent();
750  AliAODEvent *aod = dynamic_cast<AliAODEvent*>(InputEvent());
751  AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
752  // AliMultiplicity* myTracklets = NULL;
753  // AliESDPmdTrack* pmdtracks = NULL;
754  // int availableINslot=1;
755 
756  if (!(fCutsRP&&fCutsPOI&&fCutsEvent)) {
757  AliError("cuts not set");
758  return;
759  }
760 
761  Int_t RunBin=-1, bin=0;
762  Int_t RunNum = aod->GetRunNumber();
763  for(Int_t c=0;c<fCRCnRun;c++) {
764  if(fRunList[c]==RunNum) RunBin=bin;
765  else bin++;
766  }
767  if(RunBin==-1) return;
768  if(fDataSet==kAny) RunBin=0;
769 
770  //DEFAULT - automatically takes care of everything
771  if (fAnalysisType == "AUTOMATIC") {
772 
773  // get centrality
774  Double_t centrV0M=300, centrCL1=300, centrCL0=300, centrTRK=300;
775  if(fDataSet!=k2015 && fDataSet!=k2015v6) {
776  centrV0M = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("V0M");
777  centrCL1 = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("CL1");
778  centrCL0 = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("CL0");
779  centrTRK = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("TRK");
780  } else {
781  fMultSelection = (AliMultSelection*) InputEvent()->FindListObject("MultSelection");
782  if(!fMultSelection) {
783  //If you get this warning (and lPercentiles 300) please check that the AliMultSelectionTask actually ran (before your task)
784  AliWarning("AliMultSelection object not found!");
785  } else {
786  centrV0M = fMultSelection->GetMultiplicityPercentile("V0M");
787  centrCL1 = fMultSelection->GetMultiplicityPercentile("CL1");
788  centrCL0 = fMultSelection->GetMultiplicityPercentile("CL0");
789  centrTRK = fMultSelection->GetMultiplicityPercentile("TRK");
790  }
791  }
792 
793  //check event cuts
794  if (InputEvent()) {
795  if(!fCutsEvent->IsSelected(InputEvent(),MCEvent())) return;
796  if(fRejectPileUp) {
797  if(fDataSet!=k2015 && fDataSet!=k2015v6) {
798 
799  Bool_t BisPileup=kFALSE;
800 
801  // check anyway pileup
802  if (plpMV(aod)) {
803  fPileUpCount->Fill(0.5);
804  BisPileup=kTRUE;
805  }
806 
807  Short_t isPileup = aod->IsPileupFromSPD(3);
808  if (isPileup != 0) {
809  fPileUpCount->Fill(1.5);
810  BisPileup=kTRUE;
811  }
812 
813  if (((AliAODHeader*)aod->GetHeader())->GetRefMultiplicityComb08() < 0) {
814  fPileUpCount->Fill(2.5);
815  BisPileup=kTRUE;
816  }
817 
818  if (aod->IsIncompleteDAQ()) {
819  fPileUpCount->Fill(3.5);
820  BisPileup=kTRUE;
821  }
822 
823  // check vertex consistency
824  const AliAODVertex* vtTrc = aod->GetPrimaryVertex();
825  const AliAODVertex* vtSPD = aod->GetPrimaryVertexSPD();
826 
827  if (vtTrc->GetNContributors() < 2 || vtSPD->GetNContributors()<1) {
828  fPileUpCount->Fill(5.5);
829  BisPileup=kTRUE;
830  }
831 
832  double covTrc[6], covSPD[6];
833  vtTrc->GetCovarianceMatrix(covTrc);
834  vtSPD->GetCovarianceMatrix(covSPD);
835 
836  double dz = vtTrc->GetZ() - vtSPD->GetZ();
837 
838  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
839  double errTrc = TMath::Sqrt(covTrc[5]);
840  double nsigTot = dz/errTot;
841  double nsigTrc = dz/errTrc;
842 
843  if (TMath::Abs(dz)>0.2 || TMath::Abs(nsigTot)>10 || TMath::Abs(nsigTrc)>20) {
844  fPileUpCount->Fill(6.5);
845  BisPileup=kTRUE;
846  }
847 
848  if (fAnalysisUtil->IsPileUpEvent(InputEvent())) {
849  fPileUpCount->Fill(7.5);
850  BisPileup=kTRUE;
851  }
852 
853 // if(BisPileup) return;
854  } else {
855  // pileup from AliMultSelection
856  if(!fMultSelection->GetThisEventIsNotPileup()) fPileUpMultSelCount->Fill(0.5);
857  if(!fMultSelection->GetThisEventIsNotPileupMV()) fPileUpMultSelCount->Fill(1.5);
858  if(!fMultSelection->GetThisEventIsNotPileupInMultBins()) fPileUpMultSelCount->Fill(2.5);
859  if(!fMultSelection->GetThisEventHasNoInconsistentVertices()) fPileUpMultSelCount->Fill(3.5);
860  if(!fMultSelection->GetThisEventPassesTrackletVsCluster()) fPileUpMultSelCount->Fill(4.5);
861  if(!fMultSelection->GetThisEventIsNotAsymmetricInVZERO()) fPileUpMultSelCount->Fill(5.5);
862  if(!fMultSelection->GetThisEventIsNotIncompleteDAQ()) fPileUpMultSelCount->Fill(6.5);
863  if(!fMultSelection->GetThisEventHasGoodVertex2016()) fPileUpMultSelCount->Fill(7.5);
864 
865  Bool_t BisPileup=kFALSE;
866 
867  // pile-up a la Dobrin for LHC15o
868  if (plpMV(aod)) {
869  fPileUpCount->Fill(0.5);
870  BisPileup=kTRUE;
871  }
872 
873  Short_t isPileup = aod->IsPileupFromSPD(3);
874  if (isPileup != 0) {
875  fPileUpCount->Fill(1.5);
876  BisPileup=kTRUE;
877  }
878 
879  if (((AliAODHeader*)aod->GetHeader())->GetRefMultiplicityComb08() < 0) {
880  fPileUpCount->Fill(2.5);
881  BisPileup=kTRUE;
882  }
883 
884  if (aod->IsIncompleteDAQ()) {
885  fPileUpCount->Fill(3.5);
886  BisPileup=kTRUE;
887  }
888 
889  if(fabs(centrV0M-centrCL1)>7.5) {
890  fPileUpCount->Fill(4.5);
891  BisPileup=kTRUE;
892  }
893 
894  // check vertex consistency
895  const AliAODVertex* vtTrc = aod->GetPrimaryVertex();
896  const AliAODVertex* vtSPD = aod->GetPrimaryVertexSPD();
897 
898  if (vtTrc->GetNContributors() < 2 || vtSPD->GetNContributors()<1) {
899  fPileUpCount->Fill(5.5);
900  BisPileup=kTRUE;
901  }
902 
903  double covTrc[6], covSPD[6];
904  vtTrc->GetCovarianceMatrix(covTrc);
905  vtSPD->GetCovarianceMatrix(covSPD);
906 
907  double dz = vtTrc->GetZ() - vtSPD->GetZ();
908 
909  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
910  double errTrc = TMath::Sqrt(covTrc[5]);
911  double nsigTot = dz/errTot;
912  double nsigTrc = dz/errTrc;
913 
914  if (TMath::Abs(dz)>0.2 || TMath::Abs(nsigTot)>10 || TMath::Abs(nsigTrc)>20) {
915  fPileUpCount->Fill(6.5);
916  BisPileup=kTRUE;
917  }
918 
919  // cuts on tracks
920  const Int_t nTracks = aod->GetNumberOfTracks();
921  Int_t multEsd = ((AliAODHeader*)aod->GetHeader())->GetNumberOfESDTracks();
922 
923  Int_t multTrk = 0;
924  Int_t multTrkBefC = 0;
925  Int_t multTrkTOFBefC = 0;
926  Int_t multTPC = 0;
927 
928  for (Int_t it = 0; it < nTracks; it++) {
929 
930  AliAODTrack* aodTrk = (AliAODTrack*)aod->GetTrack(it);
931  if (!aodTrk){
932  delete aodTrk;
933  continue;
934  }
935 
936 // if (aodTrk->TestFilterBit(32)){
937 // multTrkBefC++;
938 //
939 // if ( TMath::Abs(aodTrk->GetTOFsignalDz()) <= 10. && aodTrk->GetTOFsignal() >= 12000. && aodTrk->GetTOFsignal() <= 25000.)
940 // multTrkTOFBefC++;
941 //
942 // if ((TMath::Abs(aodTrk->Eta()) < 0.8) && (aodTrk->GetTPCNcls() >= 70) && (aodTrk->Pt() >= 0.2) && (aodTrk->Pt() < 20.))
943 // multTrk++;
944 // }
945 
946  if (aodTrk->TestFilterBit(128))
947  multTPC++;
948 
949  } // end of for (Int_t it = 0; it < nTracks; it++)
950 
951  Double_t multTPCn = multTPC;
952  Double_t multEsdn = multEsd;
953  Double_t multESDTPCDif = multEsdn - multTPCn*3.38;
954 
955  if (multESDTPCDif > (fRejectPileUpTight?700.:15000.)) {
956  fPileUpCount->Fill(7.5);
957  BisPileup=kTRUE;
958  }
959 
960  if(fRejectPileUpTight) {
961  if(BisPileup==kFALSE) {
962  if(!fMultSelection->GetThisEventIsNotPileup()) BisPileup=kTRUE;
963  if(!fMultSelection->GetThisEventIsNotPileupMV()) BisPileup=kTRUE;
964  if(!fMultSelection->GetThisEventIsNotPileupInMultBins()) BisPileup=kTRUE;
965  if(!fMultSelection->GetThisEventHasNoInconsistentVertices()) BisPileup=kTRUE;
966  if(!fMultSelection->GetThisEventPassesTrackletVsCluster()) BisPileup=kTRUE;
967  if(!fMultSelection->GetThisEventIsNotIncompleteDAQ()) BisPileup=kTRUE;
968  if(!fMultSelection->GetThisEventHasGoodVertex2016()) BisPileup=kTRUE;
969  if(BisPileup) fPileUpCount->Fill(8.5);
970  }
971  }
972 
973  if(BisPileup) return;
974  }
975  }
976  }
977 
978  //first attach all possible information to the cuts
979  fCutsRP->SetEvent( InputEvent(), MCEvent() ); //attach event
980  fCutsPOI->SetEvent( InputEvent(), MCEvent() );
981 
982  //then make the event
984 
986 
991  fFlowEvent->SetCentralityCL1(centrCL1);
992  fFlowEvent->SetCentralityTRK(centrTRK);
993  // fFlowEvent->SetNITSCL1(((AliVAODHeader*)aod->GetHeader())->GetNumberOfITSClusters(1));
994  Double_t SumV0=0.;
995  for(Int_t i=0; i<64; i++) {
996  if(std::isfinite(aod->GetVZEROEqMultiplicity(i))) SumV0 += aod->GetVZEROEqMultiplicity(i);
997  }
998  fFlowEvent->SetNITSCL1(SumV0);
999 
1000  Double_t vtxpos[3]={0.,0.,0.};
1001  vtxpos[0] = ((AliAODVertex*)aod->GetPrimaryVertex())->GetX();
1002  vtxpos[1] = ((AliAODVertex*)aod->GetPrimaryVertex())->GetY();
1003  vtxpos[2] = ((AliAODVertex*)aod->GetPrimaryVertex())->GetZ();
1004  fFlowEvent->SetVertexPosition(vtxpos);
1005 
1006  if (McEvent && McEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(McEvent);
1007 
1008  // run-by-run QA
1009  // for(Int_t jTracks = 0; jTracks<aod->GetNumberOfTracks(); jTracks++){
1010  // AliAODTrack* track = (AliAODTrack*)aod->GetTrack(jTracks);
1011  // if(!track) continue;
1012  // // general kinematic & quality cuts
1013  // if (track->Pt() < .2 || track->Pt() > 20. || TMath::Abs(track->Eta()) > .8 || track->GetTPCNcls() < 70) continue;
1014  // if (track->TestFilterBit(128)) fPtPhiEtaRbRFB128[RunBin][CenBin]->Fill(track->Pt(),track->Phi(),track->Eta());
1015  // if (track->TestFilterBit(768)) fPtPhiEtaRbRFB768[RunBin][CenBin]->Fill(track->Pt(),track->Phi(),track->Eta());
1016  // }
1017  fCenDis->Fill(centrV0M);
1018 
1019  }
1020 
1021  if (fAnalysisType == "MCAOD") {
1022 
1023  //check event cuts
1024  if (InputEvent()) {
1025  if(!fCutsEvent->IsSelected(InputEvent(),MCEvent())) return;
1026  if(fRejectPileUp && fAnalysisUtil->IsPileUpEvent(InputEvent())) return;
1027  }
1028 
1029  fFlowEvent->ClearFast();
1030 
1031  if(!McEvent) {
1032  AliError("ERROR: Could not retrieve MCEvent");
1033  return;
1034  }
1035  fStack = (TClonesArray*)aod->FindListObject(AliAODMCParticle::StdBranchName());
1036  if(!fStack){
1037  AliError("ERROR: Could not retrieve MCStack");
1038  return;
1039  }
1040 
1041  // get centrality (from AliMultSelection or AliCentrality)
1042  Double_t centr = 300;
1043  if(fDataSet==k2015 || fDataSet==k2015v6) {
1044  fMultSelection = (AliMultSelection*)aod->FindListObject("MultSelection");
1045  if(!fMultSelection) {
1046  //If you get this warning (and lPercentiles 300) please check that the AliMultSelectionTask actually ran (before your task)
1047  AliWarning("AliMultSelection object not found!");
1048  } else {
1049  centr = fMultSelection->GetMultiplicityPercentile("V0M");
1050  }
1051  } else {
1052  centr = (((AliVAODHeader*)aod->GetHeader())->GetCentralityP())->GetCentralityPercentile("V0M");
1053  }
1054  // centrality bin
1055  if (centr<fCentrLowLim || centr>=fCentrUpLim ) return;
1056  Int_t CenBin = -1;
1057  CenBin = GetCenBin(centr);
1058  if(CenBin==-1) return;
1059 
1060  // reconstructed
1061  for(Int_t jTracks = 0; jTracks<aod->GetNumberOfTracks(); jTracks++){
1062 
1063  AliAODTrack* track = (AliAODTrack*)aod->GetTrack(jTracks);
1064  if(!track) continue;
1065 
1066  // to select primaries
1067  Int_t lp = TMath::Abs(track->GetLabel());
1068 
1069  // general kinematic & quality cuts
1070  if (track->Pt() < .2 || track->Pt() > 20. || TMath::Abs(track->Eta()) > .8 || track->GetTPCNcls() < 70) continue;
1071 
1072  // Double_t DCAxy = track->DCA();
1073  // Double_t DCAz = track->ZAtDCA();
1074  // if (std::abs((Int_t)DCAxy)==999 || std::abs((Int_t)DCAz)==999) {
1075  // // re-evaluate the dca as it seems to not be natively present
1076  // // allowed only for tracks inside the beam pipe
1077  // Double_t pos[3] = {-99., -99., -99.};
1078  // track->GetPosition(pos);
1079  // if(pos[0]*pos[0]+pos[1]*pos[1] <= 3.*3.) {
1080  // AliAODTrack copy(*track); // stack copy
1081  // Double_t b[2] = {-99., -99.};
1082  // Double_t bCov[3] = {-99., -99., -99.};
1083  // if(copy.PropagateToDCA(aod->GetPrimaryVertex(), aod->GetMagneticField(), 100., b, bCov)) {
1084  // DCAxy = b[0];
1085  // DCAz = b[1];
1086  // }
1087  // }
1088  // }
1089 
1090  // test filter bits
1091  if (((AliAODMCParticle*)fStack->At(lp))->IsPhysicalPrimary()) {
1092  if (track->TestFilterBit(32)) fPtSpecFB32[0][CenBin]->Fill(track->Pt());
1093  if (track->TestFilterBit(96)) fPtSpecFB96[0][CenBin]->Fill(track->Pt());
1094  if (track->TestFilterBit(128)) fPtSpecFB128[0][CenBin]->Fill(track->Pt());
1095  if (track->TestFilterBit(768)) fPtSpecFB768[0][CenBin]->Fill(track->Pt());
1096  } else {
1097  if (track->TestFilterBit(32)) fPtSpecFB32[1][CenBin]->Fill(track->Pt());
1098  if (track->TestFilterBit(96)) fPtSpecFB96[1][CenBin]->Fill(track->Pt());
1099  if (track->TestFilterBit(128)) fPtSpecFB128[1][CenBin]->Fill(track->Pt());
1100  if (track->TestFilterBit(768)) fPtSpecFB768[1][CenBin]->Fill(track->Pt());
1101  }
1102 
1103  fCenDis->Fill(centr);
1104  }
1105 
1106  // generated (physical primaries)
1107 
1108  for(Int_t jTracks = 0; jTracks<fStack->GetEntriesFast(); jTracks++) {
1109  AliAODMCParticle *MCpart = (AliAODMCParticle*)fStack->At(jTracks);
1110  if (!MCpart) {
1111  printf("ERROR: Could not receive MC track %d\n", jTracks);
1112  continue;
1113  }
1114  // kinematic cuts
1115  if ( MCpart->Pt() < .2 || MCpart->Pt() > 20. || TMath::Abs(MCpart->Eta()) > .8 ) continue;
1116  // select charged primaries
1117  if ( MCpart->Charge() == 0. || !MCpart->IsPhysicalPrimary()) continue;
1118 
1119  fPtSpecGen[0][CenBin]->Fill(MCpart->Pt());
1120  }
1121 
1122  // fGenHeader = McEvent->GenEventHeader();
1123  // if(fGenHeader) fPythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(fGenHeader);
1124  // printf("#reconstructed : %d (rejected from cuts %d), #MC primaries : %d (rejected from cuts %d) \n",AODPOIs,AODbads,MCPrims,MCSecos);
1125  fFlowEvent->SetReferenceMultiplicity(aod->GetNumberOfTracks());
1126  fFlowEvent->SetCentrality(centr);
1127  // if (McEvent && McEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(McEvent);
1128  fFlowEvent->SetRun(RunNum);
1129  // printf("Run : %d, RefMult : %d, Cent : %f \n",fFlowEvent->GetRun(),fFlowEvent->GetReferenceMultiplicity(),fFlowEvent->GetCentrality());
1130  }
1131 
1132  if(fAnalysisType == "MCESD") {
1133 
1134  fFlowEvent->ClearFast();
1135 
1136  if(!esd) {
1137  AliError("ERROR: Could not retrieve ESDEvent");
1138  return;
1139  }
1140  if(!McEvent) {
1141  AliError("ERROR: Could not retrieve MCEvent");
1142  return;
1143  }
1144  AliStack* fStack = fMCEvent->Stack();
1145  if(!fStack) {
1146  AliError("ERROR: Could not retrieve MCStack");
1147  return;
1148  }
1149 
1150  AliESDVertex *vertex = (AliESDVertex*) esd->GetPrimaryVertex();
1151  if (!vertex) return;
1152  if (TMath::Abs(vertex->GetZ()) > 10. ) return;
1153  if (vertex->GetNContributors() < 1 ) return;
1154  AliCentrality *centrality = esd->GetCentrality();
1155  if (!centrality) return;
1156  Double_t centr = centrality->GetCentralityPercentile("V0M");
1157  if (centr<fCentrLowLim || centr>=fCentrUpLim ) return;
1158  Int_t CenBin = -1;
1159  if (centr>0. && centr<5.) CenBin=0;
1160  if (centr>5. && centr<10.) CenBin=1;
1161  if (centr>10. && centr<20.) CenBin=2;
1162  if (centr>20. && centr<30.) CenBin=3;
1163  if (centr>30. && centr<40.) CenBin=4;
1164  if (centr>40. && centr<50.) CenBin=5;
1165  if (centr>50. && centr<60.) CenBin=6;
1166  if (centr>60. && centr<70.) CenBin=7;
1167  if (centr>70. && centr<80.) CenBin=8;
1168  if (centr>80. && centr<90.) CenBin=9;
1169  if(CenBin==-1) return;
1170 
1171  //Generated
1172  Int_t MCPrims = 0;
1173  for ( Int_t i=0 ; i<fStack->GetNtrack() ; i++ ) {
1174 
1175  //Primaries Selection
1176  TParticle *particle = (TParticle*)fStack->Particle(i);
1177  if (!particle) continue;
1178  if (!fStack->IsPhysicalPrimary(i)) continue;
1179  if ( particle->GetPDG()->Charge() == 0.) continue;
1180 
1181  //Kinematic Cuts
1182  if ( particle->Pt()<0.2 || particle->Pt()>10. ) continue;
1183  if ( TMath::Abs(particle->Eta())>0.8 ) continue;
1184 
1185  fFlowTrack->SetPhi(particle->Phi());
1186  fFlowTrack->SetEta(particle->Eta());
1187  fFlowTrack->SetPt(particle->Pt());
1189  fFlowTrack->SetForRPSelection(kTRUE);
1191  fFlowTrack->SetForPOISelection(kFALSE);
1193  MCPrims++;
1194 
1195  fPtSpecGen[0][CenBin]->Fill(particle->Pt());
1196 
1197  }
1198 
1199  //Reconstructed
1200  Int_t ESDPrims = 0;
1201  for (Int_t i=0 ; i<esd->GetNumberOfTracks() ; i++) {
1202 
1203  //Get reconstructed track
1204  AliVTrack *vtrack = static_cast<AliVTrack*>(esd->GetTrack(i));
1205  AliESDtrack *track = dynamic_cast<AliESDtrack*>(vtrack);
1206  if (!track) continue;
1207 
1208  //Primaries selection
1209  Int_t lp = TMath::Abs(track->GetLabel());
1210  if (!fStack->IsPhysicalPrimary(lp)) continue;
1211  TParticle *particle = (TParticle*)fStack->Particle(lp);
1212  if (!particle) continue;
1213  if (particle->GetPDG()->Charge() == 0.) continue;
1214 
1215  // if(!fCutsPOI->PassesESDcuts(track)) continue;
1216 
1217  Bool_t pass = kTRUE;
1218 
1219  if(fCutTPC) {
1220  // printf("******* cutting TPC ******** \n");
1221  UShort_t ntpccls = track->GetTPCNcls();
1222  Double_t tpcchi2 = track->GetTPCchi2();
1223  if (tpcchi2<0.2 || tpcchi2 >=4.) {
1224  // printf("TPCchi2 : %e %e ",tpcchi2,track->GetTPCchi2Iter1());
1225  pass=kFALSE;
1226  }
1227  if (ntpccls < 70) {
1228  // printf("#TPCcluster : %u %u %u %u ",ntpccls,track->GetTPCNclsF(),track->GetTPCNclsFIter1(),track->GetTPCNclsIter1());
1229  pass=kFALSE;
1230  }
1231  }
1232 
1233  Float_t dcaxy=0.0;
1234  Float_t dcaz=0.0;
1235  track->GetImpactParameters(dcaxy,dcaz);
1236  if (dcaxy > 0.3 || dcaz > 0.3) {
1237  // printf("DCA : %e %e ",dcaxy,dcaz);
1238  pass=kFALSE;
1239  }
1240  if(!pass) continue;
1241 
1242  //Kinematic Cuts
1243  if ( track->Pt()<0.2 || track->Pt()>10. ) continue;
1244  if ( TMath::Abs(track->Eta())>0.8 ) continue;
1245 
1246  fFlowTrack->SetPhi(track->Phi());
1247  fFlowTrack->SetEta(track->Eta());
1248  fFlowTrack->SetPt(track->Pt());
1250  fFlowTrack->SetForRPSelection(kFALSE);
1254  ESDPrims++;
1255 
1256  }
1257 
1258  // printf("#reconstructed : %d , #MC primaries : %d \n",ESDPrims,MCPrims);
1259  fFlowEvent->SetReferenceMultiplicity(esd->GetNumberOfTracks());
1260  fFlowEvent->SetCentrality(centr);
1261  if (McEvent && McEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(McEvent);
1262  fFlowEvent->SetRun(esd->GetRunNumber());
1263  // printf("Run : %d, RefMult : %d, Cent : %f \n",fFlowEvent->GetRun(),fFlowEvent->GetReferenceMultiplicity(),fFlowEvent->GetCentrality());
1264 
1265  } // end of if(fAnalysisType == "MCESD")
1266 
1267  if(fAnalysisType == "MCkine") {
1268 
1269  fFlowEvent->ClearFast();
1270 
1271  AliInputEventHandler* McHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
1272  if(!McHandler) {
1273  AliError("ERROR: Could not retrieve MCtruthEventHandler");
1274  return;
1275  }
1276  McEvent = McHandler->MCEvent();
1277  if(!McEvent) {
1278  AliError("ERROR: Could not retrieve MC event");
1279  return;
1280  }
1281 
1282  Int_t nTracks = McEvent->GetNumberOfTracks();
1283  // Int_t nPrimTr = McEvent->GetNumberOfPrimaries();
1284 
1285  //loop over tracks
1286  for (Int_t itrkN=0; itrkN<nTracks; itrkN++) {
1287  //get input particle
1288  AliMCParticle* pParticle = dynamic_cast<AliMCParticle*>(McEvent->GetTrack(itrkN));
1289  if (!pParticle) continue;
1290 
1291  //check if track passes the cuts
1292  if (McEvent->IsPhysicalPrimary(itrkN) && pParticle->Charge()!=0) {
1293  fFlowTrack->Set(pParticle);
1295  fFlowTrack->SetForRPSelection(kTRUE);
1300  }
1301  }// for all tracks
1302 
1303  // if monte carlo event get reaction plane from monte carlo (depends on generator)
1304  if (McEvent && McEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(McEvent);
1305  // set reference multiplicity
1306  fFlowEvent->SetReferenceMultiplicity(McEvent->GetNumberOfTracks());
1307  // tag subevents
1309  // set centrality from impact parameter
1310  Double_t ImpPar=0., CenPer=0.;
1311  fGenHeader = McEvent->GenEventHeader();
1312  if(fGenHeader){
1313  fPythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(fGenHeader);
1314  if(fPythiaGenHeader) ImpPar = fPythiaGenHeader->GetImpactParameter();
1315  fHijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(fGenHeader);
1316  if(fHijingGenHeader) ImpPar = fHijingGenHeader->ImpactParameter();
1317  if(ImpPar) CenPer = 0.3859796743103508*pow(ImpPar,2.);
1318  if(CenPer>0. && CenPer<100.) fFlowEvent->SetCentrality(CenPer);
1319  else return;
1320  fFlowEvent->SetRun(1);
1321  }
1322 
1323  } // end of if(fAnalysisType == "MCkine")
1324 
1325  if (!fFlowEvent) return; //shuts up coverity
1326 
1327  //check final event cuts
1328  Int_t mult = fFlowEvent->NumberOfTracks();
1329  // AliInfo(Form("FlowEvent has %i tracks",mult));
1330  if (mult<fMinMult || mult>fMaxMult) {
1331  AliWarning("FlowEvent cut on multiplicity"); return;
1332  }
1333 
1334  //define dead zone
1336 
1339  if (fAfterburnerOn)
1340  {
1341  //if reaction plane not set from elsewhere randomize it before adding flow
1343  fFlowEvent->SetMCReactionPlaneAngle(gRandom->Uniform(0.0,TMath::TwoPi()));
1344 
1345  if(fDifferentialV2)
1347  else
1348  fFlowEvent->AddFlow(fV1,fV2,fV3,fV4,fV5); //add flow
1349  fFlowEvent->CloneTracks(fNonFlowNumberOfTrackClones); //add nonflow by cloning tracks
1350  }
1352 
1353  //tag subEvents
1355 
1356  //do we want to serve shullfed tracks to everybody?
1358 
1359  // associate the mother particles to their daughters in the flow event (if any)
1361 
1362  //fListHistos->Print();
1363  //fOutputFile->WriteObject(fFlowEvent,"myFlowEventSimple");
1364 
1365  //********************************************************************************************************************************
1366 
1367  if(fAnalysisType == "AOD" || fAnalysisType == "AUTOMATIC") {
1368 
1369  // PHYSICS SELECTION
1370  AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1371  AliInputEventHandler *hdr = (AliInputEventHandler*)am->GetInputEventHandler();
1372 
1373  if(hdr->IsEventSelected() && AliVEvent::kAny) {
1374 
1375  Double_t centrperc = fFlowEvent->GetCentrality();
1376  Int_t cenb = (Int_t)centrperc;
1377 
1378  AliAODTracklets *trackl = aod->GetTracklets();
1379  Int_t nTracklets = trackl->GetNumberOfTracklets();
1380 
1381  AliAODVZERO *vzeroAOD = aod->GetVZEROData();
1382  Double_t multV0A = vzeroAOD->GetMTotV0A();
1383  Double_t multV0C = vzeroAOD->GetMTotV0C();
1384 
1385  AliAODZDC *aodZDC = aod->GetZDCData();
1386 
1387  Double_t energyZNC = (Double_t) (aodZDC->GetZNCEnergy());
1388  Double_t energyZPC = (Double_t) (aodZDC->GetZPCEnergy());
1389  Double_t energyZNA = (Double_t) (aodZDC->GetZNAEnergy());
1390  Double_t energyZPA = (Double_t) (aodZDC->GetZPAEnergy());
1391  Double_t energyZEM1 = (Double_t) (aodZDC->GetZEM1Energy());
1392  Double_t energyZEM2 = (Double_t) (aodZDC->GetZEM2Energy());
1393 
1394  const Double_t * towZNCraw = aodZDC->GetZNCTowerEnergy();
1395  const Double_t * towZNAraw = aodZDC->GetZNATowerEnergy();
1396 
1397  // Get centroid from ZDCs *******************************************************
1398 
1399  Double_t Enucl = (RunNum < 209122 ? 1380. : 2511.);
1400  Double_t xyZNC[2]={999.,999.}, xyZNA[2]={999.,999.};
1401  Double_t towZNC[5]={0.}, towZNA[5]={0.};
1402 
1403  Double_t ZNCcalib=1., ZNAcalib=1.;
1404  if(fUseTowerEq) {
1405  if(RunNum!=fCachedRunNum) {
1406  for(Int_t i=0; i<5; i++) {
1407  fTowerGainEq[0][i] = (TH1D*)(fTowerEqList->FindObject(Form("fZNCTower[%d][%d]",RunNum,i)));
1408  fTowerGainEq[1][i] = (TH1D*)(fTowerEqList->FindObject(Form("fZNATower[%d][%d]",RunNum,i)));
1409  }
1410  }
1411  for(Int_t i=0; i<5; i++) {
1412  if(fTowerGainEq[0][i]) towZNC[i] = towZNCraw[i]*fTowerGainEq[0][i]->GetBinContent(fTowerGainEq[0][i]->FindBin(centrperc));
1413  if(fTowerGainEq[1][i]) towZNA[i] = towZNAraw[i]*fTowerGainEq[1][i]->GetBinContent(fTowerGainEq[1][i]->FindBin(centrperc));
1414  }
1415  } else {
1416  for(Int_t i=0; i<5; i++) {
1417  towZNC[i] = towZNCraw[i];
1418  towZNA[i] = towZNAraw[i];
1419  fZNCTower[RunBin][i]->Fill(centrperc,towZNC[i]);
1420  fZNATower[RunBin][i]->Fill(centrperc,towZNA[i]);
1421  }
1422  }
1423 
1424  if(RunNum>=245829) towZNA[2] = 0.;
1425  Double_t zncEnergy=0., znaEnergy=0.;
1426  for(Int_t i=0; i<5; i++){
1427  zncEnergy += towZNC[i];
1428  znaEnergy += towZNA[i];
1429  }
1430  if(RunNum>=245829) znaEnergy *= 8./7.;
1431  fFlowEvent->SetZNCEnergy(towZNC[0]);
1432  fFlowEvent->SetZNAEnergy(towZNA[0]);
1433 
1434  const Double_t x[4] = {-1.75, 1.75, -1.75, 1.75};
1435  const Double_t y[4] = {-1.75, -1.75, 1.75, 1.75};
1436  Double_t numXZNC=0., numYZNC=0., denZNC=0., cZNC, wZNC, EZNC;
1437  Double_t numXZNA=0., numYZNA=0., denZNA=0., cZNA, wZNA, EZNA, BadChOr;
1438 
1439  if (fUseMCCen) {
1440  for(Int_t i=0; i<4; i++){
1441 
1442  // get energy
1443  EZNC = towZNC[i+1];
1444  fhZNSpectra->Fill(centrperc,i+0.5,EZNC);
1445  if(fUseZDCSpectraCorr) {
1446  Double_t mu1 = SpecCorMu1[i]->Interpolate(centrperc);
1447  Double_t mu2 = SpecCorMu2[i]->Interpolate(centrperc);
1448  Double_t av = SpecCorAv[i]->Interpolate(centrperc);
1449  Double_t cor1 = FitSpecCorSi[i]->Eval(centrperc);
1450  EZNC = exp( (log(EZNC) - mu1 + mu2*cor1)/cor1 ) + av;
1451  fhZNSpectraCor->Fill(centrperc,i+0.5,EZNC);
1452  }
1453 
1454  // build centroid
1455  wZNC = TMath::Power(EZNC, fZDCGainAlpha);
1456  numXZNC += x[i]*wZNC;
1457  numYZNC += y[i]*wZNC;
1458  denZNC += wZNC;
1459  fhZNSpectraPow->Fill(centrperc,i+0.5,wZNC);
1460 
1461  // get energy
1462  if(i==1) {
1463  EZNA = towZNA[0]-towZNA[1]-towZNA[3]-towZNA[4];
1464  if(fUseBadTowerCalib && fBadTowerCalibHist[cenb]) {
1465  EZNA = GetBadTowerResp(EZNA, fBadTowerCalibHist[cenb]);
1466  }
1467  } else {
1468  EZNA = towZNA[i+1];
1469  }
1470  fhZNSpectra->Fill(centrperc,i+4.5,EZNA);
1471  if(fUseZDCSpectraCorr) {
1472  Double_t mu1 = SpecCorMu1[i+4]->Interpolate(centrperc);
1473  Double_t mu2 = SpecCorMu2[i+4]->Interpolate(centrperc);
1474  Double_t av = SpecCorAv[i+4]->Interpolate(centrperc);
1475  Double_t cor1 = FitSpecCorSi[i+4]->Eval(centrperc);
1476  EZNA = exp( (log(EZNA) - mu1 + mu2*cor1)/cor1 ) + av;
1477  fhZNSpectraCor->Fill(centrperc,i+4.5,EZNA);
1478  }
1479 
1480  // build centroid
1481  wZNA = TMath::Power(EZNA, fZDCGainAlpha);
1482  numXZNA += x[i]*wZNA;
1483  numYZNA += y[i]*wZNA;
1484  denZNA += wZNA;
1485  fhZNSpectraPow->Fill(centrperc,i+4.5,wZNA);
1486  }
1487  // store distribution for unfolding
1488  if(RunNum<245829) {
1489  Double_t recoE = towZNA[0]-towZNA[1]-towZNA[3]-towZNA[4];
1490  Double_t trueE = towZNA[2];
1491  fhZNBCCorr->Fill(centrperc,trueE,recoE);
1492  }
1493  if(denZNC!=0){
1494  Double_t nSpecnC = zncEnergy/Enucl;
1495  cZNC = 1.89358-0.71262/(nSpecnC+0.71789);
1496  xyZNC[0] = cZNC*numXZNC/denZNC;
1497  xyZNC[1] = cZNC*numYZNC/denZNC;
1498  }
1499  else{
1500  xyZNC[0] = xyZNC[1] = 999.;
1501  }
1502  if(denZNA!=0){
1503  Double_t nSpecnA = znaEnergy/Enucl;
1504  cZNA = 1.89358-0.71262/(nSpecnA+0.71789);
1505  xyZNA[0] = cZNA*numXZNA/denZNA;
1506  xyZNA[1] = cZNA*numYZNA/denZNA;
1507  }
1508  else{
1509  xyZNA[0] = xyZNA[1] = 999.;
1510  }
1511  } else {
1512  for(Int_t i=0; i<4; i++) {
1513  if(towZNC[i+1]>0.) {
1514  wZNC = TMath::Power(towZNC[i+1], fZDCGainAlpha);
1515  numXZNC += x[i]*wZNC;
1516  numYZNC += y[i]*wZNC;
1517  denZNC += wZNC;
1518  }
1519  if(towZNA[i+1]>0.) {
1520  wZNA = TMath::Power(towZNA[i+1], fZDCGainAlpha);
1521  numXZNA += x[i]*wZNA;
1522  numYZNA += y[i]*wZNA;
1523  denZNA += wZNA;
1524  }
1525  }
1526  if(denZNC!=0) {
1527  xyZNC[0] = numXZNC/denZNC;
1528  xyZNC[1] = numYZNC/denZNC;
1529  }
1530  else{
1531  xyZNC[0] = xyZNC[1] = 999.;
1532  zncEnergy = 0.;
1533  }
1534  if(denZNA!=0) {
1535  xyZNA[0] = numXZNA/denZNA;
1536  xyZNA[1] = numYZNA/denZNA;
1537  }
1538  else{
1539  xyZNA[0] = xyZNA[1] = 999.;
1540  znaEnergy = 0.;
1541  }
1542  }
1543 
1544  fhZNCcentroid->Fill(xyZNC[0], xyZNC[1]);
1545  fhZNAcentroid->Fill(xyZNA[0], xyZNA[1]);
1546  fFlowEvent->SetZDC2Qsub(xyZNC,denZNC,xyZNA,denZNA);
1547 
1548  // ******************************************************************************
1549 
1550  Double_t tdcSum = aodZDC->GetZDCTimeSum();
1551  Double_t tdcDiff = aodZDC->GetZDCTimeDiff();
1552  fhDebunch->Fill(tdcDiff, tdcSum);
1553 
1554  for(int i=0; i<5; i++){
1555  fhZNCPM[i]->Fill(towZNC[i]);
1556  if((i<4) && (towZNC[0]>0.)) fhZNCPMQiPMC[i]->Fill(towZNC[i+1]/towZNC[0]);
1557  }
1558  for(int i=0; i<5; i++){
1559  fhZNAPM[i]->Fill(towZNA[i]);
1560  if(((i<4) && towZNA[0]>0.)) fhZNAPMQiPMC[i]->Fill(towZNA[i+1]/towZNA[0]);
1561  }
1562 
1563  fhZNCvsZNA->Fill(energyZNA, energyZNC);
1564  fhZDCCvsZDCCA->Fill(energyZNA+energyZPA, energyZNC+energyZPC);
1565  fhZNCvsZPC->Fill(energyZPC, energyZNC);
1566  fhZNAvsZPA->Fill(energyZPA, energyZNA);
1567  fhZNvsZP->Fill(energyZPA+energyZPC, energyZNA+energyZNC);
1568  fhZNvsVZERO->Fill(multV0A+multV0C, energyZNC+energyZNA);
1569  fhZDCvsVZERO->Fill(multV0A+multV0C, energyZNA+energyZPA+energyZNC+energyZPC);
1570  fhZDCvsTracklets->Fill((Double_t) (nTracklets), energyZNA+energyZPA+energyZNC+energyZPC);
1571 
1572  Double_t asymmetry = -999.;
1573  if((energyZNC+energyZNA)>0.) asymmetry = (energyZNC-energyZNA)/(energyZNC+energyZNA);
1574  fhAsymm->Fill(asymmetry);
1575  fhZNAvsAsymm->Fill(asymmetry, energyZNA/1000.);
1576  fhZNCvsAsymm->Fill(asymmetry, energyZNC/1000.);
1577 
1578  fhZNCvscentrality->Fill(centrperc, energyZNC/1000.);
1579  fhZNAvscentrality->Fill(centrperc, energyZNA/1000.);
1580 
1581  } // PHYSICS SELECTION
1582 
1583  }
1584 
1585  // p) cache run number
1587 
1588  // printf("debug: NoRPs %e, NoPOIs %e, RunNum %d, Cen %e \n",fFlowEvent->GetNumberOfRPs(),fFlowEvent->GetNumberOfPOIs(),fCachedRunNum,fFlowEvent->GetCentrality());
1589 
1590  PostData(1, fFlowEvent);
1591 
1592  PostData(2, fOutput);
1593 }
1594 //________________________________________________________________________
1595 
1597 {
1598  Double_t EtC = BadTowerCalibHist->ProjectionY("",BadTowerCalibHist->GetXaxis()->FindBin(Et),BadTowerCalibHist->GetXaxis()->FindBin(Et))->GetRandom();
1599  return EtC;
1600 }
1601 
1602 //________________________________________________________________________
1603 
1605 {
1606  Int_t CenBin=-1;
1607  if (Centrality>0. && Centrality<5.) CenBin=0;
1608  if (Centrality>5. && Centrality<10.) CenBin=1;
1609  if (Centrality>10. && Centrality<20.) CenBin=2;
1610  if (Centrality>20. && Centrality<30.) CenBin=3;
1611  if (Centrality>30. && Centrality<40.) CenBin=4;
1612  if (Centrality>40. && Centrality<50.) CenBin=5;
1613  if (Centrality>50. && Centrality<60.) CenBin=6;
1614  if (Centrality>60. && Centrality<70.) CenBin=7;
1615  if (Centrality>70. && Centrality<80.) CenBin=8;
1616  if (Centrality>80. && Centrality<90.) CenBin=9;
1617  if (CenBin>=fnCen) CenBin=-1;
1618  if (fnCen==1) CenBin=0;
1619  return CenBin;
1620 } // end of AliFlowAnalysisCRC::GetCRCCenBin(Double_t Centrality)
1621 //_____________________________________________________________________________
1622 
1623 Double_t AliAnalysisTaskCRCZDC::GetWDist(const AliVVertex* v0, const AliVVertex* v1)
1624 {
1625  // calculate sqrt of weighted distance to other vertex
1626  if (!v0 || !v1) {
1627  printf("One of vertices is not valid\n");
1628  return 0;
1629  }
1630  static TMatrixDSym vVb(3);
1631  double dist = -1;
1632  double dx = v0->GetX()-v1->GetX();
1633  double dy = v0->GetY()-v1->GetY();
1634  double dz = v0->GetZ()-v1->GetZ();
1635  double cov0[6],cov1[6];
1636  v0->GetCovarianceMatrix(cov0);
1637  v1->GetCovarianceMatrix(cov1);
1638  vVb(0,0) = cov0[0]+cov1[0];
1639  vVb(1,1) = cov0[2]+cov1[2];
1640  vVb(2,2) = cov0[5]+cov1[5];
1641  vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
1642  vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
1643  vVb.InvertFast();
1644  if (!vVb.IsValid()) {printf("Singular Matrix\n"); return dist;}
1645  dist = vVb(0,0)*dx*dx + vVb(1,1)*dy*dy + vVb(2,2)*dz*dz
1646  + 2*vVb(0,1)*dx*dy + 2*vVb(0,2)*dx*dz + 2*vVb(1,2)*dy*dz;
1647  return dist>0 ? TMath::Sqrt(dist) : -1;
1648 }
1649 //________________________________________________________________________
1650 
1652 {
1653  // check for multi-vertexer pile-up
1654 
1655  const int kMinPlpContrib = 5;
1656  const double kMaxPlpChi2 = 5.0;
1657  const double kMinWDist = 15;
1658 
1659  const AliVVertex* vtPrm = 0;
1660  const AliVVertex* vtPlp = 0;
1661  int nPlp = 0;
1662 
1663  if ( !(nPlp=aod->GetNumberOfPileupVerticesTracks()) ) return kFALSE;
1664  vtPrm = aod->GetPrimaryVertex();
1665  if (vtPrm == aod->GetPrimaryVertexSPD()) return kTRUE; // there are pile-up vertices but no primary
1666 
1667  //int bcPrim = vtPrm->GetBC();
1668 
1669  for (int ipl=0;ipl<nPlp;ipl++) {
1670  vtPlp = (const AliVVertex*)aod->GetPileupVertexTracks(ipl);
1671  //
1672  if (vtPlp->GetNContributors() < kMinPlpContrib) continue;
1673  if (vtPlp->GetChi2perNDF() > kMaxPlpChi2) continue;
1674  // int bcPlp = vtPlp->GetBC();
1675  // if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2) return kTRUE; // pile-up from other BC
1676  //
1677  double wDst = GetWDist(vtPrm,vtPlp);
1678  if (wDst<kMinWDist) continue;
1679  //
1680  return kTRUE; // pile-up: well separated vertices
1681  }
1682 
1683  return kFALSE;
1684 }
1685 
1686 //________________________________________________________________________
1688  fCutContainer->AddAt(cutsRP,0); fCutsRP=cutsRP; cutsRP->SetPOItype(0);
1689 }
1690 
1691 //________________________________________________________________________
1693  fCutContainer->AddAt(cutsPOI,1); fCutsPOI=cutsPOI; cutsPOI->SetPOItype(1);
1694 }
1695 
1696 //________________________________________________________________________
1698 {
1699  // Terminate analysis
1700  //
1701  /* if(fDebug > 1) printf(" **** AliAnalysisTaskCRCZDC::Terminate() \n");
1702 
1703  //fOutput = dynamic_cast<TList*> (GetOutputData(1));
1704  //if(!fOutput) printf("ERROR: fOutput not available\n");
1705  */
1706 }
TH2F * fhZDCCvsZDCCA
ZNC vs ZNA;.
TH2F * fhZNvsVZERO
ZNC+ZNA vs ZPC+ZPA;.
TH1F * fPtSpecFB128[2][10]
PtSpecRec FB96.
void FindDaughters(Bool_t keepDaughtersInRPselection=kFALSE)
TH3D * fhZNSpectraCor
ZNA vs. centrality.
const Color_t cc[]
Definition: DrawKs.C:1
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.
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)
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 SetReferenceMultiplicity(Int_t m)
virtual Int_t GetCenBin(Double_t Centrality)
void SetCutsRP(AliFlowTrackCuts *cutsRP)
TRandom * gRandom
TH2F * fhZNCvsZNA
PMQi/PMC for ZNA.
TH2F * fhDebunch
ZDC vs N_cluster layer 1;.
TH1F * fhAsymm
ZNA centroid.
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.
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
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)
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)
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 * fhZNAcentroid
ZNC centroid.
Bool_t IsSetMCReactionPlaneAngle() const
TClonesArray * fStack
ZNA tower spectra.
TList * GetQA() const
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)
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 fCutTPC
PtSpecRec FB768.
AliFlowTrackCuts * fCutsRP
void InsertTrack(AliFlowTrack *)
virtual void ClearFast()
Bool_t fUseBadTowerCalib
list for storing calib files
AliAnalysisUtils * fAnalysisUtil
Event selection.
TH2F * fhZNCcentroid
Debunch;.
TH2F * fhZNvsZP
ZNA vs ZPA;.
AliGenHijingEventHeader * fHijingGenHeader
const Double_t phimin
AliGenPythiaEventHeader * fPythiaGenHeader
Int_t NumberOfTracks() const