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