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