AliPhysics  1168478 (1168478)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskZDCGainEq.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 //
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 "TProfile2D.h"
30 #include "TList.h"
31 #include "TH1.h"
32 #include "TH2.h"
33 #include "TH3.h"
34 #include "AliMultSelection.h"
35 #include "AliVVertex.h"
36 
37 #include "AliAnalysisManager.h"
38 #include "AliFlowEvent.h"
39 #include "AliFlowEventSimple.h"
40 #include "AliFlowCommonHist.h"
42 
43 #include "AliAnalysisManager.h"
44 #include "AliInputEventHandler.h"
45 #include "AliVEvent.h"
46 #include "AliESD.h"
47 #include "AliESDEvent.h"
48 #include "AliESDHeader.h"
49 #include "AliESDInputHandler.h"
50 #include "AliESDZDC.h"
51 #include "AliMultiplicity.h"
52 #include "AliAnalysisUtils.h"
53 #include "AliAODHandler.h"
54 #include "AliAODTrack.h"
55 #include "AliAODEvent.h"
56 #include "AliAODHeader.h"
57 #include "AliAODVertex.h"
58 #include "AliAODVZERO.h"
59 #include "AliAODZDC.h"
60 #include "AliAODMCHeader.h"
61 #include "AliMCEventHandler.h"
62 #include "AliMCEvent.h"
63 #include "AliHeader.h"
64 #include "AliVParticle.h"
65 #include "AliStack.h"
66 #include "AliAODMCParticle.h"
67 #include "AliAnalysisTaskSE.h"
68 #include "AliGenEventHeader.h"
69 #include "AliPhysicsSelectionTask.h"
70 #include "AliPhysicsSelection.h"
71 #include "AliBackgroundSelection.h"
72 #include "AliTriggerAnalysis.h"
73 #include "AliCentrality.h"
74 #include "AliLog.h"
76 #include "AliFlowVector.h"
77 
78 using std::endl;
79 using std::cout;
80 
81 
83 
84 //________________________________________________________________________
86  AliAnalysisTaskSE(name),
87  fEvent(NULL),
88  fMultSelection(NULL),
89  fAnalysisUtil(NULL),
90  fListHistos(NULL),
91  fListZDCQxy(NULL),
92  fListZDCWgt(NULL),
93  fListDummy1(NULL),
94  fListHijing(NULL),
95  fListSubRun(NULL),
96  fRejectPileUpTight(kFALSE),
97  fRejectPileUp(kFALSE),
98  bFillCosSin(kFALSE),
99  bFillZDCQAon(kFALSE),
100  bRunAveragedQn(kFALSE),
101  bApplyRecent(kFALSE),
102  fHarmonic(2),
103  frunflag(0),
104  fievent(0),
105  vxBin(0),
106  vyBin(0),
107  vzBin(0),
108  fcheckOnce(0),
109  fOldRunNum(0),
110  fHist_Event_count(NULL),
111  fPileUpMultSelCount(NULL),
112  fPileUpCount(NULL),
113  fHist_ChanWgt_ZDCC(NULL),
114  fHist_ChanWgt_ZDCA(NULL),
115  fHist_Vx_ArrayFinder(NULL),
116  fHist_Vy_ArrayFinder(NULL),
117  fHist_Vz_ArrayFinder(NULL),
118  fHist_Task_config(NULL),
119  fHist_Cent_woZDCcut(NULL),
120  fHist_Cent_wiZDCcut(NULL),
121  fHist_CutParameters(NULL),
122  fHist_Psi1_ZDCC_wGainCorr(NULL),
123  fHist_Psi1_ZDCA_wGainCorr(NULL),
124  fHist_Psi1_ZDCC_wRectCorr(NULL),
125  fHist_Psi1_ZDCA_wRectCorr(NULL),
126  fHist_v2xV1_ZDN_Norm_All(NULL),
127  fHist_v2xV1_ZDN_Refm_All(NULL),
128  fHist_v2xV1_ZDN_Cent_All(NULL),
129  fHist_ZDN_resol_Norm_All(NULL),
130  fHist_ZDN_resol_Cent_All(NULL),
131  fHist_ZDN_resol_Refm_All(NULL),
132  fHist_Vx_vs_runnum(NULL),
133  fHist_Vy_vs_runnum(NULL),
134  fHist_Vz_vs_runnum(NULL),
135  fWeight_Cent(NULL),
136  fDataSet("2010"),
137  fAnalysisSet("DoGainEq"),
138  sCentEstimator("V0")
139 {
140  for(int i=0;i<90;i++){
141  runNums[i] = 0;
142  fHist_ZDCA_En_Run[i] = NULL;
143  fHist_ZDCC_En_Run[i] = NULL;
144 
145  for(int j=0;j<10;j++){
146  fHist_znCx_V0_VxVy[i][j] = NULL;
147  fHist_znCy_V0_VxVy[i][j] = NULL;
148  fHist_znAx_V0_VxVy[i][j] = NULL;
149  fHist_znAy_V0_VxVy[i][j] = NULL;
150  }
151  }
152 
153  for(int i=0;i<4;i++){
154  for(int j=0;j<5;j++){
155  fHist_Qx_vs_Obs_woCorr[i][j] = NULL;
156  fHist_XX_vs_Obs_woCorr[i][j] = NULL;
157  fHist_Qx_vs_Obs_wiCorr[i][j] = NULL;
158  fHist_XX_vs_Obs_wiCorr[i][j] = NULL;
159  }
160  }
161  for(int i=0;i<20;i++){
162  fHist_Recenter_ZDCCx[i] = NULL;
163  fHist_Recenter_ZDCCy[i] = NULL;
164  fHist_Recenter_ZDCAx[i] = NULL;
165  fHist_Recenter_ZDCAy[i] = NULL;
166  }
167  for(int i=0;i<20;i++){
168  fHist_ZDCC_En_CommonCh[i] = NULL;
169  fHist_ZDCA_En_CommonCh[i] = NULL;
170  }
171 
172  for(int i=0;i<2;i++){
173  VxCut[i] = 0;
174  VyCut[i] = 0;
175  VzCut[i] = 0;
176  }
177 
178  for(int i=0;i<10;i++){
179  fFB_Efficiency_Cent[i] = NULL;
180  fHist_v2xV1_ZDN_pTDiff_All[i] = NULL;
181 
182  for(int j=0;j<4;j++){
183  fHist_v1xV1_ZDN_EtaDiff[j][i] = NULL;
184  fHist_v1xV1_ZDN_pTDiff[j][i] = NULL;
185  }
186  }
187 
188  for(int i=0;i<4;i++){
189  fHist_v2xV1_ZDN_Norm_Sep[i] = NULL;
190  fHist_v2xV1_ZDN_Cent_Sep[i] = NULL;
191  }
192  for(int i=0;i<2;i++){
193  fHist_ZDN_resol_Norm_Sep[i] = NULL;
194  fHist_ZDN_resol_Cent_Sep[i] = NULL;
195  }
196 
197  DefineInput(1, AliFlowEventSimple::Class()); // Input slot #1 works with an AliFlowEventSimple
198  DefineInput(2, AliFlowEventSimple::Class()); // Input slot #2 for ZDC flow event
199 
200 
201  DefineOutput(1,TList::Class());
202  DefineOutput(2,TList::Class());
203 
204  //fDataSet="2010";
205  //fAnalysisSet="DoGainEq";
206  //sCentEstimator="V0";
207 
208  //fTotalQvector = new TString("QaQb"); // "QaQb" (means Qa+Qb), "Qa" or "Qb"
209 
210 }//-------------constructor-----------------
211 
212 //________________________________________________
215  fEvent(NULL),
216  fMultSelection(NULL),
217  fAnalysisUtil(NULL),
218  fListHistos(NULL),
219  fListZDCQxy(NULL),
220  fListZDCWgt(NULL),
221  fListDummy1(NULL),
222  fListHijing(NULL),
223  fListSubRun(NULL),
224  fRejectPileUpTight(kFALSE),
225  fRejectPileUp(kFALSE),
226  bFillCosSin(kFALSE),
227  bFillZDCQAon(kFALSE),
228  bRunAveragedQn(kFALSE),
229  bApplyRecent(kFALSE),
230  fHarmonic(2),
231  frunflag(0),
232  fievent(0),
233  vxBin(0),
234  vyBin(0),
235  vzBin(0),
236  fcheckOnce(0),
237  fOldRunNum(0),
238  fHist_Event_count(NULL),
239  fPileUpMultSelCount(NULL),
240  fPileUpCount(NULL),
241  fHist_ChanWgt_ZDCC(NULL),
242  fHist_ChanWgt_ZDCA(NULL),
243  fHist_Vx_ArrayFinder(NULL),
244  fHist_Vy_ArrayFinder(NULL),
245  fHist_Vz_ArrayFinder(NULL),
246  fHist_Task_config(NULL),
247  fHist_Cent_woZDCcut(NULL),
248  fHist_Cent_wiZDCcut(NULL),
249  fHist_CutParameters(NULL),
250  fHist_Psi1_ZDCC_wGainCorr(NULL),
251  fHist_Psi1_ZDCA_wGainCorr(NULL),
252  fHist_Psi1_ZDCC_wRectCorr(NULL),
253  fHist_Psi1_ZDCA_wRectCorr(NULL),
254  fHist_v2xV1_ZDN_Norm_All(NULL),
255  fHist_v2xV1_ZDN_Refm_All(NULL),
256  fHist_v2xV1_ZDN_Cent_All(NULL),
257  fHist_ZDN_resol_Norm_All(NULL),
258  fHist_ZDN_resol_Cent_All(NULL),
259  fHist_ZDN_resol_Refm_All(NULL),
260  fHist_Vx_vs_runnum(NULL),
261  fHist_Vy_vs_runnum(NULL),
262  fHist_Vz_vs_runnum(NULL),
263  fWeight_Cent(NULL),
264  fDataSet("2010"),
265  fAnalysisSet("DoGainEq"),
266  sCentEstimator("V0")
267 {
268  for(int i=0;i<90;i++){
269  runNums[i] = 0;
270  fHist_ZDCA_En_Run[i] = NULL;
271  fHist_ZDCC_En_Run[i] = NULL;
272 
273  for(int j=0;j<10;j++){
274  fHist_znCx_V0_VxVy[i][j] = NULL;
275  fHist_znCy_V0_VxVy[i][j] = NULL;
276  fHist_znAx_V0_VxVy[i][j] = NULL;
277  fHist_znAy_V0_VxVy[i][j] = NULL;
278  }
279  }
280 
281  for(int i=0;i<4;i++){
282  for(int j=0;j<5;j++){
283  fHist_Qx_vs_Obs_woCorr[i][j] = NULL;
284  fHist_XX_vs_Obs_woCorr[i][j] = NULL;
285  fHist_Qx_vs_Obs_wiCorr[i][j] = NULL;
286  fHist_XX_vs_Obs_wiCorr[i][j] = NULL;
287  }
288  }
289  for(int i=0;i<20;i++){
290  fHist_Recenter_ZDCCx[i] = NULL;
291  fHist_Recenter_ZDCCy[i] = NULL;
292  fHist_Recenter_ZDCAx[i] = NULL;
293  fHist_Recenter_ZDCAy[i] = NULL;
294  }
295  for(int i=0;i<20;i++){
296  fHist_ZDCC_En_CommonCh[i] = NULL;
297  fHist_ZDCA_En_CommonCh[i] = NULL;
298  }
299 
300  for(int i=0;i<2;i++){
301  VxCut[i] = 0;
302  VyCut[i] = 0;
303  VzCut[i] = 0;
304  }
305 
306  for(int i=0;i<10;i++){
307  fFB_Efficiency_Cent[i] = NULL;
308  fHist_v2xV1_ZDN_pTDiff_All[i] = NULL;
309 
310  for(int j=0;j<4;j++){
311  fHist_v1xV1_ZDN_EtaDiff[j][i] = NULL;
312  fHist_v1xV1_ZDN_pTDiff[j][i] = NULL;
313  }
314  }
315 
316  for(int i=0;i<4;i++){
317  fHist_v2xV1_ZDN_Norm_Sep[i] = NULL;
318  fHist_v2xV1_ZDN_Cent_Sep[i] = NULL;
319  }
320  for(int i=0;i<2;i++){
321  fHist_ZDN_resol_Norm_Sep[i] = NULL;
322  fHist_ZDN_resol_Cent_Sep[i] = NULL;
323  }
324 
325  //fDataSet="2010";
326  //fAnalysisSet="DoGainEq";
327  //sCentEstimator="V0";
328 }
329 
330 
331 
332 
333 
334 //________________________________________________________________________
336 {
337  delete fListHistos;
338  delete fListZDCQxy;
339  delete fListZDCWgt;
340  delete fListDummy1;
341  delete fListHijing;
342 
343  delete fMultSelection;
344  delete fAnalysisUtil; // it is '= new' !!!
345 
346 
347  //these histograms are not in any list:
348  //therefore, deleted manually.
349  //delete *fHist_ChanWgt_ZDCC; //can't delete, Execption thrown.!! why?
350  //delete *fHist_ChanWgt_ZDCA;
351  //delete *fHist_Vx_ArrayFinder;
352  //delete *fHist_Vy_ArrayFinder;
353  //delete *fHist_Vz_ArrayFinder;
354 
355 
356  //printf("\n\n ~AliAnalysisTaskZDCGainEq::Destructor is called..!!\n\n");
357 }
358 
359 //________________________________________________________________________
361 {
362 
363  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};
364 
365  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};
366 
367  Int_t runArray_2015[90] = {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, 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};
368 
369 
370  if(fDataSet=="2010"){
371  frunflag = 89;
372  for(int i=0;i<frunflag;i++)
373  runNums[i] = runArray_2010[i];
374  }
375  if(fDataSet=="2011"){
376  frunflag = 10; //<--------- 2011 runnumbers have to be checked and put...
377  for(int i=0;i<frunflag;i++)
378  runNums[i] = runArray_2011[i];
379  }
380  if(fDataSet=="2015"){
381  frunflag = 90;
382  for(int i=0;i<frunflag;i++)
383  runNums[i] = runArray_2015[i];
384  }
385 
386  fListHistos = new TList();
387  fListHistos->SetOwner(kTRUE);
388 
389  fHist_Event_count = new TH1F("fHist_Event_count"," ",20,0,20);
390  fHist_Event_count->GetXaxis()->SetBinLabel(1,"Called Exec()");
391  fHist_Event_count->GetXaxis()->SetBinLabel(2,"AOD Exist");
392  fHist_Event_count->GetXaxis()->SetBinLabel(3,"Pass VzCut");
393  fHist_Event_count->GetXaxis()->SetBinLabel(4,"Pass VxCut");
394  fHist_Event_count->GetXaxis()->SetBinLabel(5,"Pass VyCut");
395  fHist_Event_count->GetXaxis()->SetBinLabel(6,"NotPileUp");
396  fHist_Event_count->GetXaxis()->SetBinLabel(7,"ZNC Ei>=0");
397  fHist_Event_count->GetXaxis()->SetBinLabel(8,"ZNA Ei>=0");
398  fHist_Event_count->GetXaxis()->SetBinLabel(9,"|QnC| < 1.5");
399  fHist_Event_count->GetXaxis()->SetBinLabel(10,"|QnA| < 1.5");
400  fHist_Event_count->GetXaxis()->SetBinLabel(11,"#Psi_{A}=0 && #Psi_{C}=0");
401  fHist_Event_count->GetXaxis()->SetBinLabel(12,"..TBA..");
403 
404  fPileUpCount = new TH1F("fPileUpCount", "fPileUpCount", 9, 0., 9.);
405  fPileUpCount->GetXaxis()->SetBinLabel(1,"plpMV");
406  fPileUpCount->GetXaxis()->SetBinLabel(2,"fromSPD");
407  fPileUpCount->GetXaxis()->SetBinLabel(3,"RefMultiplicityComb08");
408  fPileUpCount->GetXaxis()->SetBinLabel(4,"IncompleteDAQ");
409  fPileUpCount->GetXaxis()->SetBinLabel(5,"abs(V0M-CL1)>7.5");
410  fPileUpCount->GetXaxis()->SetBinLabel(6,"missingVtx");
411  fPileUpCount->GetXaxis()->SetBinLabel(7,"inconsistentVtx");
412  fPileUpCount->GetXaxis()->SetBinLabel(8,"multESDTPCDif");
413  fPileUpCount->GetXaxis()->SetBinLabel(9,"extraPileUpMultSel");
415 
416  fPileUpMultSelCount = new TH1F("fPileUpMultSelCount", "fPileUpMultSelCount", 8, 0., 8.);
417  fPileUpMultSelCount->GetXaxis()->SetBinLabel(1,"IsNotPileup");
418  fPileUpMultSelCount->GetXaxis()->SetBinLabel(2,"IsNotPileupMV");
419  fPileUpMultSelCount->GetXaxis()->SetBinLabel(3,"IsNotPileupInMultBins");
420  fPileUpMultSelCount->GetXaxis()->SetBinLabel(4,"InconsistentVertices");
421  fPileUpMultSelCount->GetXaxis()->SetBinLabel(5,"TrackletVsCluster");
422  fPileUpMultSelCount->GetXaxis()->SetBinLabel(6,"AsymmetricInVZERO");
423  fPileUpMultSelCount->GetXaxis()->SetBinLabel(7,"IncompleteDAQ");
424  fPileUpMultSelCount->GetXaxis()->SetBinLabel(8,"GoodVertex2016");
426 
427  fHist_Task_config = new TH1F("fTask_Configuration", "Task Configuration", 20, 0., 20.);
428  fHist_Task_config->GetXaxis()->SetBinLabel(1,"IsZDCQAon");
429  fHist_Task_config->GetXaxis()->SetBinLabel(2,"IsCentV0M");
430  fHist_Task_config->GetXaxis()->SetBinLabel(3,"IsCentTPC");
431  fHist_Task_config->GetXaxis()->SetBinLabel(4,"IsCentCL1");
432  fHist_Task_config->GetXaxis()->SetBinLabel(5,"IsPileUpOn");
433  fHist_Task_config->GetXaxis()->SetBinLabel(6,"IsPileUpTightOn");
434  fHist_Task_config->GetXaxis()->SetBinLabel(7,"IsFillGainEq");
435  fHist_Task_config->GetXaxis()->SetBinLabel(8,"IsDoGainEq");
436  fHist_Task_config->GetXaxis()->SetBinLabel(9,"IsAvailGainFile");
437  fHist_Task_config->GetXaxis()->SetBinLabel(10,"IsFillCosSin");
438  fHist_Task_config->GetXaxis()->SetBinLabel(11,"IsDoRecenter");
439  fHist_Task_config->GetXaxis()->SetBinLabel(12,"IsAvailRecntFile");
440  fHist_Task_config->GetXaxis()->SetBinLabel(13,"IsRunByRun");
442 
443 
444  fHist_CutParameters = new TH1F("fTask_CutParameters", "Variable Cut Values", 10, 0., 10.);
445  fHist_CutParameters->GetXaxis()->SetBinLabel(1,"VzCutLowValue");
446  fHist_CutParameters->GetXaxis()->SetBinLabel(2,"VzCutHighValue");
447  fHist_CutParameters->GetXaxis()->SetBinLabel(3,"NBins_Vz");
448  fHist_CutParameters->GetXaxis()->SetBinLabel(4,"VxCutHighValue");
449  fHist_CutParameters->GetXaxis()->SetBinLabel(5,"VxCutLowValue");
450  fHist_CutParameters->GetXaxis()->SetBinLabel(6,"NBins_Vx");
451  fHist_CutParameters->GetXaxis()->SetBinLabel(7,"VyCutHighValue");
452  fHist_CutParameters->GetXaxis()->SetBinLabel(8,"VyCutLowValue");
453  fHist_CutParameters->GetXaxis()->SetBinLabel(9,"NBins_Vy");
455 
456  fHist_Cent_woZDCcut = new TH1F("fHist_Cent_before_ZDCcut"," ",100,0,100);
458  fHist_Cent_wiZDCcut = new TH1F("fHist_Cent_afterr_ZDCcut"," ",100,0,100);
460 
461 
462 
463 
464 
465  //-------- define and add the recenter histograms ----------------
466  if(fDataSet=="2010") {
467  //vxBin = 9;
468  //vyBin = 9;
469  //vzBin = 10;
470 
471  VxCut[0] = -0.030;
472  VxCut[1] = 0.015;
473  VyCut[0] = 0.150;
474  VyCut[1] = 0.204;
475  VzCut[0] = -10.0;
476  VzCut[1] = 9.4;
477  }
478 
479  if(fDataSet=="2015"){
480  //vxBin = 8;
481  //vyBin = 8;
482  //vzBin = 10;
483 
484  VxCut[0] = 0.060;
485  VxCut[1] = 0.086;
486  VyCut[0] = 0.321;
487  VyCut[1] = 0.345;
488  VzCut[0] = -10.0;
489  VzCut[1] = 10.0;
490  }
491 
492  if(fDataSet=="2011"){
493  AliDebug(2,"\n\n !!** WARNING ***!! \n cuts not defined for DATASET: %s\n ...EXIT...\n\n)");
494  exit(1);
495  }
496 
497  const int NbinVt = (vyBin-1)*vxBin + vxBin;
498 
499  fHist_Vx_ArrayFinder = new TH1F("fHist_Vx_ArrayFinder","",vxBin,VxCut[0],VxCut[1]);
500  fHist_Vy_ArrayFinder = new TH1F("fHist_Vy_ArrayFinder","",vyBin,VyCut[0],VyCut[1]);
501  fHist_Vz_ArrayFinder = new TH1F("fHist_Vz_ArrayFinder","",vzBin,VzCut[0],VzCut[1]);
502 
503 
504 
505  char name[100];
506 
507 
508  if(fAnalysisSet=="FillGainEq") {
509 
510  fHist_Vx_vs_runnum = new TProfile("fHist_Vx_vs_runnum","<Vx>_vs_runnum",frunflag,0,frunflag,"s");
512  fHist_Vy_vs_runnum = new TProfile("fHist_Vy_vs_runnum","<Vy>_vs_runnum",frunflag,0,frunflag,"s");
514  fHist_Vz_vs_runnum = new TProfile("fHist_Vz_vs_runnum","<Vz>_vs_runnum",frunflag,0,frunflag,"s");
516 
517  for(int i=0;i<vzBin;i++) {
518  fHist_ZDCC_En_CommonCh[i] = new TProfile2D(Form("fHist_ZDCC_En_CommonCh_Vz%d",i+1),"",100,0,100,NbinVt,0,NbinVt,"");
519  fHist_ZDCA_En_CommonCh[i] = new TProfile2D(Form("fHist_ZDCA_En_CommonCh_Vz%d",i+1),"",100,0,100,NbinVt,0,NbinVt,"");
522  }
523  for(int i=0;i<frunflag;i++) {
524  //store: ZDC energy for gain calibration:
525  fHist_ZDCA_En_Run[i] = new TProfile2D(Form("fHist_ZDCA_En_Run%d",runNums[i]),"",100,0,100,5,0,5,"");
527  fHist_ZDCC_En_Run[i] = new TProfile2D(Form("fHist_ZDCC_En_Run%d",runNums[i]),"",100,0,100,5,0,5,"");
529  }
530  }
531 
532 
533  if(bFillZDCQAon){
534 
535  fHist_Task_config->Fill(0.5);
536 
537  if(fAnalysisSet=="FillGainEq"){
538  fHist_Psi1_ZDCC_wGainCorr = new TH1F("fHist_Psi1_ZDCC_woGainCorr","", 200, 0,2.*TMath::Pi());
539  fHist_Psi1_ZDCA_wGainCorr = new TH1F("fHist_Psi1_ZDCA_woGainCorr","", 200, 0,2.*TMath::Pi());
540  }
541  else{
542  fHist_Psi1_ZDCC_wGainCorr = new TH1F("fHist_Psi1_ZDCC_wiGainCorr","", 200, 0,2.*TMath::Pi());
543  fHist_Psi1_ZDCA_wGainCorr = new TH1F("fHist_Psi1_ZDCA_wiGainCorr","", 200, 0,2.*TMath::Pi());
544  }
545 
548 
549  fHist_Psi1_ZDCC_wRectCorr = new TH1F("fHist_Psi1_ZDCC_wiRectCorr","", 200, 0,2.*TMath::Pi());
551  fHist_Psi1_ZDCA_wRectCorr = new TH1F("fHist_Psi1_ZDCA_wiRectCorr","", 200, 0,2.*TMath::Pi());
553 
554 
555  TString sNameQn[4] = {"Xa","Xc","Ya","Yc"};
556  TString sNameQn2[4] = {"XaXc","YaYc","XcYa","YcXa"};
557  TString sNameVs[5] = {"Cent","Mult","Vx","Vy","Vz"};
558  Int_t nBinNumVs[5] = {100, 100, 50, 50, 400}; // number of bins for "Cent", "Mult", "Vx","Vy","Vz"
559  Double_t lBinLowVs[5] = { 0, 0, VxCut[0], VyCut[0], -10}; // bin low value: "Cent", "Mult", "Vx", "Vy", "Vz"
560  Double_t lBinHighVs[5] = {100, 4000, VxCut[1], VyCut[1], 10}; // bin high value: "Cent", "Mult", "Vx", "Vy", "Vz"
561 
562 
563  for(int i=0;i<4;i++) {
564  for(int j=0;j<5;j++) {//fHist_Qx_vs_Obs_woCorr
565  //store: X,Y position for recenter:
566  sprintf(name,"fHist_%s_vs_%s_woCorr",static_cast<const char*>(sNameQn[i]),static_cast<const char*>(sNameVs[j]));
567  fHist_Qx_vs_Obs_woCorr[i][j] = new TProfile(name,"", nBinNumVs[j], lBinLowVs[j], lBinHighVs[j],"");
569 
570  sprintf(name,"fHist_%s_vs_%s_woCorr",static_cast<const char*>(sNameQn2[i]),static_cast<const char*>(sNameVs[j]));
571  fHist_XX_vs_Obs_woCorr[i][j] = new TProfile(name,"", nBinNumVs[j], lBinLowVs[j], lBinHighVs[j],"");
573 
574  if(bFillCosSin && fAnalysisSet=="FillGainEq") {
575  sprintf(name,"fHist_%s_vs_%s_woCorrCosSin",static_cast<const char*>(sNameQn[i]),static_cast<const char*>(sNameVs[j]));
576  }
577  else{
578  sprintf(name,"fHist_%s_vs_%s_wiCorr",static_cast<const char*>(sNameQn[i]),static_cast<const char*>(sNameVs[j]));
579  }
580  fHist_Qx_vs_Obs_wiCorr[i][j] = new TProfile(name,"", nBinNumVs[j], lBinLowVs[j], lBinHighVs[j],"");
582 
583  if(bFillCosSin && fAnalysisSet=="FillGainEq") {
584  sprintf(name,"fHist_%s_vs_%s_woCorrCosSin",static_cast<const char*>(sNameQn2[i]),static_cast<const char*>(sNameVs[j]));
585  }
586  else{
587  sprintf(name,"fHist_%s_vs_%s_wiCorr",static_cast<const char*>(sNameQn2[i]),static_cast<const char*>(sNameVs[j]));
588  }
589 
590  fHist_XX_vs_Obs_wiCorr[i][j] = new TProfile(name,"", nBinNumVs[j], lBinLowVs[j], lBinHighVs[j],"");
592  }
593  }
594  }
595 
596 
597  if(vzBin>20){
598  printf("\n\n::UserCreateOutPutObject()\n Vz Binning more than 20 not allowed\n ==> Exit <== \n\n");
599  exit(1);
600  }
601 
602  const int VzBinIter = (const int) vzBin;
603 
604 
605  //filling the Q vectors for recenter:
606  if(fAnalysisSet=="DoGainEq") {
607  if(!bApplyRecent) {
608 
609  fListZDCQxy = new TList();
610  fListZDCQxy->SetOwner(kTRUE);
611 
612  if(!bRunAveragedQn) {
613  for(int i=0;i<frunflag;i++) {
614  for(int j=0;j<VzBinIter;j++) {
615  fHist_znCx_V0_VxVy[i][j] = new TProfile2D(Form("fHist_znCx_V0_Run%d_Vz%d",runNums[i],j+1),"",NbinVt,0,NbinVt,90,0,90);
616  fHist_znCy_V0_VxVy[i][j] = new TProfile2D(Form("fHist_znCy_V0_Run%d_Vz%d",runNums[i],j+1),"",NbinVt,0,NbinVt,90,0,90);
617  fHist_znAx_V0_VxVy[i][j] = new TProfile2D(Form("fHist_znAx_V0_Run%d_Vz%d",runNums[i],j+1),"",NbinVt,0,NbinVt,90,0,90);
618  fHist_znAy_V0_VxVy[i][j] = new TProfile2D(Form("fHist_znAy_V0_Run%d_Vz%d",runNums[i],j+1),"",NbinVt,0,NbinVt,90,0,90);
619 
620  fListZDCQxy->Add(fHist_znCx_V0_VxVy[i][j]);
621  fListZDCQxy->Add(fHist_znCy_V0_VxVy[i][j]);
622  fListZDCQxy->Add(fHist_znAx_V0_VxVy[i][j]);
623  fListZDCQxy->Add(fHist_znAy_V0_VxVy[i][j]);
624  }
625  }
626  }
627  else{
628  for(int j=0;j<VzBinIter;j++) {
629  fHist_znCx_V0_VxVy[0][j] = new TProfile2D(Form("fHist_znCx_V0_Run%d_Vz%d",0,j+1),"",NbinVt,0,NbinVt,90,0,90);
630  fHist_znCy_V0_VxVy[0][j] = new TProfile2D(Form("fHist_znCy_V0_Run%d_Vz%d",0,j+1),"",NbinVt,0,NbinVt,90,0,90);
631  fHist_znAx_V0_VxVy[0][j] = new TProfile2D(Form("fHist_znAx_V0_Run%d_Vz%d",0,j+1),"",NbinVt,0,NbinVt,90,0,90);
632  fHist_znAy_V0_VxVy[0][j] = new TProfile2D(Form("fHist_znAy_V0_Run%d_Vz%d",0,j+1),"",NbinVt,0,NbinVt,90,0,90);
633 
634  fListZDCQxy->Add(fHist_znCx_V0_VxVy[0][j]);
635  fListZDCQxy->Add(fHist_znCy_V0_VxVy[0][j]);
636  fListZDCQxy->Add(fHist_znAx_V0_VxVy[0][j]);
637  fListZDCQxy->Add(fHist_znAy_V0_VxVy[0][j]);
638  }
639  }
640  }
641  }
642 
643 
644 
645 
646 
647 
648  Double_t centRange[12] = {0,5,10,20,30,40,50,60,70,80,90,100};
649  Double_t pTRange[21] = {0,0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3.0,3.3,3.6,4.0,4.5,5.0};
650  TString sNameComp[4] = {"uxCx","uyCy","uxAx","uyAy"};
651 
652 
653  if(bApplyRecent) {
654  //v2 integrated:
655  fHist_v2xV1_ZDN_Norm_All = new TProfile("fHist_v2xV1_ZDN_Norm_All","v2 X V1^2 (ZDC) all terms",11,centRange,"");
656  fHist_v2xV1_ZDN_Norm_All->Sumw2();
658  fHist_v2xV1_ZDN_Cent_All = new TProfile("fHist_v2xV1_ZDN_Cent_All","v2 X V1^2 (ZDC) all terms",100, 0, 100,"");
659  fHist_v2xV1_ZDN_Cent_All->Sumw2();
661  fHist_v2xV1_ZDN_Refm_All = new TProfile("fHist_v2xV1_ZDN_Refm_All","v2 X V1^2 (ZDC) all terms", 40, 0, 4000,"");
662  fHist_v2xV1_ZDN_Refm_All->Sumw2();
664 
665  fHist_ZDN_resol_Norm_All = new TProfile("fHist_ZDN_resol_Norm_All","Resol. V1^2(ZDC) All",11,centRange,"");
666  fHist_ZDN_resol_Norm_All->Sumw2();
668  fHist_ZDN_resol_Cent_All = new TProfile("fHist_ZDN_resol_Cent_All","Resol. V1^2(ZDC) All",100, 0, 100,"");
669  fHist_ZDN_resol_Cent_All->Sumw2();
671  fHist_ZDN_resol_Refm_All = new TProfile("fHist_ZDN_resol_Refm_All","Resol. V1^2(ZDC) All",40, 0, 4000,"");
672  fHist_ZDN_resol_Refm_All->Sumw2();
674 
675 
676  //v1 integrated, 4 terms, normal bin, fine Centr bin
677  for(int i=0; i<4; i++) {
678  sprintf(name,"fHist_v1_%s_ZDN_Norm",static_cast<const char*>(sNameComp[i]));
679  fHist_v2xV1_ZDN_Norm_Sep[i] = new TProfile(name,"v1 vs Cent",11,centRange,"");
680  fHist_v2xV1_ZDN_Norm_Sep[i]->Sumw2();
682 
683  sprintf(name,"fHist_v1_%s_ZDN_Cent",static_cast<const char*>(sNameComp[i]));
684  fHist_v2xV1_ZDN_Cent_Sep[i] = new TProfile(name,"v1 vs Cent",100,0,100,"");
685  fHist_v2xV1_ZDN_Cent_Sep[i]->Sumw2();
687  }
688  //v1 resulution, 2 terms, normal bin, fine Centr bin
689  fHist_ZDN_resol_Norm_Sep[0] = new TProfile("fHist_ZDN_resol_XaXc_Norm"," XaXc vs Cent",11,centRange,"");
690  fHist_ZDN_resol_Norm_Sep[0]->Sumw2();
692  fHist_ZDN_resol_Norm_Sep[1] = new TProfile("fHist_ZDN_resol_YaYc_Norm"," YaYc vs Cent",11,centRange,"");
693  fHist_ZDN_resol_Norm_Sep[1]->Sumw2();
695 
696  fHist_ZDN_resol_Cent_Sep[0] = new TProfile("fHist_ZDN_resol_XaXc_Cent"," XaXc vs Cent",100,0,100,"");
697  fHist_ZDN_resol_Cent_Sep[0]->Sumw2();
699  fHist_ZDN_resol_Cent_Sep[1] = new TProfile("fHist_ZDN_resol_YaYc_Cent"," YaYc vs Cent",100,0,100,"");
700  fHist_ZDN_resol_Cent_Sep[1]->Sumw2();
702 
703 
704 
705  //v2 differential pT:
706  for(int i=0; i<10; i++) {
707  sprintf(name,"fHist_v2xV1_ZDN_pTDiff_Cent%d",i);
708  fHist_v2xV1_ZDN_pTDiff_All[i] = new TProfile(name,"v2 X V1^2 pTdiff",20,pTRange);
709  fHist_v2xV1_ZDN_pTDiff_All[i]->Sumw2();
711  }
712 
713  //v1 differential Eta, pT:
714  for(int i=0; i<10; i++) {
715  for(int j=0; j<4; j++) {
716  sprintf(name,"fHist_%s_ZDN_EtaDiff_Cent%d",static_cast<const char*>(sNameComp[j]),i);
717  fHist_v1xV1_ZDN_EtaDiff[j][i] = new TProfile(name,"v1 Eta-diff", 5, -0.8, 0.8);
718  fHist_v1xV1_ZDN_EtaDiff[j][i]->Sumw2();
720 
721  sprintf(name,"fHist_%s_ZDN_pTDiff_Cent%d",static_cast<const char*>(sNameComp[j]),i);
722  fHist_v1xV1_ZDN_pTDiff[j][i] = new TProfile(name,"v1 Eta-diff", 20, pTRange);
723  fHist_v1xV1_ZDN_pTDiff[j][i]->Sumw2();
725  }
726  }
727 
728  //------------ calculate centrality weight: ------------
729  TH1F *fCent_fromDATA;
730 
731  if(fListZDCQxy) {
732  fCent_fromDATA = (TH1F *) fListZDCQxy->FindObject("fHist_Cent_afterr_ZDCcut");
733  }
734  else{
735  printf("\n\n ******** Running without Centrality weight !!******** \n\n");
736  fCent_fromDATA = new TH1F("fCent_DATA","Centrality distribution",100,0,100);
737  for(int i=1;i<=fCent_fromDATA->GetNbinsX();i++){
738  fCent_fromDATA->SetBinContent(i,100);
739  }
740  }
741 
742  Double_t Content,error;
743  Double_t maxCount = -1;
744  Int_t maxCent = -1;
745  Double_t weight = 0.;
746 
747  for(int i=1;i<=fCent_fromDATA->GetNbinsX();i++){
748  Content = fCent_fromDATA->GetBinContent(i);
749  if(maxCount < Content){
750  maxCount = Content;
751  maxCent = i;
752  }
753  }
754 
755  fWeight_Cent = new TH1F("fWeight_Cent","Weight for centrality",100,0,100);
756 
757  for(int i=1;i<=fWeight_Cent->GetNbinsX();i++)
758  {
759  Content = fCent_fromDATA->GetBinContent(i);
760  error = fCent_fromDATA->GetBinError(i);
761  if(Content){
762  weight = maxCount/Content;
763  }
764  else{
765  weight = 0;
766  }
767  fWeight_Cent->SetBinContent(i,weight);
768  fWeight_Cent->SetBinError(i,0.0);
769  //cout<<"cent "<<i-1<<"-"<<i<<" \t wgt = "<<fWeight_Cent->GetBinContent(i)<<" error = "<<fWeight_Cent->GetBinError(i)<<endl;
770  }
771 
773  //delete fCent_fromDATA;
774 
775 
776 
777  //---------Filter bit efficiency----------
778  if(fListHijing) {
779  for(int i=0;i<10;i++) {
780  fFB_Efficiency_Cent[i] = (TH1D *) fListHijing->FindObject(Form("eff_unbiased_%d",i));
781  //cout<<"cent = "<<i<<" name = "<<fFB_Efficiency_Cent[i]->GetName()<<endl;
782  }
783  }
784  else{ // if MC efficiency not found then use weight = 1. Code won't crash.
785  printf("\n\n !!***** Warning *****!! \n TList for FilterBit efficiency not found !!\n\n");
786  for(int i=0;i<10;i++){
787  fFB_Efficiency_Cent[i] = new TH1D(Form("eff_unbiased_%d",i),"",100,0,20);
788  for(int j=1;j<=fFB_Efficiency_Cent[i]->GetNbinsX();j++){
789  fFB_Efficiency_Cent[i]->SetBinContent(j,1.000);
790  }
791  }
792  }
793 
794  } //if(bApplyRecent)
795 
796 
797 
798 
799 
800 
801 
802 
803 
804  if(sCentEstimator=="V0"){
805  fHist_Task_config->Fill(1.5);
806  }
807  if(sCentEstimator=="TPC"){
808  fHist_Task_config->Fill(2.5);
809  }
810  if(sCentEstimator=="CL1"){
811  fHist_Task_config->Fill(3.5);
812  }
813  if(fRejectPileUp){
814  fHist_Task_config->Fill(4.5);
815  }
816  if(fRejectPileUpTight){
817  fHist_Task_config->Fill(5.5);
818  }
819  if(fAnalysisSet=="FillGainEq"){
820  fHist_Task_config->Fill(6.5);
821  }
822  if(fAnalysisSet=="DoGainEq"){
823  fHist_Task_config->Fill(7.5);
824  if(fListZDCWgt){
825  fHist_Task_config->Fill(8.5);
826  }
827  }
828  if(bFillCosSin){
829  fHist_Task_config->Fill(9.5);
830  }
831  if(bApplyRecent){
832  fHist_Task_config->Fill(10.5);
833  if(fListZDCQxy){
834  fHist_Task_config->Fill(11.5);
835  }
836  }
837  if(!bRunAveragedQn){
838  fHist_Task_config->Fill(12.5);
839  }
840 
841 
842 
843  fHist_CutParameters->SetBinContent(1,VzCut[0]);
844  fHist_CutParameters->SetBinContent(2,VzCut[1]);
845  fHist_CutParameters->SetBinContent(3,vzBin);
846  fHist_CutParameters->SetBinContent(4,VxCut[0]);
847  fHist_CutParameters->SetBinContent(5,VxCut[1]);
848  fHist_CutParameters->SetBinContent(6,vxBin);
849  fHist_CutParameters->SetBinContent(7,VyCut[0]);
850  fHist_CutParameters->SetBinContent(8,VyCut[1]);
851  fHist_CutParameters->SetBinContent(9,vyBin);
852 
853 
854 
855  PostData(1,fListHistos);
856 
857  if(!bApplyRecent && fAnalysisSet=="DoGainEq") {
858  PostData(2,fListZDCQxy);
859  }
860  else{
861  fListDummy1 = new TList();
862  fListDummy1->SetOwner(kTRUE);
864  PostData(2,fListDummy1);
865  }
866 
867 
868  fAnalysisUtil = new AliAnalysisUtils();
869  fAnalysisUtil->SetUseMVPlpSelection(kTRUE);
870  fAnalysisUtil->SetUseOutOfBunchPileUp(kTRUE);
871 
872 
873  fcheckOnce = 0;
874  fOldRunNum = 0;
875  printf("\n\n::UserCreateOutPutObject()\n Runflag= %d, dataset: %s, VxyBin = %d, Analysis= %s\n\n",frunflag,fDataSet.Data(),NbinVt,fAnalysisSet.Data());
876 
877 }
878 
879 //________________________________________________________________________
881 {
882 
883  float stepCount = 0.5;
884 
885  //printf("\n ... ::UserExec() is being called. 1 Step %d... \n",stepCount);
886 
887  fHist_Event_count->Fill(stepCount);
888  stepCount++;
889 
890  AliAODEvent *aod = dynamic_cast<AliAODEvent*>(InputEvent());
891  fEvent = dynamic_cast<AliFlowEventSimple*>(GetInputData(1));
892 
893  AliFlowEvent* anEvent = dynamic_cast<AliFlowEvent*>(GetInputData(2));
894 
895  AliFlowVector vQarray[2];
896 
897  if(anEvent) {
898  // Get Q vectors for the subevents
899  anEvent->GetZDC2Qsub(vQarray);
900  }
901  //I left here. I need to print ZDC q vectors from vQarray. See the class structure.
902  //how does it takes care of run by run events?
903 
904 
905 
906 
907  if(!aod){
908  printf("\n ... ::UserExec = no aod found..... \n");
909  return;
910  }
911 
912  fHist_Event_count->Fill(stepCount);
913  stepCount++;
914 
915 
916  AliAODVertex *pVertex = aod->GetPrimaryVertex();
917  Double_t Vxyz[3] = {0,0,0};
918  Vxyz[0] = pVertex->GetX();
919  Vxyz[1] = pVertex->GetY();
920  Vxyz[2] = pVertex->GetZ();
921 
922  //------- Apply Necessary event cuts ---------
923  if(Vxyz[2] >= VzCut[1] || Vxyz[2] <= VzCut[0]) return;
924 
925  fHist_Event_count->Fill(stepCount);
926  stepCount++;
927 
928  if(Vxyz[0] >= VxCut[1] || Vxyz[0] <= VxCut[0]) return;
929 
930  fHist_Event_count->Fill(stepCount);
931  stepCount++;
932 
933  if(Vxyz[1] >= VyCut[1] || Vxyz[1] <= VyCut[0]) return;
934 
935  fHist_Event_count->Fill(stepCount);
936  stepCount++;
937 
938 
939  Double_t EvtCent = fEvent->GetCentrality();
940 
941 
942  //---------- get runindex: --------------
943  Int_t runNumber = aod->GetRunNumber();
944  Int_t runindex = -111;
945 
946  for(int i=0;i<frunflag;i++){
947  if(runNumber==runNums[i])
948  {
949  runindex = i;
950  break;
951  }
952  }
953  if(runindex<0) {
954  printf("\n ... **WARNING** \n::UserExec() runnumber not listed.\n EXIT..\n");
955  //exit(1);
956  }
957  //-----------------------------------------
958 
959 
960 
961 
962 
963 
964 
965 
966 
967 
968 
969  //--------- starting pileup rejection work: --------
970  Double_t centrV0M=300;
971  Double_t centrCL1=300;
972  Double_t centrCL0=300;
973  Double_t centrTRK=300;
974 
975  if(fDataSet=="2010"||fDataSet=="2011"){
976  centrV0M = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("V0M");
977  centrCL1 = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("CL1");
978  centrCL0 = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("CL0");
979  centrTRK = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("TRK");
980  }
981  else{
982  fMultSelection = (AliMultSelection*) InputEvent()->FindListObject("MultSelection");
983  if(!fMultSelection) {
984  printf("\n\n **WARNING** ::UserExec() AliMultSelection object not found.\n\n");
985  exit(1);
986  }
987  centrV0M = fMultSelection->GetMultiplicityPercentile("V0M");
988  centrCL1 = fMultSelection->GetMultiplicityPercentile("CL1");
989  centrCL0 = fMultSelection->GetMultiplicityPercentile("CL0");
990  centrTRK = fMultSelection->GetMultiplicityPercentile("TRK");
991  }// 2015
992 
993  Bool_t BisPileup=kFALSE;
994 
995  if(fRejectPileUp && InputEvent()) {
996  //if(!fCutsEvent->IsSelected(InputEvent(),MCEvent())) return;
997  if(fDataSet!="2015") {
998  if(plpMV(aod)) {
999  fPileUpCount->Fill(0.5);
1000  BisPileup=kTRUE;
1001  }
1002  Int_t isPileup = aod->IsPileupFromSPD(3);
1003  if(isPileup != 0) {
1004  fPileUpCount->Fill(1.5);
1005  //BisPileup=kTRUE; // ----- Rihan: pileup from SPD not used for 2010
1006  }
1007  if(((AliAODHeader*)aod->GetHeader())->GetRefMultiplicityComb08() < 0) {
1008  fPileUpCount->Fill(2.5);
1009  BisPileup=kTRUE;
1010  }
1011  if(aod->IsIncompleteDAQ()) {
1012  fPileUpCount->Fill(3.5);
1013  BisPileup=kTRUE;
1014  }
1015 
1016  //check vertex consistency
1017  const AliAODVertex* vtTrc = aod->GetPrimaryVertex();
1018  const AliAODVertex* vtSPD = aod->GetPrimaryVertexSPD();
1019 
1020  if(vtTrc->GetNContributors() < 2 || vtSPD->GetNContributors()<1) {
1021  fPileUpCount->Fill(5.5);
1022  BisPileup=kTRUE;
1023  }
1024 
1025  double covTrc[6], covSPD[6];
1026  vtTrc->GetCovarianceMatrix(covTrc);
1027  vtSPD->GetCovarianceMatrix(covSPD);
1028 
1029  double dz = vtTrc->GetZ() - vtSPD->GetZ();
1030 
1031  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
1032  double errTrc = TMath::Sqrt(covTrc[5]);
1033  double nsigTot = dz/errTot;
1034  double nsigTrc = dz/errTrc;
1035 
1036  if(TMath::Abs(dz)>0.2 || TMath::Abs(nsigTot)>10 || TMath::Abs(nsigTrc)>20) {
1037  fPileUpCount->Fill(6.5);
1038  BisPileup=kTRUE;
1039  }
1040  if(fAnalysisUtil->IsPileUpEvent(InputEvent())) {
1041  fPileUpCount->Fill(7.5);
1042  BisPileup=kTRUE;
1043  }
1044  }
1045 
1046  else {
1047  // pileup from AliMultSelection for 2015
1048  if(!fMultSelection->GetThisEventIsNotPileup())
1049  fPileUpMultSelCount->Fill(0.5);
1050  if(!fMultSelection->GetThisEventIsNotPileupMV())
1051  fPileUpMultSelCount->Fill(1.5);
1052  if(!fMultSelection->GetThisEventIsNotPileupInMultBins())
1053  fPileUpMultSelCount->Fill(2.5);
1054  if(!fMultSelection->GetThisEventHasNoInconsistentVertices())
1055  fPileUpMultSelCount->Fill(3.5);
1056  if(!fMultSelection->GetThisEventPassesTrackletVsCluster())
1057  fPileUpMultSelCount->Fill(4.5);
1058  if(!fMultSelection->GetThisEventIsNotAsymmetricInVZERO())
1059  fPileUpMultSelCount->Fill(5.5);
1060  if(!fMultSelection->GetThisEventIsNotIncompleteDAQ())
1061  fPileUpMultSelCount->Fill(6.5);
1062  if(!fMultSelection->GetThisEventHasGoodVertex2016())
1063  fPileUpMultSelCount->Fill(7.5);
1064 
1065  BisPileup=kFALSE;
1066 
1067  // pile-up a la Dobrin for LHC15o
1068  if(plpMV(aod)) {
1069  fPileUpCount->Fill(0.5);
1070  BisPileup=kTRUE;
1071  }
1072 
1073  Int_t isPileup = aod->IsPileupFromSPD(3);
1074  if(isPileup != 0) {
1075  fPileUpCount->Fill(1.5);
1076  BisPileup=kTRUE;
1077  }
1078 
1079  if(((AliAODHeader*)aod->GetHeader())->GetRefMultiplicityComb08() < 0) {
1080  fPileUpCount->Fill(2.5);
1081  BisPileup=kTRUE;
1082  }
1083 
1084  if(aod->IsIncompleteDAQ()) {
1085  fPileUpCount->Fill(3.5);
1086  BisPileup=kTRUE;
1087  }
1088 
1089  if(fabs(centrV0M-centrCL1)>7.5) {
1090  fPileUpCount->Fill(4.5);
1091  BisPileup=kTRUE;
1092  }
1093 
1094  // check vertex consistency
1095  const AliAODVertex* vtTrc = aod->GetPrimaryVertex();
1096  const AliAODVertex* vtSPD = aod->GetPrimaryVertexSPD();
1097 
1098  if (vtTrc->GetNContributors() < 2 || vtSPD->GetNContributors()<1) {
1099  fPileUpCount->Fill(5.5);
1100  BisPileup=kTRUE;
1101  }
1102 
1103  double covTrc[6], covSPD[6];
1104  vtTrc->GetCovarianceMatrix(covTrc);
1105  vtSPD->GetCovarianceMatrix(covSPD);
1106 
1107  double dz = vtTrc->GetZ() - vtSPD->GetZ();
1108 
1109  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
1110  double errTrc = TMath::Sqrt(covTrc[5]);
1111  double nsigTot = dz/errTot;
1112  double nsigTrc = dz/errTrc;
1113 
1114  if(TMath::Abs(dz)>0.2 || TMath::Abs(nsigTot)>10 || TMath::Abs(nsigTrc)>20) {
1115  fPileUpCount->Fill(6.5);
1116  BisPileup=kTRUE;
1117  }
1118 
1119  // cuts on tracks
1120  const Int_t nTracks = aod->GetNumberOfTracks();
1121  Int_t multEsd = ((AliAODHeader*)aod->GetHeader())->GetNumberOfESDTracks();
1122 
1123  Int_t multTrk = 0;
1124  Int_t multTrkBefC = 0;
1125  Int_t multTrkTOFBefC = 0;
1126  Int_t multTPC = 0;
1127 
1128  for(Int_t it = 0; it < nTracks; it++) {
1129  AliAODTrack* aodTrk = (AliAODTrack*)aod->GetTrack(it);
1130  if(!aodTrk) {
1131  delete aodTrk;
1132  continue;
1133  }
1134 // if(aodTrk->TestFilterBit(32)){
1135 // multTrkBefC++;
1136 // if(TMath::Abs(aodTrk->GetTOFsignalDz()) <= 10. && aodTrk->GetTOFsignal() >= 12000. && aodTrk->GetTOFsignal() <= 25000.)
1137 // multTrkTOFBefC++;
1138 // if((TMath::Abs(aodTrk->Eta()) < 0.8) && (aodTrk->GetTPCNcls() >= 70) && (aodTrk->Pt() >= 0.2) && (aodTrk->Pt() < 20.))
1139 // multTrk++;
1140 // }
1141  if(aodTrk->TestFilterBit(128))
1142  multTPC++;
1143  } // end of for (Int_t it = 0; it < nTracks; it++)
1144 
1145  Double_t multTPCn = multTPC;
1146  Double_t multEsdn = multEsd;
1147  Double_t multESDTPCDif = multEsdn - multTPCn*3.38;
1148 
1149  if(multESDTPCDif > (fRejectPileUpTight?700.:15000.)) {
1150  fPileUpCount->Fill(7.5);
1151  BisPileup=kTRUE;
1152  }
1153 
1154  if(fRejectPileUpTight) {
1155  if(BisPileup==kFALSE) {
1156  if(!fMultSelection->GetThisEventIsNotPileup()) BisPileup=kTRUE;
1157  if(!fMultSelection->GetThisEventIsNotPileupMV()) BisPileup=kTRUE;
1158  if(!fMultSelection->GetThisEventIsNotPileupInMultBins()) BisPileup=kTRUE;
1159  if(!fMultSelection->GetThisEventHasNoInconsistentVertices()) BisPileup=kTRUE;
1160  if(!fMultSelection->GetThisEventPassesTrackletVsCluster()) BisPileup=kTRUE;
1161  if(!fMultSelection->GetThisEventIsNotIncompleteDAQ()) BisPileup=kTRUE;
1162  if(!fMultSelection->GetThisEventHasGoodVertex2016()) BisPileup=kTRUE;
1163  if(BisPileup) fPileUpCount->Fill(8.5);
1164  }
1165  }
1166  }
1167  }
1168 
1169  //----------- pile up rejection done ---------
1170  if(BisPileup) {
1171  return;
1172  }
1173 
1174  fHist_Event_count->Fill(stepCount);
1175  stepCount++;
1176 
1177 
1178 
1179 
1180 
1181 
1182 
1183  if(runNumber!=fOldRunNum){ //if runNumber changed, re-read new list.
1184  fcheckOnce = 0;
1185  }
1186 
1187 
1188  if(!fcheckOnce && fAnalysisSet=="DoGainEq") {
1189  fHist_ChanWgt_ZDCC = (TH1F *) fListZDCWgt->FindObject(Form("fHist1F_ZDCC_ChannelWgt_Run%d",runNumber));
1190  fHist_ChanWgt_ZDCA = (TH1F *) fListZDCWgt->FindObject(Form("fHist1F_ZDCA_ChannelWgt_Run%d",runNumber));
1191 
1192  fcheckOnce++;
1193  fOldRunNum = runNumber;
1194  }
1195  /*
1196  if(!fcheckOnce && fAnalysisSet=="DoGainEq") {
1197  fHist_ChanWgt_ZDCC = (TH1F *) fListZDCWgt->FindObject(Form("fHist1F_ZDCC_ChannelWgt_Run%d",runNumber));
1198  fHist_ChanWgt_ZDCA = (TH1F *) fListZDCWgt->FindObject(Form("fHist1F_ZDCA_ChannelWgt_Run%d",runNumber));
1199 
1200  if(bApplyRecent) {
1201  if(!bRunAveragedQn) {
1202  if(fListZDCQxy) {
1203  for(int i=0; i<vzBin; i++) {
1204  fHist_Recenter_ZDCCx[i] = (TH2F *) fListZDCQxy->FindObject(Form("fHist2F_znCx_V0_Run%d_Vz%d",runNumber,i+1));
1205  fHist_Recenter_ZDCCy[i] = (TH2F *) fListZDCQxy->FindObject(Form("fHist2F_znCy_V0_Run%d_Vz%d",runNumber,i+1));
1206  fHist_Recenter_ZDCAx[i] = (TH2F *) fListZDCQxy->FindObject(Form("fHist2F_znAx_V0_Run%d_Vz%d",runNumber,i+1));
1207  fHist_Recenter_ZDCAy[i] = (TH2F *) fListZDCQxy->FindObject(Form("fHist2F_znAy_V0_Run%d_Vz%d",runNumber,i+1));
1208  }
1209  }
1210  }
1211  else {
1212  if(fListZDCQxy) {
1213  for(int i=0; i<vzBin; i++) {
1214  fHist_Recenter_ZDCCx[i] = (TH2F *) fListZDCQxy->FindObject(Form("fHist2F_znCx_V0_Run%d_Vz%d",0,i+1));
1215  fHist_Recenter_ZDCCy[i] = (TH2F *) fListZDCQxy->FindObject(Form("fHist2F_znCy_V0_Run%d_Vz%d",0,i+1));
1216  fHist_Recenter_ZDCAx[i] = (TH2F *) fListZDCQxy->FindObject(Form("fHist2F_znAx_V0_Run%d_Vz%d",0,i+1));
1217  fHist_Recenter_ZDCAy[i] = (TH2F *) fListZDCQxy->FindObject(Form("fHist2F_znAy_V0_Run%d_Vz%d",0,i+1));
1218  }
1219  }
1220  }
1221  }
1222 
1223  fcheckOnce++;
1224  fOldRunNum = runNumber;
1225  } */
1226 
1227 
1228 
1229 
1230  Double_t ChanWgtZDCC[5] = {1.,1.,1.,1.,1.};
1231  Double_t ChanWgtZDCA[5] = {1.,1.,1.,1.,1.};
1232 
1233  Int_t iCentBin = abs(EvtCent) + 1;
1234  Int_t iWgtBin = -1;
1235  Int_t iCommon = -1;
1236 
1237 
1238  if(fAnalysisSet=="DoGainEq") {
1240 
1241  for(int ich=1; ich<=5; ich++){
1242  iWgtBin = 5*(iCentBin-1) + ich;
1243  ChanWgtZDCC[ich-1] = fHist_ChanWgt_ZDCC->GetBinContent(iWgtBin);
1244  ChanWgtZDCA[ich-1] = fHist_ChanWgt_ZDCA->GetBinContent(iWgtBin);
1245  }
1246  }
1247  else{
1248  //printf("\n\n **WARNING**\n ZDC Channel Weights not found. Using weights = 1.0 \n\n");
1249  //exit(1);
1250  }
1251  }
1252 
1253 
1254 
1255  fHist_Cent_woZDCcut->Fill(EvtCent);
1256 
1257 
1258  //----------- Read ZDC information ----------
1259  AliAODZDC *aodZDC = aod->GetZDCData();
1260  Float_t fZDCGainAlpha = 0.500; // fZDCGainAlpha : Jacopo uses 0.35 ??
1261  Float_t energyZNC = (Float_t) (aodZDC->GetZNCEnergy());
1262  Float_t energyZPC = (Float_t) (aodZDC->GetZPCEnergy());
1263  Float_t energyZNA = (Float_t) (aodZDC->GetZNAEnergy());
1264  Float_t energyZPA = (Float_t) (aodZDC->GetZPAEnergy());
1265  Float_t energyZEM1 = (Float_t) (aodZDC->GetZEM1Energy());
1266  Float_t energyZEM2 = (Float_t) (aodZDC->GetZEM2Energy());
1267 
1268  const Double_t * towZNC = aodZDC->GetZNCTowerEnergy();
1269  const Double_t * towZPC = aodZDC->GetZPCTowerEnergy();
1270  const Double_t * towZNA = aodZDC->GetZNATowerEnergy();
1271  const Double_t * towZPA = aodZDC->GetZPATowerEnergy();
1272 
1273  const Double_t * towZNClg = aodZDC->GetZNCTowerEnergyLR(); // Low gain something, should not be used.
1274  const Double_t * towZNAlg = aodZDC->GetZNATowerEnergyLR();
1275 
1276  Double_t towZPClg[5] = {0.,};
1277  Double_t towZPAlg[5] = {0.,};
1278 
1279  for(Int_t it=0; it<5; it++) {
1280  towZPClg[it] = 8*towZPC[it];
1281  towZPAlg[it] = 8*towZNA[it];
1282  }
1283 
1284  Int_t BadChannel = 0;
1285 
1286  //sanity: remove if any of ZDC_C_A has negetive Energy:
1287  //This was a bad choice of cut. I need to see the Physics effect later..!!
1288  //if(towZNC[1]<0 || towZNC[2]<0 || towZNC[3]<0 || towZNC[4]<0) return;
1289 
1290  for(int i=0; i<4; i++) {
1291  if(towZNC[i] < 0) {
1292  BadChannel++;
1293  }
1294  }
1295 
1296  if(BadChannel>=2) return; // Remove Event if more than one channel has Energy < 0
1297 
1298  fHist_Event_count->Fill(stepCount);
1299  stepCount++;
1300 
1301  //if(towZNA[1]<0 || towZNA[2]<0 || towZNA[3]<0 || towZNA[4]<0) return;
1302  BadChannel = 0;
1303 
1304  for(int i=0; i<4; i++) {
1305  if(towZNA[i] < 0) {
1306  BadChannel++;
1307  }
1308  }
1309 
1310  if(BadChannel>=2) return;
1311 
1312  fHist_Event_count->Fill(stepCount);
1313  stepCount++;
1314 
1315 
1316 
1317 
1318 
1319 
1320 
1321 
1322 
1323 
1324 
1325 
1326 //********** Get centroid from ZDCs **************
1327 
1328  Double_t xyZNC[2]={0.,0.};
1329  Double_t xyZNA[2]={0.,0.};
1330 
1331  Float_t zncEnergy=0., znaEnergy=0.;
1332 
1333 
1334 
1335 /*-----------------------Not used---------------------------------
1336  Int_t CenBin = GetCenBin(centrperc);
1337  Double_t zvtxpos[3]={0.,0.,0.};
1338  fFlowEvent->GetVertexPosition(zvtxpos);
1339  Int_t RunNum=fFlowEvent->GetRun();
1340  if(fTowerEqList) {
1341  if(RunNum!=fCachedRunNum) {
1342  for(Int_t i=0; i<8; i++) {
1343  fTowerGainEq[i] = (TH1D*)(fTowerEqList->FindObject(Form("Run %d",RunNum))->FindObject(Form("fhnTowerGainEqFactor[%d][%d]",RunNum,i)));
1344  }
1345  }
1346  }
1347  Bool_t fUseMCCen = kFALSE; //rihan:hardcoded
1348  if (fUseMCCen) {
1349  if(aod->GetRunNumber() < 209122)
1350  aodZDC->GetZNCentroidInPbPb(1380., xyZNC, xyZNA);
1351  else
1352  aodZDC->GetZNCentroidInPbPb(2510., xyZNC, xyZNA);
1353  }
1354  else {
1355  //set tower gain equalization, if available
1356  if(fTowerEqList) {
1357  for(Int_t i=0; i<8; i++)
1358  {
1359  if(fTowerGainEq[i])
1360  AvTowerGain[i] = fTowerGainEq[i]->GetBinContent(fTowerGainEq[i]->FindBin(centrperc));
1361  }
1362  }//--------------------------------------------------------------- */
1363 
1364 
1365 
1366  Double_t towCalibZNC[5] = {0,};
1367  Double_t towCalibZNA[5] = {0,};
1368 
1369  // towZNC[] is constant; so Need to make a copy
1370  for(int i=0;i<5;i++){
1371  towCalibZNC[i] = towZNC[i];
1372  towCalibZNA[i] = towZNA[i];
1373  }
1374 
1375  // If Ei < 0 , Ei = 0 (Maxim's idea)
1376  for(int i=1;i<5;i++){
1377  if(towCalibZNC[i] < 0) towCalibZNC[i] = 0;
1378  if(towCalibZNA[i] < 0) towCalibZNA[i] = 0;
1379  }
1380 
1381 
1382 
1383  Int_t indexVx = fHist_Vx_ArrayFinder->FindBin(Vxyz[0]);
1384  Int_t indexVy = fHist_Vy_ArrayFinder->FindBin(Vxyz[1]);
1385  Int_t indexVz = fHist_Vz_ArrayFinder->FindBin(Vxyz[2]);
1386 
1387 
1388  Double_t tVertexBin1 = 0;
1389  tVertexBin1 = (Double_t) (indexVy-1)*vxBin + (Double_t)indexVx - 0.5 ;
1390 
1391 
1392  if(fAnalysisSet=="FillGainEq") {
1393  fHist_Vx_vs_runnum->Fill(runindex,Vxyz[0]);
1394  fHist_Vy_vs_runnum->Fill(runindex,Vxyz[1]);
1395  fHist_Vz_vs_runnum->Fill(runindex,Vxyz[2]);
1396 
1397  fHist_ZDCC_En_CommonCh[indexVz-1]->Fill(EvtCent,tVertexBin1,towCalibZNC[0]);
1398  fHist_ZDCA_En_CommonCh[indexVz-1]->Fill(EvtCent,tVertexBin1,towCalibZNA[0]);
1399 
1400  for(int ich=0; ich<5; ich++) {
1401  fHist_ZDCC_En_Run[runindex]->Fill(EvtCent,ich,towCalibZNC[ich]);
1402  fHist_ZDCA_En_Run[runindex]->Fill(EvtCent,ich,towCalibZNA[ich]);
1403  }
1404  }
1405 
1406 
1407 
1408 
1409 
1410 
1411 
1412  // Now calibrate the energy of channel [0-4]:
1413  for(int i=0;i<5;i++){
1414  towCalibZNC[i] = ChanWgtZDCC[i]*towCalibZNC[i];
1415 
1416  if(ChanWgtZDCA[i] < 4.){
1417  towCalibZNA[i] = ChanWgtZDCA[i]*towCalibZNA[i];
1418  }
1419  }
1420 
1421  //manually get Energy in ZDC-A channel [2]:
1422  if(ChanWgtZDCA[2] >= 4.0){
1423  towCalibZNA[2] = towCalibZNA[0] - towCalibZNA[1] - towCalibZNA[3] - towCalibZNA[4];
1424  }
1425 
1426 
1427  for(Int_t i=0; i<5; i++){
1428  zncEnergy += towCalibZNC[i];
1429  znaEnergy += towCalibZNA[i];
1430  }
1431 
1432 
1433  Double_t AvTowerGain[8] = {1., 1., 1., 1., 1., 1., 1., 1.};
1434 
1435  const Float_t x[4] = {-1.75, 1.75,-1.75, 1.75};
1436  const Float_t y[4] = {-1.75, -1.75, 1.75, 1.75};
1437 
1438  Float_t numXZNC=0., numYZNC=0., denZNC=0., wZNC;
1439  Float_t numXZNA=0., numYZNA=0., denZNA=0., wZNA;
1440 
1441  for(Int_t i=0; i<4; i++) {
1442  if(towCalibZNC[i+1]>0.) {
1443  wZNC = TMath::Power(towCalibZNC[i+1], fZDCGainAlpha)*AvTowerGain[i];
1444  numXZNC += x[i]*wZNC;
1445  numYZNC += y[i]*wZNC;
1446  denZNC += wZNC;
1447  }
1448 
1449  if(towCalibZNA[i+1]>0.) {
1450  wZNA = TMath::Power(towCalibZNA[i+1], fZDCGainAlpha)*AvTowerGain[i+4];
1451  numXZNA += x[i]*wZNA;
1452  numYZNA += y[i]*wZNA;
1453  denZNA += wZNA;
1454  }
1455  }
1456 
1457  if(denZNC!=0) {
1458  xyZNC[0] = numXZNC/denZNC;
1459  xyZNC[1] = numYZNC/denZNC;
1460  }
1461  else{
1462  xyZNC[0] = 999.;
1463  xyZNC[1] = 999.;
1464  zncEnergy = 0.;
1465  }
1466  if(denZNA!=0) {
1467  xyZNA[0] = numXZNA/denZNA;
1468  xyZNA[1] = numYZNA/denZNA;
1469  }
1470  else{
1471  xyZNA[0] = 999.;
1472  xyZNA[1] = 999.;
1473  znaEnergy = 0.;
1474  }
1475 
1476 
1477 
1478  xyZNA[0] = -1.*xyZNA[0]; //----- Important: zdcA_X = -zdcA_X ---------
1479 
1480 
1481 
1482  if(sqrt(xyZNC[0]*xyZNC[0] + xyZNC[1]*xyZNC[1])>1.5) return;
1483 
1484  fHist_Event_count->Fill(stepCount);
1485  stepCount++;
1486 
1487  if(sqrt(xyZNA[0]*xyZNA[0] + xyZNA[1]*xyZNA[1])>1.5) return;
1488 
1489  fHist_Event_count->Fill(stepCount);
1490  stepCount++;
1491 
1492 
1493 
1494  Double_t psi1C = TMath::ATan2(xyZNC[1],xyZNC[0]);
1495  if(psi1C<0){
1496  psi1C += 2.*TMath::Pi();
1497  }
1498  Double_t psi1A = TMath::ATan2(xyZNA[1],xyZNA[0]);
1499  if(psi1A<0){
1500  psi1A += 2.*TMath::Pi();
1501  }
1502 
1503  fHist_Psi1_ZDCC_wGainCorr->Fill(psi1C);
1504  fHist_Psi1_ZDCA_wGainCorr->Fill(psi1A);
1505 
1506 
1507  if(fAnalysisSet=="DoGainEq") {
1508 
1509  tVertexBin1 = (Double_t) (indexVy-1)*vxBin + (Double_t)indexVx - 0.5 ;
1510 
1511  if(!bApplyRecent) {
1512  if(bFillCosSin) {
1513  if(!bRunAveragedQn){
1514  fHist_znCx_V0_VxVy[runindex][indexVz-1]->Fill(tVertexBin1,EvtCent,TMath::Cos(psi1C));
1515  fHist_znCy_V0_VxVy[runindex][indexVz-1]->Fill(tVertexBin1,EvtCent,TMath::Sin(psi1C));
1516  fHist_znAx_V0_VxVy[runindex][indexVz-1]->Fill(tVertexBin1,EvtCent,TMath::Cos(psi1A));
1517  fHist_znAy_V0_VxVy[runindex][indexVz-1]->Fill(tVertexBin1,EvtCent,TMath::Sin(psi1A));
1518  }
1519  else{
1520  fHist_znCx_V0_VxVy[0][indexVz-1]->Fill(tVertexBin1,EvtCent,TMath::Cos(psi1C));
1521  fHist_znCy_V0_VxVy[0][indexVz-1]->Fill(tVertexBin1,EvtCent,TMath::Sin(psi1C));
1522  fHist_znAx_V0_VxVy[0][indexVz-1]->Fill(tVertexBin1,EvtCent,TMath::Cos(psi1A));
1523  fHist_znAy_V0_VxVy[0][indexVz-1]->Fill(tVertexBin1,EvtCent,TMath::Sin(psi1A));
1524  }
1525  }
1526  else{
1527  if(!bRunAveragedQn){
1528  fHist_znCx_V0_VxVy[runindex][indexVz-1]->Fill(tVertexBin1,EvtCent,xyZNC[0]);
1529  fHist_znCy_V0_VxVy[runindex][indexVz-1]->Fill(tVertexBin1,EvtCent,xyZNC[1]);
1530  fHist_znAx_V0_VxVy[runindex][indexVz-1]->Fill(tVertexBin1,EvtCent,xyZNA[0]);
1531  fHist_znAy_V0_VxVy[runindex][indexVz-1]->Fill(tVertexBin1,EvtCent,xyZNA[1]);
1532  }
1533  else{
1534  fHist_znCx_V0_VxVy[0][indexVz-1]->Fill(tVertexBin1,EvtCent,xyZNC[0]);
1535  fHist_znCy_V0_VxVy[0][indexVz-1]->Fill(tVertexBin1,EvtCent,xyZNC[1]);
1536  fHist_znAx_V0_VxVy[0][indexVz-1]->Fill(tVertexBin1,EvtCent,xyZNA[0]);
1537  fHist_znAy_V0_VxVy[0][indexVz-1]->Fill(tVertexBin1,EvtCent,xyZNA[1]);
1538  }
1539  }
1540  }
1541  }
1542 
1543 
1545 
1546  if(fAnalysisSet!="FillGainEq" && bFillCosSin) {
1547  xyZNC[0] = TMath::Cos(psi1C);
1548  xyZNC[1] = TMath::Sin(psi1C);
1549  xyZNA[0] = TMath::Cos(psi1A);
1550  xyZNA[1] = TMath::Sin(psi1A);
1551  }
1552 
1553 
1554 
1555  if(bFillZDCQAon){
1556 
1557  //Double_t FillVsWith[5] = {EvtCent,static_cast<Double_t>(nRefMult), Vxyz[0], Vxyz[1], Vxyz[2]};
1558  Double_t FillVsWith[5] = {EvtCent, fRefMult, Vxyz[0], Vxyz[1], Vxyz[2]};
1559  Double_t FillValueQx[4] = {xyZNA[0],xyZNC[0],xyZNA[1],xyZNC[1]};
1560  Double_t FillValueXX[4] = {xyZNA[0]*xyZNC[0],xyZNA[1]*xyZNC[1],xyZNC[0]*xyZNA[1],xyZNC[1]*xyZNA[0]}; //XaXc,YaYc,XcYa,YcXa
1561 
1562  for(int i=0;i<4;i++){
1563  for(int j=0;j<5;j++){
1564  fHist_Qx_vs_Obs_woCorr[i][j]->Fill(FillVsWith[j],FillValueQx[i]);
1565  fHist_XX_vs_Obs_woCorr[i][j]->Fill(FillVsWith[j],FillValueXX[i]);
1566  }
1567  }
1568  }
1569 
1570 
1571  if(fAnalysisSet=="FillGainEq" && bFillCosSin) {
1572  xyZNC[0] = TMath::Cos(psi1C);
1573  xyZNC[1] = TMath::Sin(psi1C);
1574  xyZNA[0] = TMath::Cos(psi1A);
1575  xyZNA[1] = TMath::Sin(psi1A);
1576  }
1577 
1578 
1579  Double_t meanCx = 0.;
1580  Double_t meanCy = 0.;
1581  Double_t meanAx = 0.;
1582  Double_t meanAy = 0.;
1583 
1584 
1585  //Apply the recentering:
1586  /*
1587  if(bApplyRecent) {
1588  Int_t tVertexBin2 = (indexVy-1)*vxBin + indexVx;
1589  if(fListZDCQxy) {
1590  meanCx = fHist_Recenter_ZDCCx[indexVz-1]->GetBinContent(tVertexBin2,iCentBin);
1591  meanCy = fHist_Recenter_ZDCCy[indexVz-1]->GetBinContent(tVertexBin2,iCentBin);
1592  meanAx = fHist_Recenter_ZDCAx[indexVz-1]->GetBinContent(tVertexBin2,iCentBin);
1593  meanAy = fHist_Recenter_ZDCAy[indexVz-1]->GetBinContent(tVertexBin2,iCentBin);
1594  }
1595  }
1596 
1597  xyZNC[0] = xyZNC[0] - meanCx;
1598  xyZNC[1] = xyZNC[1] - meanCy;
1599  xyZNA[0] = xyZNA[0] - meanAx;
1600  xyZNA[1] = xyZNA[1] - meanAy;
1601  */
1602 
1603  //use Jacopo's flowEvent ZDC Q-vectors:
1604  xyZNC[0] = vQarray[0].Px();
1605  xyZNC[1] = vQarray[0].Py();
1606  xyZNA[0] = vQarray[1].Px();
1607  xyZNA[1] = vQarray[1].Py();
1608 
1609 
1610  double Psi1C = TMath::ATan2(xyZNC[1],xyZNC[0]);
1611  if(Psi1C<0){
1612  Psi1C += 2.*TMath::Pi();
1613  }
1614  double Psi1A = TMath::ATan2(xyZNA[1],xyZNA[0]);
1615  if(Psi1A<0){
1616  Psi1A += 2.*TMath::Pi();
1617  }
1618 
1619  if(Psi1C==0 && Psi1A==0) return;
1620 
1621  fHist_Event_count->Fill(stepCount);
1622  stepCount++;
1623 
1624 
1625  if(bFillZDCQAon){
1626  //Double_t FillVsWithNew[5] = {EvtCent,static_cast<Double_t>(nRefMult), Vxyz[0], Vxyz[1], Vxyz[2]};
1627  Double_t FillVsWithNew[5] = {EvtCent, fRefMult, Vxyz[0], Vxyz[1], Vxyz[2]};
1628  Double_t FillValueQxNew[4] = {xyZNA[0],xyZNC[0],xyZNA[1],xyZNC[1]};
1629  Double_t FillValueXXNew[4] = {xyZNA[0]*xyZNC[0],xyZNA[1]*xyZNC[1],xyZNC[0]*xyZNA[1],xyZNC[1]*xyZNA[0]}; //XaXc,YaYc,XcYa,YcXa
1630 
1631  for(int i=0;i<4;i++){
1632  for(int j=0;j<5;j++){
1633  fHist_Qx_vs_Obs_wiCorr[i][j]->Fill(FillVsWithNew[j],FillValueQxNew[i]);
1634  fHist_XX_vs_Obs_wiCorr[i][j]->Fill(FillVsWithNew[j],FillValueXXNew[i]);
1635  }
1636  }
1637  }
1638 
1639 
1640 
1641 
1642  Int_t iTracks = fEvent->NumberOfTracks();
1643  AliFlowTrackSimple* pTrack = NULL;
1644  Double_t Qnx_TPC[3] = {0,};
1645  Double_t Qny_TPC[3] = {0,};
1646  Double_t psi2,dPhi,dPt,dEta;
1647  Double_t pTwgt = 1.0;
1648  Double_t npoiMult = 0;
1649  Double_t dUx = 0;
1650  Double_t dUy = 0;
1651  Int_t ipTBin = 1;
1652  Int_t cIndex = -1;
1653  Double_t nRefMult = 0;
1654 
1655  if(EvtCent<5.0) { cIndex = 0; }
1656  else if(EvtCent>=5.0 && EvtCent<10){
1657  cIndex = 1;
1658  }
1659  else{
1660  cIndex = abs(EvtCent/10.0) + 1;
1661  }
1662 
1663  Double_t fullTerm = 0;
1664  Double_t fullReso = 0;
1665  Double_t CentWgt = 1.0;
1666 
1667  if(bApplyRecent) {
1668  for(int i=0; i<iTracks; i++) {
1669  pTrack = fEvent->GetTrack(i);
1670  if(!pTrack) continue;
1671  dPhi = pTrack->Phi();
1672  dPt = pTrack-> Pt();
1673  dEta = pTrack->Eta();
1674  //dChrg = pTrack->Charge();
1675  if(fabs(dEta) > 0.8) continue;
1676  if(dPt<0.20 || dPt>50.0) continue;
1677  if(!pTrack->IsPOItype(1)) continue;
1678  nRefMult++;
1679 
1680  ipTBin = fFB_Efficiency_Cent[cIndex]->FindBin(dPt);
1681  pTwgt = 1.0/fFB_Efficiency_Cent[cIndex]->GetBinContent(ipTBin);
1682 
1683 
1684  dUx = cos(dPhi);
1685  dUy = sin(dPhi);
1686 
1687  fHist_v1xV1_ZDN_EtaDiff[0][cIndex]->Fill(dEta, dUx*xyZNC[0], pTwgt); //uxCx
1688  fHist_v1xV1_ZDN_EtaDiff[1][cIndex]->Fill(dEta, dUy*xyZNC[1], pTwgt); //uyCy
1689  fHist_v1xV1_ZDN_EtaDiff[2][cIndex]->Fill(dEta, dUx*xyZNA[0], pTwgt); //uxAx
1690  fHist_v1xV1_ZDN_EtaDiff[3][cIndex]->Fill(dEta, dUy*xyZNA[1], pTwgt); //uyAy
1691 
1692  fHist_v1xV1_ZDN_pTDiff[0][cIndex]->Fill(dPt, dUx*xyZNC[0], pTwgt); //uxCx
1693  fHist_v1xV1_ZDN_pTDiff[1][cIndex]->Fill(dPt, dUy*xyZNC[1], pTwgt); //uyCy
1694  fHist_v1xV1_ZDN_pTDiff[2][cIndex]->Fill(dPt, dUx*xyZNA[0], pTwgt); //uxAx
1695  fHist_v1xV1_ZDN_pTDiff[3][cIndex]->Fill(dPt, dUy*xyZNA[1], pTwgt); //uyAy
1696 
1697  Qnx_TPC[0] += dUx*pTwgt;
1698  Qny_TPC[0] += dUy*pTwgt;
1699 
1700 
1701  dUx = cos(2.*dPhi);
1702  dUy = sin(2.*dPhi);
1703 
1704  fullTerm = dUx*xyZNA[0]*xyZNC[0]-dUx*xyZNA[1]*xyZNC[1]+dUy*xyZNA[0]*xyZNC[1]+dUy*xyZNA[1]*xyZNC[0];
1705  fHist_v2xV1_ZDN_pTDiff_All[cIndex]->Fill(dPt,fullTerm,pTwgt);
1706 
1707  Qnx_TPC[1] += dUx*pTwgt;
1708  Qny_TPC[1] += dUy*pTwgt;
1709 
1710  //cout<<"i = "<<i<<" iCent = "<<cIndex<<" pTBin "<<ipTBin<<" pt = "<<dPt<<"\t wgt = "<<pTwgt<<endl;
1711  npoiMult += pTwgt;
1712  }
1713  }
1714 
1715 
1716 
1717  if(npoiMult>0){
1718  for(int i=0;i<2;i++){
1719  Qnx_TPC[i] = Qnx_TPC[i]/npoiMult;
1720  Qny_TPC[i] = Qny_TPC[i]/npoiMult;
1721  }
1722  }
1723  else{
1724  for(int i=0;i<2;i++){
1725  Qnx_TPC[i] = 0;
1726  Qny_TPC[i] = 0;
1727  }
1728  }
1729 
1730  fullTerm = 0;
1731  fullReso = 0;
1732 
1733  if(bApplyRecent) {
1734 
1735  fHist_Psi1_ZDCC_wRectCorr->Fill(Psi1C);
1736  fHist_Psi1_ZDCA_wRectCorr->Fill(Psi1A);
1737 
1738  CentWgt = fWeight_Cent->GetBinContent(iCentBin);
1739 
1740  fHist_ZDN_resol_Norm_Sep[0]->Fill(EvtCent, xyZNA[0]*xyZNC[0]);
1741  fHist_ZDN_resol_Cent_Sep[0]->Fill(EvtCent, xyZNA[0]*xyZNC[0]);
1742 
1743  fHist_ZDN_resol_Norm_Sep[1]->Fill(EvtCent, xyZNA[1]*xyZNC[1]);
1744  fHist_ZDN_resol_Cent_Sep[1]->Fill(EvtCent, xyZNA[1]*xyZNC[1]);
1745 
1746  Double_t v1FillTerms[4] = {Qnx_TPC[0]*xyZNC[0], Qny_TPC[0]*xyZNC[1], Qnx_TPC[0]*xyZNA[0], Qny_TPC[0]*xyZNA[1]};
1747 
1748 
1749  for(int i=0; i<4; i++){
1750  fHist_v2xV1_ZDN_Norm_Sep[i]->Fill(EvtCent, v1FillTerms[i]);
1751  fHist_v2xV1_ZDN_Cent_Sep[i]->Fill(EvtCent, v1FillTerms[i]);
1752  }
1753 
1754 
1755  fullTerm = Qnx_TPC[1]*(xyZNA[0]*xyZNC[0] - xyZNA[1]*xyZNC[1]) + Qny_TPC[1]*(xyZNA[0]*xyZNC[1] + xyZNA[1]*xyZNC[0]);
1756  fullReso = xyZNA[0]*xyZNC[0] + xyZNA[1]*xyZNC[1];
1757 
1758  fHist_v2xV1_ZDN_Norm_All->Fill( EvtCent, fullTerm, CentWgt);
1759  fHist_v2xV1_ZDN_Cent_All->Fill( EvtCent, fullTerm, CentWgt);
1760  fHist_v2xV1_ZDN_Refm_All->Fill(nRefMult, fullTerm, CentWgt);
1761 
1762  fHist_ZDN_resol_Norm_All->Fill( EvtCent, fullReso, CentWgt);
1763  fHist_ZDN_resol_Cent_All->Fill( EvtCent, fullReso, CentWgt);
1764  fHist_ZDN_resol_Refm_All->Fill(nRefMult, fullReso, CentWgt);
1765  }
1766 
1767 
1768 
1769  fHist_Event_count->Fill(stepCount);
1770  stepCount++;
1771 
1772 
1773  fHist_Cent_wiZDCcut->Fill(EvtCent);
1774 
1775 
1776 
1777  //if(fievent%10==0) {
1778  //std::cout<<fievent<<" cTPC= "<<EvtCent<<"\t wZDA1 = "<<ChanWgtZDCA[1]<<"\t wZDA2 = "<<ChanWgtZDCA[2]<<"\tRefMult = "<<nRefMult<<std::endl;
1779  //std::cout<<" Cx = "<<vQarray[0].Px()<<"\t Cy = "<<vQarray[0].Py()<<"\t Ax = "<<vQarray[1].Px()<<"\t Ay = "<<vQarray[1].Py()<<std::endl;
1780  //}
1781 
1782 
1783  PostData(1,fListHistos);
1784 
1785  if(!bApplyRecent && fAnalysisSet=="DoGainEq") {
1786  PostData(2,fListZDCQxy);
1787  }
1788  else{
1789  PostData(2,fListDummy1); //default if no 'fAnalysisSet' found.
1790  }
1791 
1792 
1793  fievent++;
1794 
1795  //printf("\n ... ::UserExec() is being called. Step last %d... Event %d \n",stepCount,fievent);
1796 
1797 }// UserExec ends
1798 
1799 
1800 
1801 
1803 {
1804  AliDebug(2,"\n ... AliAnalysisTaskZDCGainEq::Terminate() is being called ... \n");
1805 }
1806 
1807 
1808 double AliAnalysisTaskZDCGainEq::GetWDist(const AliVVertex* v0, const AliVVertex* v1)
1809 {
1810  // calculate sqrt of weighted distance to other vertex
1811  if (!v0 || !v1) {
1812  AliDebug(2,"\n\n ::GetWDist => One of vertices is not valid\n\n");
1813  return 0;
1814  }
1815  static TMatrixDSym vVb(3);
1816  double dist = -1;
1817  double dx = v0->GetX()-v1->GetX();
1818  double dy = v0->GetY()-v1->GetY();
1819  double dz = v0->GetZ()-v1->GetZ();
1820  double cov0[6],cov1[6];
1821  v0->GetCovarianceMatrix(cov0);
1822  v1->GetCovarianceMatrix(cov1);
1823  vVb(0,0) = cov0[0]+cov1[0];
1824  vVb(1,1) = cov0[2]+cov1[2];
1825  vVb(2,2) = cov0[5]+cov1[5];
1826  vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
1827  vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
1828  vVb.InvertFast();
1829  if (!vVb.IsValid()) {
1830  AliDebug(2,"Singular Matrix\n");
1831  return dist;
1832  }
1833  dist = vVb(0,0)*dx*dx + vVb(1,1)*dy*dy + vVb(2,2)*dz*dz
1834  + 2*vVb(0,1)*dx*dy + 2*vVb(0,2)*dx*dz + 2*vVb(1,2)*dy*dz;
1835  return dist>0 ? TMath::Sqrt(dist) : -1;
1836 }
1837 
1839  { // check for multi-vertexer pile-up
1840  const int kMinPlpContrib = 5;
1841  const double kMaxPlpChi2 = 5.0;
1842  const double kMinWDist = 15;
1843 
1844  const AliVVertex* vtPrm = 0;
1845  const AliVVertex* vtPlp = 0;
1846 
1847  int nPlp = 0;
1848 
1849  if(!(nPlp=aod->GetNumberOfPileupVerticesTracks()))
1850  return kFALSE;
1851 
1852  vtPrm = aod->GetPrimaryVertex();
1853  if(vtPrm == aod->GetPrimaryVertexSPD())
1854  return kTRUE; // there are pile-up vertices but no primary
1855 
1856  //int bcPrim = vtPrm->GetBC();
1857 
1858  for(int ipl=0;ipl<nPlp;ipl++) {
1859  vtPlp = (const AliVVertex*)aod->GetPileupVertexTracks(ipl);
1860  if (vtPlp->GetNContributors() < kMinPlpContrib) continue;
1861  if (vtPlp->GetChi2perNDF() > kMaxPlpChi2) continue;
1862  //int bcPlp = vtPlp->GetBC();
1863  //if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2)
1864  // return kTRUE; // pile-up from other BC
1865 
1866  double wDst = GetWDist(vtPrm,vtPlp);
1867  if (wDst<kMinWDist) continue;
1868 
1869  return kTRUE; // pile-up: well separated vertices
1870  }
1871 
1872  return kFALSE;
1873  }
1874 
1875 
1876 
1877 
virtual void UserExec(Option_t *option)
double Double_t
Definition: External.C:58
virtual void Terminate(Option_t *)
AliFlowTrackSimple * GetTrack(Int_t i)
AliMultSelection * fMultSelection
input event
double GetWDist(const AliVVertex *v0, const AliVVertex *v1)
TProfile * fHist_v1xV1_ZDN_pTDiff[4][10]
TProfile2D * fHist_znAy_V0_VxVy[90][10]
Double_t Phi() const
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
TProfile2D * fHist_znAx_V0_VxVy[90][10]
Definition: External.C:212
Bool_t plpMV(const AliAODEvent *aod)
Double_t GetCentrality() const
Int_t GetReferenceMultiplicity() const
virtual void GetZDC2Qsub(AliFlowVector *Qarray)
TList * fListZDCQxy
collection of output
TProfile * fHist_v1xV1_ZDN_EtaDiff[4][10]
Bool_t IsPOItype(Int_t poiType) const
const char Option_t
Definition: External.C:48
Double_t Eta() const
bool Bool_t
Definition: External.C:53
TProfile2D * fHist_znCy_V0_VxVy[90][10]
ClassImp(AliAnalysisTaskZDCGainEq) AliAnalysisTaskZDCGainEq
TProfile2D * fHist_znCx_V0_VxVy[90][10]
Int_t NumberOfTracks() const