AliPhysics  6cf2591 (6cf2591)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskVnZDC.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 the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
17 // AliAnalysisTaskVnZDC:
18 // Author: Rihan Haque (mhaque@cern.ch)
20 
21 
22 #include "Riostream.h" //needed as include
23 
24 class TFile;
25 class TList;
26 class AliAnalysisTaskSE;
27 
28 #include "TProfile.h" //needed as include
29 #include "TNtuple.h"
30 #include "TProfile2D.h"
31 #include "TApplication.h"
32 #include "TList.h"
33 #include "TH1.h"
34 #include "TH2.h"
35 #include "AliMultSelection.h"
36 #include "AliVVertex.h"
37 
38 #include "AliAnalysisManager.h"
39 #include "AliFlowEventSimple.h"
40 #include "AliAnalysisTaskVnZDC.h"
41 #include "AliFlowCommonHist.h"
43 
44 #include "AliAnalysisManager.h"
45 #include "AliInputEventHandler.h"
46 #include "AliVEvent.h"
47 #include "AliESD.h"
48 #include "AliESDEvent.h"
49 #include "AliESDHeader.h"
50 #include "AliESDInputHandler.h"
51 #include "AliESDZDC.h"
52 #include "AliMultiplicity.h"
53 #include "AliAnalysisUtils.h"
54 #include "AliAODHandler.h"
55 #include "AliAODTrack.h"
56 #include "AliAODEvent.h"
57 #include "AliAODHeader.h"
58 #include "AliAODVertex.h"
59 #include "AliAODVZERO.h"
60 #include "AliAODZDC.h"
61 #include "AliAODMCHeader.h"
62 #include "AliMCEventHandler.h"
63 #include "AliMCEvent.h"
64 #include "AliHeader.h"
65 #include "AliVParticle.h"
66 #include "AliStack.h"
67 #include "AliAODMCParticle.h"
68 #include "AliAnalysisTaskSE.h"
69 #include "AliGenEventHeader.h"
70 #include "AliPhysicsSelectionTask.h"
71 #include "AliPhysicsSelection.h"
72 #include "AliBackgroundSelection.h"
73 #include "AliTriggerAnalysis.h"
74 #include "AliCentrality.h"
75 #include "AliLog.h"
76 
77 using std::endl;
78 using std::cout;
79 
80 
82 
83 //________________________________________________________________________
84 AliAnalysisTaskVnZDC::AliAnalysisTaskVnZDC(const char *name, Bool_t usePhiWeights) :
85 AliAnalysisTaskSE(name),
86 fEvent(NULL),
87 fListHistos(NULL),
88 fMinimalBook(kFALSE),
89 fUsePhiWeights(usePhiWeights),
90 fListWeights(NULL),
91 fRelDiffMsub(1.0),
92 fApplyCorrectionForNUA(kFALSE),
93 fHarmonic(2),
94 fNormalizationType(1),
95 fNCentBins(9),
96 fTotalQvector(NULL),
97 fievent(0),
98 frunflag(0),
99 fHist_Vertex_XY(NULL),
100 fHist_Vx_vs_runnum(NULL),
101 fHist_Vy_vs_runnum(NULL),
102 fHist_Vz_vs_runnum(NULL),
103 fHist_tracks_vs_runnum(NULL),
104 fHist_Vx_ArrayFinder(NULL),
105 fHist_Vy_ArrayFinder(NULL),
106 fHist_Vz_ArrayFinder(NULL),
107 fHist_Vertex_Vz(NULL),
108 fHist_Event_count(NULL),
109 fHist_Cent_count1(NULL),
110 fHist_Cent_count2(NULL),
111 fHist_EventperRun(NULL),
112 fHist_Psi1_zdnA(NULL),
113 fHist_Psi1_zdnC(NULL),
114 fHist_vBincount(NULL),
115 fZDCESEList(NULL),
116 fPileUpMultSelCount(NULL),
117 fPileUpCount(NULL),
118 fMultSelection(NULL),
119 fAnalysisUtil(NULL),
120 fRejectPileUpTight(kFALSE),
121 fRejectPileUp(kFALSE),
122 checkOnce(0),
123 vxBin(0),
124 vyBin(0),
125 vzBin(0)
126 {
127  for(int i=0;i<89;i++)
128  runNums[i] = 0;
129 
130  for(int i=0;i<89;i++){
131  fHist_znCx_V0_Cent[i] = NULL;
132  fHist_znCy_V0_Cent[i] = NULL;
133  fHist_znAx_V0_Cent[i] = NULL;
134  fHist_znAy_V0_Cent[i] = NULL;
135 
136  fHist_znCx_V0_VxVy[i] = NULL;
137  fHist_znCy_V0_VxVy[i] = NULL;
138  fHist_znAx_V0_VxVy[i] = NULL;
139  fHist_znAy_V0_VxVy[i] = NULL;
140 
141  fHist_ZDCn_A_XY[i] = NULL;
142  fHist_ZDCn_C_XY[i] = NULL;
143  fHist_ZDCA_En_Run[i] = NULL;
144  fHist_ZDCC_En_Run[i] = NULL;
145  }
146 
147  for(int i=0;i<2;i++)
148  {
149  VxCut[i] = 0;
150  VyCut[i] = 0;
151  VzCut[i] = 0;
152  }
153 
154  DefineInput(1, AliFlowEventSimple::Class()); // Input slot #0 works with an AliFlowEventSimple
155 
156  DefineOutput(1,TList::Class());
157 
158  if(usePhiWeights) {
159  DefineInput(2, TList::Class()); // Input slot #2 is needed for the weights input file
160  }
161 
162  fDataSet="2010";
163  fAnalysisSet="recenter1";
164 
165  fTotalQvector = new TString("QaQb"); // Total Q-vector is: "QaQb" (means Qa+Qb), "Qa" or "Qb"
166  AliDebug(2,"AliAnalysisTaskVnZDC::Constructor is called..!!)");
167 
168 }//-------------constructor-----------------
169 
170 //________________________________________________
173 fEvent(NULL),
174 fListHistos(NULL),
175 fMinimalBook(kFALSE),
176 fUsePhiWeights(kFALSE),
177 fListWeights(NULL),
178 fRelDiffMsub(1.0),
179 fApplyCorrectionForNUA(kFALSE),
180 fHarmonic(2),
181 fNormalizationType(1),
182 fNCentBins(9),
183 fTotalQvector(NULL),
184 fievent(0),
185 frunflag(0),
186 fHist_Vertex_XY(NULL),
187 fHist_Vx_vs_runnum(NULL),
188 fHist_Vy_vs_runnum(NULL),
189 fHist_Vz_vs_runnum(NULL),
190 fHist_tracks_vs_runnum(NULL),
191 fHist_Vx_ArrayFinder(NULL),
192 fHist_Vy_ArrayFinder(NULL),
193 fHist_Vz_ArrayFinder(NULL),
194 fHist_Vertex_Vz(NULL),
195 fHist_Event_count(NULL),
196 fHist_Cent_count1(NULL),
197 fHist_Cent_count2(NULL),
198 fHist_EventperRun(NULL),
199 fHist_Psi1_zdnA(NULL),
200 fHist_Psi1_zdnC(NULL),
201 fHist_vBincount(NULL),
202 fZDCESEList(NULL),
203 fPileUpMultSelCount(NULL),
204 fPileUpCount(NULL),
205 fMultSelection(NULL),
206 fAnalysisUtil(NULL),
207 fRejectPileUpTight(kFALSE),
208 fRejectPileUp(kFALSE),
209 checkOnce(0),
210 vxBin(0),
211 vyBin(0),
212 vzBin(0)
213 {
214  for(int i=0;i<89;i++)
215  runNums[i] = 0;
216 
217  for(int i=0;i<89;i++){
218  fHist_znCx_V0_Cent[i] = NULL;
219  fHist_znCy_V0_Cent[i] = NULL;
220  fHist_znAx_V0_Cent[i] = NULL;
221  fHist_znAy_V0_Cent[i] = NULL;
222 
223  fHist_znCx_V0_VxVy[i] = NULL;
224  fHist_znCy_V0_VxVy[i] = NULL;
225  fHist_znAx_V0_VxVy[i] = NULL;
226  fHist_znAy_V0_VxVy[i] = NULL;
227 
228  fHist_ZDCA_En_Run[i] = NULL;
229  fHist_ZDCC_En_Run[i] = NULL;
230  fHist_ZDCn_A_XY[i] = NULL;
231  fHist_ZDCn_C_XY[i] = NULL;
232  }
233 
234  for(int i=0;i<2;i++)
235  {
236  VxCut[i] = 0;
237  VyCut[i] = 0;
238  VzCut[i] = 0;
239  }
240 
241  // AliDebug(2,"AliAnalysisTaskVnZDC::AliAnalysisTaskVnZDC()");
242 }
243 
244 
245 
246 
247 
248 //________________________________________________________________________
250 {
251  //delete [] fSP; //is this not deleting all fSP[i] ?
252  //delete fListHistos;
253  delete fHist_Vx_ArrayFinder;
254  delete fHist_Vy_ArrayFinder;
255  delete fHist_Vz_ArrayFinder;
256  delete fListHistos;
257  delete fZDCESEList;
258  delete fMultSelection;
259  delete fAnalysisUtil;
260  delete fEvent;
261 
262 
263  if(fAnalysisSet=="recenter2"){
264  for(int j=0;j<89;j++){
265  delete fHist_znCx_V0_VxVy[j];
266  delete fHist_znCy_V0_VxVy[j];
267  delete fHist_znAx_V0_VxVy[j];
268  delete fHist_znAy_V0_VxVy[j];
269  }
270  }
271  if(fAnalysisSet=="recenter1"){
272  for(int j=0;j<89;j++){
273  delete fHist_znCx_V0_Cent[j];
274  delete fHist_znCy_V0_Cent[j];
275  delete fHist_znAx_V0_Cent[j];
276  delete fHist_znAy_V0_Cent[j];
277  }
278  }
279 
280  AliDebug(2,"AliAnalysisTaskVnZDC::Constructor is called..!!)");
281 
282 }
283 
284 //________________________________________________________________________
286 {
287 
288  AliDebug(2,"AliAnalysisTaskVnZDC::CreateOutputObjects() is called....");
289 
290 
291  Int_t runArray_2010[89] = {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, 137243, 137236, 137235, 137232, 137231, 137162, 137161};
292 
293  Int_t runArray_2011[68] = {167915, 168115, 168460, 169035, 169238, 169859, 170228, 167920, 168310, 168464, 169091, 169411, 169923, 170230, 167985, 168311, 168467, 169094, 169415, 170027, 170268, 167987, 168322, 168511, 169138, 169417, 170081, 170269, 167988, 168325, 168512, 169144, 169835, 170155, 170270, 168069, 168341, 168514, 169145, 169837, 170159, 170306, 168076, 168342, 168777, 169148, 169838, 170163, 170308, 168105, 168361, 168826, 169156, 169846, 170193, 170309, 168107, 168362, 168988, 169160, 169855, 170203, 168108, 168458, 168992, 169167, 169858, 170204};
294 
295  Int_t runArray_2015[77] = {246994, 246991, 246989, 246984, 246982, 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, 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, 245692, 245683};
296 
297 
298  if(fDataSet=="2010"){
299  frunflag = 89;
300  for(int i=0;i<frunflag;i++)
301  runNums[i] = runArray_2010[i];
302  }
303  if(fDataSet=="2011"){
304  frunflag = 10; //<--------- 2011 runnumbers have to be checked and put...
305  for(int i=0;i<frunflag;i++)
306  runNums[i] = runArray_2011[i];
307  }
308  if(fDataSet=="2015"){
309  frunflag = 77;
310  for(int i=0;i<frunflag;i++)
311  runNums[i] = runArray_2015[i];
312  }
313 
314 
315  Int_t totalQ = 0;
316  if(fTotalQvector->Contains("Qa")) totalQ += 1;
317  if(fTotalQvector->Contains("Qb")) totalQ += 2;
318 
319  fListHistos = new TList();
320  fListHistos->SetOwner(kTRUE);
321 
322  fHist_Event_count = new TH1F("fHist_Event_count"," ",20,0,20);
324  fHist_vBincount = new TH1F("fHist_vBincount"," ",500,0,500);
326  fHist_EventperRun = new TH1F("fHist_EventperRun","", frunflag,0,frunflag);
328  fHist_Psi1_zdnA = new TH1F("fHist_Psi1_zdnA","", 200, 0,2.*TMath::Pi());
330  fHist_Psi1_zdnC = new TH1F("fHist_Psi1_zdnC","", 200, 0,2.*TMath::Pi());
332 
333  fHist_Vertex_Vz = new TH1F("fHist_Vertex_Vz_dist","Vertex Z Distribution",200, -15, 15);
335  fHist_Vertex_XY = new TH2F("fHist_Vertex_XY_dist","Vertex XY Distributn.",100,-0.2,0.2,100,0.0,0.4);
337 
338  // Vx,Vy,Vz vs runnumber:
339  fHist_Vx_vs_runnum = new TProfile("fHist_Vx_vs_runnum","<Vx>_vs_runnum",frunflag,0,frunflag);
341  fHist_Vy_vs_runnum = new TProfile("fHist_Vy_vs_runnum","<Vy>_vs_runnum",frunflag,0,frunflag);
343  fHist_Vz_vs_runnum = new TProfile("fHist_Vz_vs_runnum","<Vz>_vs_runnum",frunflag,0,frunflag);
345  fHist_tracks_vs_runnum = new TProfile("fHist_tracks_vs_runnum","<nTracks>_vs_runnum",frunflag,0,frunflag);
347 
348 
349  Double_t centRange[12] = {0,5,10,20,30,40,50,60,70,80,90,100};
350  fHist_Cent_count1 = new TH1F("fHist_Cent_before_ZDCcut"," ",100,0,100);
352  fHist_Cent_count2 = new TH1F("fHist_Cent_afterr_ZDCcut"," ",100,0,100);
354 
355  fPileUpCount = new TH1F("fPileUpCount", "fPileUpCount", 9, 0., 9.);
356  fPileUpCount->GetXaxis()->SetBinLabel(1,"plpMV");
357  fPileUpCount->GetXaxis()->SetBinLabel(2,"fromSPD");
358  fPileUpCount->GetXaxis()->SetBinLabel(3,"RefMultiplicityComb08");
359  fPileUpCount->GetXaxis()->SetBinLabel(4,"IncompleteDAQ");
360  fPileUpCount->GetXaxis()->SetBinLabel(5,"abs(V0M-CL1)>7.5");
361  fPileUpCount->GetXaxis()->SetBinLabel(6,"missingVtx");
362  fPileUpCount->GetXaxis()->SetBinLabel(7,"inconsistentVtx");
363  fPileUpCount->GetXaxis()->SetBinLabel(8,"multESDTPCDif");
364  fPileUpCount->GetXaxis()->SetBinLabel(9,"extraPileUpMultSel");
366 
367  fPileUpMultSelCount = new TH1F("fPileUpMultSelCount", "fPileUpMultSelCount", 8, 0., 8.);
368  fPileUpMultSelCount->GetXaxis()->SetBinLabel(1,"IsNotPileup");
369  fPileUpMultSelCount->GetXaxis()->SetBinLabel(2,"IsNotPileupMV");
370  fPileUpMultSelCount->GetXaxis()->SetBinLabel(3,"IsNotPileupInMultBins");
371  fPileUpMultSelCount->GetXaxis()->SetBinLabel(4,"InconsistentVertices");
372  fPileUpMultSelCount->GetXaxis()->SetBinLabel(5,"TrackletVsCluster");
373  fPileUpMultSelCount->GetXaxis()->SetBinLabel(6,"AsymmetricInVZERO");
374  fPileUpMultSelCount->GetXaxis()->SetBinLabel(7,"IncompleteDAQ");
375  fPileUpMultSelCount->GetXaxis()->SetBinLabel(8,"GoodVertex2016");
377 
378 
379 
380 
381  //-------- define and add the recenter histograms ----------------
382 if(fDataSet=="2010") {
383  vxBin = 11;
384  vyBin = 14;
385  vzBin = 20;
386 
387  VxCut[0] = -0.035;
388  VxCut[1] = 0.020;
389  VyCut[0] = 0.144;
390  VyCut[1] = 0.214;
391  VzCut[0] = -10.0;
392  VzCut[1] = 10.0;
393  }
394 
395  if(fDataSet=="2015"){
396  vxBin = 13;
397  vyBin = 12;
398  vzBin = 20;
399 
400  VxCut[0] = 0.060;
401  VxCut[1] = 0.086;
402  VyCut[0] = 0.321;
403  VyCut[1] = 0.345;
404  VzCut[0] = -10.0;
405  VzCut[1] = 10.0;
406  }
407  if(fDataSet=="2011"){
408  AliDebug(2,"\n\n !!** WARNING ***!! \n cuts not defined for DATASET: %s\n ...EXIT...\n\n)");
409  exit(0);
410  }
411 
412 
413  const int NbinVt = (vyBin-1)*vxBin + vxBin;
414 
415  fHist_Vx_ArrayFinder = new TH1F("fHist_Vx_ArrayFinder","",vxBin,VxCut[0],VxCut[1]);
416  fHist_Vy_ArrayFinder = new TH1F("fHist_Vy_ArrayFinder","",vyBin,VyCut[0],VyCut[1]);
417  fHist_Vz_ArrayFinder = new TH1F("fHist_Vz_ArrayFinder","",vzBin,VzCut[0],VzCut[1]);
418 
419 
420  //std::cout<<"\n\n Recenter Histograms found..\n dataset: "<<fDataSet.Data()<<std::endl;
421  //fZDCESEList->ls();
422  //std::cout<<"\n\n";
423 
424 
425  if(fAnalysisSet=="recenter2") {
426  if(!fZDCESEList) {
427  AliDebug(2,"\n\n !!** ERROR ***!! \n \n TList fZDCESEList not found !!\n .......EXIT...... \n\n)");
428  exit(0);
429  }
430  if(fZDCESEList) {
431  for(int i=0;i<frunflag;i++) { //**** very carefull with the array index.. you can always misplace and curse youself..***
432  fHist_znCx_V0_VxVy[i] = (TProfile2D *) fZDCESEList->FindObject(Form("fHist_znCx_V0_Run%d_Vz%d",runNums[i],1));
433  fHist_znCy_V0_VxVy[i] = (TProfile2D *) fZDCESEList->FindObject(Form("fHist_znCy_V0_Run%d_Vz%d",runNums[i],1));
434  fHist_znAx_V0_VxVy[i] = (TProfile2D *) fZDCESEList->FindObject(Form("fHist_znAx_V0_Run%d_Vz%d",runNums[i],1));
435  fHist_znAy_V0_VxVy[i] = (TProfile2D *) fZDCESEList->FindObject(Form("fHist_znAy_V0_Run%d_Vz%d",runNums[i],1));
437  AliDebug(2,"\n\n !!** WARNING *** One/more Recenter histograms NOT FOUND\n ......EXIT.....!! \n\n");
438  exit(0);
439  }
440  }
441  }
442  else{ // if recenter file not found then define empty profile histograms so that code doesn't crash
443  AliDebug(2,"\n\n !!** WARNING *** <Qx> = 0, <Qy> = 0, No Recentering....!! \n\n");
444  for(int i=0;i<frunflag;i++){
445  fHist_znCx_V0_VxVy[i] = new TProfile2D(Form("fHist_znCx_V0_Run%d_Vz%d",runNums[i],1),"",NbinVt,0,NbinVt,vzBin,VzCut[0],VzCut[1],"");
446  fHist_znCy_V0_VxVy[i] = new TProfile2D(Form("fHist_znCy_V0_Run%d_Vz%d",runNums[i],1),"",NbinVt,0,NbinVt,vzBin,VzCut[0],VzCut[1],"");
447  fHist_znAx_V0_VxVy[i] = new TProfile2D(Form("fHist_znAx_V0_Run%d_Vz%d",runNums[i],1),"",NbinVt,0,NbinVt,vzBin,VzCut[0],VzCut[1],"");
448  fHist_znAy_V0_VxVy[i] = new TProfile2D(Form("fHist_znAy_V0_Run%d_Vz%d",runNums[i],1),"",NbinVt,0,NbinVt,vzBin,VzCut[0],VzCut[1],"");
449  }
450  }
451 
452  for(int i=0;i<frunflag;i++){
453  //store: X,Y position for recenter:
454  fHist_znCx_V0_Cent[i] = new TProfile(Form("fHist_znCx_V0_Run%d_Cent%d",runNums[i],1),"",200,0,100,"");
455  fHist_znCy_V0_Cent[i] = new TProfile(Form("fHist_znCy_V0_Run%d_Cent%d",runNums[i],1),"",200,0,100,"");
456  fHist_znAx_V0_Cent[i] = new TProfile(Form("fHist_znAx_V0_Run%d_Cent%d",runNums[i],1),"",200,0,100,"");
457  fHist_znAy_V0_Cent[i] = new TProfile(Form("fHist_znAy_V0_Run%d_Cent%d",runNums[i],1),"",200,0,100,"");
458  }
459  for(int i=0;i<frunflag;i++){
464  }
465  }//--------------"recenter2"------------------
466 
467 
468 
469 
470  if(fAnalysisSet=="recenter1") {
471  for(int i=0;i<frunflag;i++){
472  //store: X,Y position run by run:
473  fHist_ZDCn_A_XY[i] = new TH2F(Form("fHist_zdcA_XY_Run%d",runNums[i]),Form("ZDC_A X-Y Run%d",runNums[i]),80,-2,2,80,-2,2);
474  fHist_ZDCn_C_XY[i] = new TH2F(Form("fHist_zdcC_XY_Run%d",runNums[i]),Form("ZDC_C X-Y Run%d",runNums[i]),80,-2,2,80,-2,2);
475  //store: ZDC energy for gain calibration:
476  fHist_ZDCA_En_Run[i] = new TProfile2D(Form("fHist_ZDCA_En_Run%d",runNums[i]),"",100,0,100,5,0,5,"");
477  fHist_ZDCC_En_Run[i] = new TProfile2D(Form("fHist_ZDCC_En_Run%d",runNums[i]),"",100,0,100,5,0,5,"");
478  //store: <X>,<Y> run by run for recenter:
479  fHist_znCx_V0_VxVy[i] = new TProfile2D(Form("fHist_znCx_V0_Run%d_Vz%d",runNums[i],1),"",NbinVt,0,NbinVt,vzBin,VzCut[0],VzCut[1],"");
480  fHist_znCy_V0_VxVy[i] = new TProfile2D(Form("fHist_znCy_V0_Run%d_Vz%d",runNums[i],1),"",NbinVt,0,NbinVt,vzBin,VzCut[0],VzCut[1],"");
481  fHist_znAx_V0_VxVy[i] = new TProfile2D(Form("fHist_znAx_V0_Run%d_Vz%d",runNums[i],1),"",NbinVt,0,NbinVt,vzBin,VzCut[0],VzCut[1],"");
482  fHist_znAy_V0_VxVy[i] = new TProfile2D(Form("fHist_znAy_V0_Run%d_Vz%d",runNums[i],1),"",NbinVt,0,NbinVt,vzBin,VzCut[0],VzCut[1],"");
483  }
484  for(int i=0;i<frunflag;i++){
485  fListHistos->Add(fHist_ZDCn_A_XY[i]);
486  fListHistos->Add(fHist_ZDCn_C_XY[i]);
493  }
494  }//---------- recenter1 ------------
495 
496 
497  fAnalysisUtil = new AliAnalysisUtils();
498  fAnalysisUtil->SetUseMVPlpSelection(kTRUE);
499  fAnalysisUtil->SetUseOutOfBunchPileUp(kTRUE);
500 
501  PostData(1,fListHistos); //posting in slot #1,
502  AliDebug(2,Form("\n::UserCreateOutPutObject(). NbinVt= %d, frunflag= %d, dataset: %s, Analysis= %s\n\n",NbinVt,frunflag,fDataSet.Data(),fAnalysisSet.Data()));
503 
504 }
505 
506 //________________________________________________________________________
508 {
509 
510  int stepCount = 0;
511  AliDebug(2,Form("\n ... ::UserExec() is being called. Step %d... \n",stepCount));
512 
513  fHist_Event_count->Fill(stepCount);
514  stepCount++;
515 
516  /*if(!checkOnce){
517  checkOnce++;
518  }*/
519 
520  AliAODEvent *aod = dynamic_cast<AliAODEvent*>(InputEvent());
521  fEvent = dynamic_cast<AliFlowEventSimple*>(GetInputData(1));
522 
523 
524  if(!aod){
525  AliDebug(2,"\n ... ::UserExec no aod found..... \n");
526  return;
527  }
528 
529  fHist_Event_count->Fill(stepCount);
530  stepCount++;
531 
532 
533  AliAODVertex *pVertex = aod->GetPrimaryVertex();
534  Double_t Vxyz[3] = {0,0,0};
535  Vxyz[0] = pVertex->GetX();
536  Vxyz[1] = pVertex->GetY();
537  Vxyz[2] = pVertex->GetZ();
538 
539  //------- Apply Necessary event cuts ---------
540  if(Vxyz[2] >= VzCut[1] || Vxyz[2] <= VzCut[0]) return;
541  fHist_Event_count->Fill(stepCount);
542  stepCount++;
543 
544  if(Vxyz[0] >= VxCut[1] || Vxyz[0] <= VxCut[0]) return;
545  fHist_Event_count->Fill(stepCount);
546  stepCount++;
547 
548  if(Vxyz[1] >= VyCut[1] || Vxyz[1] <= VyCut[0]) return;
549  fHist_Event_count->Fill(stepCount);
550  stepCount++;
551 
552 
553  Double_t EvtCent = fEvent->GetCentrality();
554  fHist_Cent_count1-> Fill(EvtCent);
555 
556 
557  //---------- get runindex: --------------
558  Int_t runNumber = aod->GetRunNumber();
559  Int_t runindex = -111;
560 
561  for(int i=0;i<frunflag;i++){
562  if(runNumber==runNums[i])
563  {
564  runindex = i;
565  break;
566  }
567  }
568  if(runindex<0) {
569  AliDebug(2,"\n ::UserExec Runnumber is not listed ..... \n");
570  exit(1);
571  }
572  //-----------------------------------------
573 
574  //--------- starting pileup rejection work: --------
575  Double_t centrV0M=300; Double_t centrCL1=300;
576  Double_t centrCL0=300; Double_t centrTRK=300;
577 
578  if(fDataSet=="2010"||fDataSet=="2011"){
579  centrV0M = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("V0M");
580  centrCL1 = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("CL1");
581  centrCL0 = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("CL0");
582  centrTRK = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("TRK");
583  }// 2010/2011
584 
585  else{
586  fMultSelection = (AliMultSelection*) InputEvent()->FindListObject("MultSelection");
587  if(!fMultSelection) {
588  //If you get this warning, please check if AliMultSelectionTask actually ran (before your task)
589  AliDebug(2,Form("\n **WARNING** ::UserExec() AliMultSelection object not found. Step# %d\n",stepCount));
590  return;
591  }
592  centrV0M = fMultSelection->GetMultiplicityPercentile("V0M");
593  centrCL1 = fMultSelection->GetMultiplicityPercentile("CL1");
594  centrCL0 = fMultSelection->GetMultiplicityPercentile("CL0");
595  centrTRK = fMultSelection->GetMultiplicityPercentile("TRK");
596  }// 2015
597 
598 
599  Bool_t BisPileup=kFALSE;
600 
601  if(fRejectPileUp && InputEvent()) {
602  //if(!fCutsEvent->IsSelected(InputEvent(),MCEvent())) return;
603  //if(fDataSet!=k2015 && fDataSet!=k2015v6) { //jacopo
604  if(fDataSet!="2015") {
605  if(plpMV(aod)) {
606  fPileUpCount->Fill(0.5);
607  BisPileup=kTRUE;
608  }
609  Int_t isPileup = aod->IsPileupFromSPD(3);
610  if(isPileup != 0) {
611  fPileUpCount->Fill(1.5);
612  BisPileup=kTRUE;
613  }
614  if(((AliAODHeader*)aod->GetHeader())->GetRefMultiplicityComb08() < 0) {
615  fPileUpCount->Fill(2.5);
616  BisPileup=kTRUE;
617  }
618  if(aod->IsIncompleteDAQ()) {
619  fPileUpCount->Fill(3.5);
620  BisPileup=kTRUE;
621  }
622 
623  //check vertex consistency
624  const AliAODVertex* vtTrc = aod->GetPrimaryVertex();
625  const AliAODVertex* vtSPD = aod->GetPrimaryVertexSPD();
626 
627  if(vtTrc->GetNContributors() < 2 || vtSPD->GetNContributors()<1) {
628  fPileUpCount->Fill(5.5);
629  BisPileup=kTRUE;
630  }
631 
632  double covTrc[6], covSPD[6];
633  vtTrc->GetCovarianceMatrix(covTrc);
634  vtSPD->GetCovarianceMatrix(covSPD);
635 
636  double dz = vtTrc->GetZ() - vtSPD->GetZ();
637 
638  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
639  double errTrc = TMath::Sqrt(covTrc[5]);
640  double nsigTot = dz/errTot;
641  double nsigTrc = dz/errTrc;
642 
643  if(TMath::Abs(dz)>0.2 || TMath::Abs(nsigTot)>10 || TMath::Abs(nsigTrc)>20) {
644  fPileUpCount->Fill(6.5);
645  BisPileup=kTRUE;
646  }
647  if(fAnalysisUtil->IsPileUpEvent(InputEvent())) {
648  fPileUpCount->Fill(7.5);
649  BisPileup=kTRUE;
650  }
651  }
652 
653  else {
654  // pileup from AliMultSelection
655  if(!fMultSelection->GetThisEventIsNotPileup())
656  fPileUpMultSelCount->Fill(0.5);
657  if(!fMultSelection->GetThisEventIsNotPileupMV())
658  fPileUpMultSelCount->Fill(1.5);
659  if(!fMultSelection->GetThisEventIsNotPileupInMultBins())
660  fPileUpMultSelCount->Fill(2.5);
661  if(!fMultSelection->GetThisEventHasNoInconsistentVertices())
662  fPileUpMultSelCount->Fill(3.5);
663  if(!fMultSelection->GetThisEventPassesTrackletVsCluster())
664  fPileUpMultSelCount->Fill(4.5);
665  if(!fMultSelection->GetThisEventIsNotAsymmetricInVZERO())
666  fPileUpMultSelCount->Fill(5.5);
667  if(!fMultSelection->GetThisEventIsNotIncompleteDAQ())
668  fPileUpMultSelCount->Fill(6.5);
669  if(!fMultSelection->GetThisEventHasGoodVertex2016())
670  fPileUpMultSelCount->Fill(7.5);
671 
672  BisPileup=kFALSE;
673 
674  // pile-up a la Dobrin for LHC15o
675  if(plpMV(aod)) {
676  fPileUpCount->Fill(0.5);
677  BisPileup=kTRUE;
678  }
679 
680  Int_t isPileup = aod->IsPileupFromSPD(3);
681  if(isPileup != 0) {
682  fPileUpCount->Fill(1.5);
683  BisPileup=kTRUE;
684  }
685 
686  if(((AliAODHeader*)aod->GetHeader())->GetRefMultiplicityComb08() < 0) {
687  fPileUpCount->Fill(2.5);
688  BisPileup=kTRUE;
689  }
690 
691  if(aod->IsIncompleteDAQ()) {
692  fPileUpCount->Fill(3.5);
693  BisPileup=kTRUE;
694  }
695 
696  if(fabs(centrV0M-centrCL1)>7.5) {
697  fPileUpCount->Fill(4.5);
698  BisPileup=kTRUE;
699  }
700 
701  // check vertex consistency
702  const AliAODVertex* vtTrc = aod->GetPrimaryVertex();
703  const AliAODVertex* vtSPD = aod->GetPrimaryVertexSPD();
704 
705  if (vtTrc->GetNContributors() < 2 || vtSPD->GetNContributors()<1) {
706  fPileUpCount->Fill(5.5);
707  BisPileup=kTRUE;
708  }
709 
710  double covTrc[6], covSPD[6];
711  vtTrc->GetCovarianceMatrix(covTrc);
712  vtSPD->GetCovarianceMatrix(covSPD);
713 
714  double dz = vtTrc->GetZ() - vtSPD->GetZ();
715 
716  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
717  double errTrc = TMath::Sqrt(covTrc[5]);
718  double nsigTot = dz/errTot;
719  double nsigTrc = dz/errTrc;
720 
721  if(TMath::Abs(dz)>0.2 || TMath::Abs(nsigTot)>10 || TMath::Abs(nsigTrc)>20) {
722  fPileUpCount->Fill(6.5);
723  BisPileup=kTRUE;
724  }
725 
726  // cuts on tracks
727  const Int_t nTracks = aod->GetNumberOfTracks();
728  Int_t multEsd = ((AliAODHeader*)aod->GetHeader())->GetNumberOfESDTracks();
729 
730  Int_t multTrk = 0;
731  Int_t multTrkBefC = 0;
732  Int_t multTrkTOFBefC = 0;
733  Int_t multTPC = 0;
734 
735  for(Int_t it = 0; it < nTracks; it++) {
736  AliAODTrack* aodTrk = (AliAODTrack*)aod->GetTrack(it);
737  if(!aodTrk) {
738  delete aodTrk;
739  continue;
740  }
741 // if(aodTrk->TestFilterBit(32)){
742 // multTrkBefC++;
743 // if(TMath::Abs(aodTrk->GetTOFsignalDz()) <= 10. && aodTrk->GetTOFsignal() >= 12000. && aodTrk->GetTOFsignal() <= 25000.)
744 // multTrkTOFBefC++;
745 // if((TMath::Abs(aodTrk->Eta()) < 0.8) && (aodTrk->GetTPCNcls() >= 70) && (aodTrk->Pt() >= 0.2) && (aodTrk->Pt() < 20.))
746 // multTrk++;
747 // }
748  if(aodTrk->TestFilterBit(128))
749  multTPC++;
750  } // end of for (Int_t it = 0; it < nTracks; it++)
751 
752  Double_t multTPCn = multTPC;
753  Double_t multEsdn = multEsd;
754  Double_t multESDTPCDif = multEsdn - multTPCn*3.38;
755 
756  if(multESDTPCDif > (fRejectPileUpTight?700.:15000.)) {
757  fPileUpCount->Fill(7.5);
758  BisPileup=kTRUE;
759  }
760 
761  if(fRejectPileUpTight) {
762  if(BisPileup==kFALSE) {
763  if(!fMultSelection->GetThisEventIsNotPileup()) BisPileup=kTRUE;
764  if(!fMultSelection->GetThisEventIsNotPileupMV()) BisPileup=kTRUE;
765  if(!fMultSelection->GetThisEventIsNotPileupInMultBins()) BisPileup=kTRUE;
766  if(!fMultSelection->GetThisEventHasNoInconsistentVertices()) BisPileup=kTRUE;
767  if(!fMultSelection->GetThisEventPassesTrackletVsCluster()) BisPileup=kTRUE;
768  if(!fMultSelection->GetThisEventIsNotIncompleteDAQ()) BisPileup=kTRUE;
769  if(!fMultSelection->GetThisEventHasGoodVertex2016()) BisPileup=kTRUE;
770  if(BisPileup) fPileUpCount->Fill(8.5);
771  }
772  }
773  }//-> 2015
774  }
775 
776  //----------- pile up rejection done ---------
777 
778 
779 
780 
781 
782 
783 
784  if(BisPileup){
785  AliDebug(2,Form("\n ::UserExec Pileup event found, skipping.\n Dataset: %s \n",fDataSet.Data()));
786  return;
787  }
788 
789  fHist_Event_count->Fill(stepCount);
790  stepCount++;
791 
792 
793 
794 
795 
796 
797  // ********** fZDCgain alpha = 0.50 instead of 0.35 *********
798  AliAODZDC *aodZDC = aod->GetZDCData();
799  Float_t fZDCGainAlpha = 0.500;
800  Float_t energyZNC = (Float_t) (aodZDC->GetZNCEnergy());
801  Float_t energyZPC = (Float_t) (aodZDC->GetZPCEnergy());
802  Float_t energyZNA = (Float_t) (aodZDC->GetZNAEnergy());
803  Float_t energyZPA = (Float_t) (aodZDC->GetZPAEnergy());
804  Float_t energyZEM1 = (Float_t) (aodZDC->GetZEM1Energy());
805  Float_t energyZEM2 = (Float_t) (aodZDC->GetZEM2Energy());
806 
807  const Double_t * towZNC = aodZDC->GetZNCTowerEnergy();
808  const Double_t * towZPC = aodZDC->GetZPCTowerEnergy();
809  const Double_t * towZNA = aodZDC->GetZNATowerEnergy();
810  const Double_t * towZPA = aodZDC->GetZPATowerEnergy();
811 
812  const Double_t * towZNClg = aodZDC->GetZNCTowerEnergyLR(); // Low gain something, should not be used.
813  const Double_t * towZNAlg = aodZDC->GetZNATowerEnergyLR();
814 
815  Double_t towZPClg[5] = {0.,};
816  Double_t towZPAlg[5] = {0.,};
817 
818  for(Int_t it=0; it<5; it++) {
819  towZPClg[it] = 8*towZPC[it];
820  towZPAlg[it] = 8*towZNA[it];
821  }
822 
823  //sanity: remove if any of ZDC_C_A has negetive Energy:
824  if(towZNC[1]<0 || towZNC[2]<0 || towZNC[3]<0 || towZNC[4]<0) return;
825  fHist_Event_count->Fill(stepCount);
826  stepCount++;
827 
828  if(towZNA[1]<0 || towZNA[2]<0 || towZNA[3]<0 || towZNA[4]<0) return;
829  fHist_Event_count->Fill(stepCount);
830  stepCount++;
831 
832  if(fAnalysisSet=="recenter1"){
833  for(Int_t it=0; it<5; it++) {
834  fHist_ZDCA_En_Run[runindex]->Fill(EvtCent,it,towZNA[it]);
835  fHist_ZDCC_En_Run[runindex]->Fill(EvtCent,it,towZNC[it]);
836  }
837  }
838 
839 //********** Get centroid from ZDCs **************
840 
841  Double_t xyZNC[2]={999.,999.};
842  Double_t xyZNA[2]={999.,999.};
843 
844  Float_t zncEnergy=0., znaEnergy=0.;
845 
846  for(Int_t i=0; i<5; i++){
847  zncEnergy += towZNC[i];
848  znaEnergy += towZNA[i];
849  }
850 
851  //fEvent->SetZNCEnergy(towZNC[0]);
852  //fEvent->SetZNAEnergy(towZNA[0]);
853 
854  Double_t AvTowerGain[8] = {1., 1., 1., 1., 1., 1., 1., 1.};
855 
856 /*---------------------------------------------------------------
857  Int_t CenBin = GetCenBin(centrperc);
858  Double_t zvtxpos[3]={0.,0.,0.};
859  fFlowEvent->GetVertexPosition(zvtxpos);
860  Int_t RunNum=fFlowEvent->GetRun();
861  if(fTowerEqList) {
862  if(RunNum!=fCachedRunNum) {
863  for(Int_t i=0; i<8; i++) {
864  fTowerGainEq[i] = (TH1D*)(fTowerEqList->FindObject(Form("Run %d",RunNum))->FindObject(Form("fhnTowerGainEqFactor[%d][%d]",RunNum,i)));
865  }
866  }
867  }
868  Bool_t fUseMCCen = kFALSE; //rihan:hardcoded
869  if (fUseMCCen) {
870  if(aod->GetRunNumber() < 209122)
871  aodZDC->GetZNCentroidInPbPb(1380., xyZNC, xyZNA);
872  else
873  aodZDC->GetZNCentroidInPbPb(2510., xyZNC, xyZNA);
874  }
875  else {
876  //set tower gain equalization, if available
877  if(fTowerEqList) {
878  for(Int_t i=0; i<8; i++)
879  {
880  if(fTowerGainEq[i])
881  AvTowerGain[i] = fTowerGainEq[i]->GetBinContent(fTowerGainEq[i]->FindBin(centrperc));
882  }
883  }//--------------------------------------------------------------- */
884 
885 
886 
887  const Float_t x[4] = {-1.75, 1.75,-1.75, 1.75};
888  const Float_t y[4] = {-1.75, -1.75, 1.75, 1.75};
889 
890  Float_t numXZNC=0., numYZNC=0., denZNC=0., wZNC;
891  Float_t numXZNA=0., numYZNA=0., denZNA=0., wZNA;
892 
893  for(Int_t i=0; i<4; i++)
894  {
895  if(towZNC[i+1]>0.)
896  {
897  wZNC = TMath::Power(towZNC[i+1], fZDCGainAlpha)*AvTowerGain[i];
898  numXZNC += x[i]*wZNC;
899  numYZNC += y[i]*wZNC;
900  denZNC += wZNC;
901  }
902 
903  if(towZNA[i+1]>0.) {
904  wZNA = TMath::Power(towZNA[i+1], fZDCGainAlpha)*AvTowerGain[i+4];
905  numXZNA += x[i]*wZNA;
906  numYZNA += y[i]*wZNA;
907  denZNA += wZNA;
908  }
909  }
910 
911  if(denZNC!=0) {
912  xyZNC[0] = numXZNC/denZNC;
913  xyZNC[1] = numYZNC/denZNC;
914  }
915  else{
916  xyZNC[0] = 999.;
917  xyZNC[1] = 999.;
918  zncEnergy = 0.;
919  }
920  if(denZNA!=0) {
921  xyZNA[0] = numXZNA/denZNA;
922  xyZNA[1] = numYZNA/denZNA;
923  }
924  else{
925  xyZNA[0] = 999.;
926  xyZNA[1] = 999.;
927  znaEnergy = 0.;
928  }
929 
930 
931  //----- Important: zdcA_X = -zdcA_X ---------
932  xyZNA[0] = -1.*xyZNA[0];
933  //===========================================
934 
935 
936  if(sqrt(xyZNC[0]*xyZNC[0] + xyZNC[1]*xyZNC[1])>1.5) return;
937  fHist_Event_count->Fill(stepCount);
938  stepCount++;
939 
940  if(sqrt(xyZNA[0]*xyZNA[0] + xyZNA[1]*xyZNA[1])>1.5) return;
941  fHist_Event_count->Fill(stepCount);
942  stepCount++;
943 
944  // ---------------- ZDC multiplicities posted at the ende of this code-----------------------
945 
946 
947 
948 
949 
950 
951  fHist_EventperRun->Fill(runindex);
952 
953  if(fAnalysisSet=="recenter1"){
954  fHist_ZDCn_A_XY[runindex]->Fill(xyZNA[0],xyZNA[1]);
955  fHist_ZDCn_C_XY[runindex]->Fill(xyZNC[0],xyZNC[1]);
956  }
957 
958 
959  //----- calculate RefMult for event ---------
960  AliFlowTrackSimple* pTrack = NULL;
961  Int_t iTracks = fEvent->NumberOfTracks();
962 
963  Double_t Qnx_TPC[3] = {0,};
964  Double_t Qny_TPC[3] = {0,};
965  Double_t psi2,dPhi,dPt,dEta;
966 
967  Int_t nRefMult = 0;
968  Int_t npoiMult = 0;
969 
970 /*
971  for(Int_t i=0; i<iTracks; i++)
972  {
973  pTrack = fEvent->GetTrack(i);
974  if (!pTrack) continue;
975 
976  dPhi = pTrack->Phi();
977  dPt = pTrack-> Pt();
978  dEta = pTrack->Eta();
979  //dCharge = pTrack->Charge();
980 
981  if(dPt<0.2 || dPt>10.0) continue;
982  if(fabs(dEta)>0.8) continue;
983 
984  //fHistQA_etaphi->Fill(dPhi,dEta);
985 
986  Qnx_TPC[0] += TMath::Cos(2.*dPhi);
987  Qny_TPC[0] += TMath::Sin(2.*dPhi);
988 
989  npoiMult++;
990 
991 
992  if(fabs(dEta)<0.5)
993  nRefMult++;
994  }//track loop ends
995 */
996 
997  //double psi2 = 0.5*TMath::ATan2(Qny_TPC[0],Qnx_TPC[0]);
998  //if(psi2<0) psi2 += TMath::Pi();
999  //fHist_EventPlane2->Fill(psi2);
1000  //fill q vectors for TPC recentering
1001  //fHist_QnxRecent->Fill(EvtCent,runindex,(Qnx_TPC[0]/(double)npoiMult));
1002  //fHist_QnyRecent->Fill(EvtCent,runindex,(Qny_TPC[0]/(double)npoiMult));
1003  //---------------------------------------------------------
1004 
1005 
1006  psi2 = TMath::ATan2(xyZNC[1],xyZNC[0]);
1007  if(psi2<0) psi2 += 2.*TMath::Pi();
1008  fHist_Psi1_zdnC->Fill(psi2);
1009 
1010  psi2 = TMath::ATan2(xyZNA[1],xyZNA[0]);
1011  if(psi2<0) psi2 += 2.*TMath::Pi();
1012  fHist_Psi1_zdnA->Fill(psi2);
1013 
1014  fHist_Cent_count2-> Fill(EvtCent);
1015 
1016  Int_t nTracks = aod->GetNumberOfTracks(); //number of AOD tracks
1017 
1018 
1019 
1020  fHist_Vertex_Vz->Fill(Vxyz[2]);
1021  fHist_Vertex_XY->Fill(Vxyz[0],Vxyz[1]);
1022 
1023  fHist_Vx_vs_runnum->Fill(runindex,Vxyz[0]);
1024  fHist_Vy_vs_runnum->Fill(runindex,Vxyz[1]);
1025  fHist_Vz_vs_runnum->Fill(runindex,Vxyz[2]);
1026 
1027  fHist_tracks_vs_runnum->Fill(runindex,nTracks);
1028 
1029  Int_t indexVx = fHist_Vx_ArrayFinder->FindBin(Vxyz[0]);
1030  Int_t indexVy = fHist_Vy_ArrayFinder->FindBin(Vxyz[1]);
1031  Int_t indexVz = fHist_Vz_ArrayFinder->FindBin(Vxyz[2]);
1032 
1033  Double_t tVertexBin1 = (Double_t) (indexVy-1)*vxBin + indexVx - 0.5 ; // 'Int_t' because bin starts with 1.
1034 
1035 
1036  if(fAnalysisSet=="recenter1"){
1037  fHist_vBincount->Fill(tVertexBin1);
1038  fHist_znCx_V0_VxVy[runindex]->Fill(tVertexBin1,Vxyz[2],xyZNC[0]);
1039  fHist_znCy_V0_VxVy[runindex]->Fill(tVertexBin1,Vxyz[2],xyZNC[1]);
1040  fHist_znAx_V0_VxVy[runindex]->Fill(tVertexBin1,Vxyz[2],xyZNA[0]);
1041  fHist_znAy_V0_VxVy[runindex]->Fill(tVertexBin1,Vxyz[2],xyZNA[1]);
1042  }
1043 
1044 
1045  Double_t meanCx, meanCy, meanAx, meanAy;
1046  meanCx = 0; meanCy = 0; meanAx = 0; meanAy = 0;
1047 
1048  Int_t tVertexBin2 = (indexVy-1)*vxBin + indexVx; // 'Int_t' because bin starts with 1.
1049 
1050 
1051  if(fAnalysisSet=="recenter2"){
1052  Double_t tVertexBinD = (Double_t)tVertexBin2 - 0.5;
1053  fHist_vBincount->Fill(tVertexBinD);
1054 
1055  meanCx = fHist_znCx_V0_VxVy[runindex]->GetBinContent(tVertexBin2,indexVz);
1056  meanCy = fHist_znCy_V0_VxVy[runindex]->GetBinContent(tVertexBin2,indexVz);
1057  meanAx = fHist_znAx_V0_VxVy[runindex]->GetBinContent(tVertexBin2,indexVz);
1058  meanAy = fHist_znAy_V0_VxVy[runindex]->GetBinContent(tVertexBin2,indexVz);
1059 
1060  xyZNC[0] = xyZNC[0] - meanCx;
1061  xyZNC[1] = xyZNC[1] - meanCy;
1062  xyZNA[0] = xyZNA[0] - meanAx;
1063  xyZNA[1] = xyZNA[1] - meanAy;
1064 
1065  fHist_znCx_V0_Cent[runindex]->Fill(EvtCent,xyZNC[0]);
1066  fHist_znCy_V0_Cent[runindex]->Fill(EvtCent,xyZNC[1]);
1067  fHist_znAx_V0_Cent[runindex]->Fill(EvtCent,xyZNA[0]);
1068  fHist_znAy_V0_Cent[runindex]->Fill(EvtCent,xyZNA[1]);
1069  }
1070 
1071 
1072 
1073 
1074  if(fievent%100==0){
1075  AliDebug(2,Form("cent= %f run= %d <Cx> = %f \t <Ax>= %f \t <Cy>= %f \t<Ay>= %f",EvtCent,runNumber,meanCx,meanAx,meanCy,meanAy));
1076  }
1077  //std::cout<<fievent<<" cent = "<<EvtCent<<" run = "<<runNumber<<" <Cx> = "<<meanCx<<"\t <Ax> "<<meanAx<<"\t <Cy> "<<meanCy<<"\t <Ay> "<<meanAy<<std::endl;
1078 
1079 
1080 
1081 
1082  PostData(1,fListHistos);
1083  fievent++;
1084 
1085 }
1086 
1087 
1088 
1089 
1091 {
1092  // Called once at the end of the query
1093  /* AliFlowAnalysisIDCSP *fSPTerm = new AliFlowAnalysisIDCSP();
1094  fListHistos = (TList*) GetOutputData(1);
1095  if(fListHistos)
1096  {
1097  fSPTerm->GetOutputHistograms(fListHistos);
1098  fSPTerm->Finish();
1099  PostData(1,fListHistos);
1100  }
1101  else
1102  {
1103  std::cout << "histgram list pointer is empty in Scalar Product" << endl;
1104  } */
1105  AliDebug(2,"\n ... AliAnalysisTaskVnZDC::Terminate() is being called ... \n");
1106 }
1107 
1108 
1109 double AliAnalysisTaskVnZDC::GetWDist(const AliVVertex* v0, const AliVVertex* v1)
1110 {
1111  // calculate sqrt of weighted distance to other vertex
1112  if (!v0 || !v1) {
1113  AliDebug(2,"\n\n ::GetWDist => One of vertices is not valid\n\n");
1114  return 0;
1115  }
1116  static TMatrixDSym vVb(3);
1117  double dist = -1;
1118  double dx = v0->GetX()-v1->GetX();
1119  double dy = v0->GetY()-v1->GetY();
1120  double dz = v0->GetZ()-v1->GetZ();
1121  double cov0[6],cov1[6];
1122  v0->GetCovarianceMatrix(cov0);
1123  v1->GetCovarianceMatrix(cov1);
1124  vVb(0,0) = cov0[0]+cov1[0];
1125  vVb(1,1) = cov0[2]+cov1[2];
1126  vVb(2,2) = cov0[5]+cov1[5];
1127  vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
1128  vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
1129  vVb.InvertFast();
1130  if (!vVb.IsValid()) {
1131  AliDebug(2,"Singular Matrix\n");
1132  return dist;}
1133  dist = vVb(0,0)*dx*dx + vVb(1,1)*dy*dy + vVb(2,2)*dz*dz
1134  + 2*vVb(0,1)*dx*dy + 2*vVb(0,2)*dx*dz + 2*vVb(1,2)*dy*dz;
1135  return dist>0 ? TMath::Sqrt(dist) : -1;
1136 }
1137 
1139  { // check for multi-vertexer pile-up
1140  const int kMinPlpContrib = 5;
1141  const double kMaxPlpChi2 = 5.0;
1142  const double kMinWDist = 15;
1143 
1144  const AliVVertex* vtPrm = 0;
1145  const AliVVertex* vtPlp = 0;
1146 
1147  int nPlp = 0;
1148 
1149  if(!(nPlp=aod->GetNumberOfPileupVerticesTracks()))
1150  return kFALSE;
1151 
1152  vtPrm = aod->GetPrimaryVertex();
1153  if(vtPrm == aod->GetPrimaryVertexSPD())
1154  return kTRUE; // there are pile-up vertices but no primary
1155 
1156  //int bcPrim = vtPrm->GetBC();
1157 
1158  for(int ipl=0;ipl<nPlp;ipl++) {
1159  vtPlp = (const AliVVertex*)aod->GetPileupVertexTracks(ipl);
1160  if (vtPlp->GetNContributors() < kMinPlpContrib) continue;
1161  if (vtPlp->GetChi2perNDF() > kMaxPlpChi2) continue;
1162  //int bcPlp = vtPlp->GetBC();
1163  //if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2)
1164  // return kTRUE; // pile-up from other BC
1165 
1166  double wDst = GetWDist(vtPrm,vtPlp);
1167  if (wDst<kMinWDist) continue;
1168 
1169  return kTRUE; // pile-up: well separated vertices
1170  }
1171 
1172  return kFALSE;
1173  }
1174 
1175 
1176 
1177 
1178 
1179 
1180 
1181 
1182 
1183 
1184 /*// ---------------- ZDC multiplicities -----------------------
1185  Float_t MulA=0., MulC=0.;
1186  for(Int_t i=0; i<4; i++) {
1187  if(towZNC[i+1]>0.) {
1188  MulC += TMath::Power(towZNC[i+1], fZDCGainAlpha)*AvTowerGain[i+4];
1189  }
1190  if(towZNA[i+1]>0.) {
1191  MulA += TMath::Power(towZNA[i+1], fZDCGainAlpha)*AvTowerGain[i+4];
1192  }
1193  }
1194 //fhZNCcentroid->Fill(xyZNC[0], xyZNC[1]);
1195 //fhZNAcentroid->Fill(xyZNA[0], xyZNA[1]);
1196  Double_t testC[2] = {1,0};
1197  Double_t testA[2] = {1,0};
1198 //fEvent->SetZDC2Qsub(xyZNC,MulC,xyZNA,MulA);
1199  Float_t tdcSum = aodZDC->GetZDCTimeSum();
1200  Float_t tdcDiff = aodZDC->GetZDCTimeDiff();
1201  Double_t asymmetry = -999.;
1202  if((energyZNC+energyZNA)>0.)
1203  asymmetry = (energyZNC-energyZNA)/(energyZNC+energyZNA);
1204 //----------------------- ZDC Mutliplicities ----------------------- */
TProfile * fHist_znAx_V0_Cent[89]
ClassImp(AliAnalysisTaskVnZDC) AliAnalysisTaskVnZDC
TProfile2D * fHist_ZDCC_En_Run[89]
AliFlowEventSimple * fEvent
double Double_t
Definition: External.C:58
Definition: External.C:236
AliMultSelection * fMultSelection
TProfile * fHist_znCy_V0_Cent[89]
TProfile2D * fHist_znCy_V0_VxVy[89]
TProfile2D * fHist_znCx_V0_VxVy[89]
virtual void Terminate(Option_t *)
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
Bool_t plpMV(const AliAODEvent *aod)
TProfile2D * fHist_ZDCA_En_Run[89]
virtual void UserCreateOutputObjects()
TProfile * fHist_znAy_V0_Cent[89]
TProfile * fHist_znCx_V0_Cent[89]
Double_t GetCentrality() const
AliAnalysisUtils * fAnalysisUtil
TProfile2D * fHist_znAx_V0_VxVy[89]
double GetWDist(const AliVVertex *v0, const AliVVertex *v1)
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
virtual void UserExec(Option_t *option)
TProfile2D * fHist_znAy_V0_VxVy[89]
Int_t NumberOfTracks() const