AliPhysics  b555aef (b555aef)
AliAnalysisTaskMixedHarmonics.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 
16 /**************************************
17  * analysis task for mixed harmomics *
18  * *
19  * authors: Naomi van der Kolk *
20  * (kolk@nikhef.nl) *
21  * Raimond Snellings *
22  * (snelling@nikhef.nl) *
23  * Ante Bilandzic *
24  * (anteb@nikhef.nl) *
25  * ***********************************/
26 
27 class TFile;
28 class TString;
29 class TList;
30 class AliAnalysisTaskSE;
31 
32 #include "Riostream.h"
33 #include "AliFlowEventSimple.h"
36 
37 #include "AliMultSelection.h"
38 #include "AliAnalysisUtils.h"
39 #include "AliVVertex.h"
40 #include "AliAODTrack.h"
41 #include "AliAODEvent.h"
42 #include "AliAODHeader.h"
43 #include "AliAODVertex.h"
44 #include "TMatrixDSym.h"
45 
46 
47 using std::cout;
48 using std::endl;
50 
51 //================================================================================================================
52 
53 AliAnalysisTaskMixedHarmonics::AliAnalysisTaskMixedHarmonics(const char *name, Bool_t useParticleWeights):
54 AliAnalysisTaskSE(name),
55 fMultSelection(NULL),
56 fAnalysisUtil(NULL),
57 fEvent(NULL),
58 fMH(NULL),
59 fListHistos(NULL),
60 fHarmonic(1),
61 fNoOfMultipicityBins(100),
62 fMultipicityBinWidth(1.),
63 fMinMultiplicity(3.),
64 fOppositeChargesPOI(kFALSE),
65 fEvaluateDifferential3pCorrelator(kFALSE),
66 fCorrectForDetectorEffects(kFALSE),
67 fPrintOnTheScreen(kTRUE),
68 fCalculateVsM(kFALSE),
69 fShowBinLabelsVsM(kFALSE),
70 fUseParticleWeights(useParticleWeights),
71 fUsePhiWeights(kFALSE),
72 fUsePtWeights(kFALSE),
73 fUseEtaWeights(kFALSE),
74 fRejectPileUp(kFALSE),
75 fRejectPileUpTight(kFALSE),
76 fWeightsList(NULL)
77 {
78  // constructor
79  cout<<"AliAnalysisTaskMixedHarmonics::AliAnalysisTaskMixedHarmonics(const char *name, Bool_t useParticleWeights)"<<endl;
80 
81  // Define input and output slots here
82  // Input slot #0 works with an AliFlowEventSimple
83  DefineInput(0, AliFlowEventSimple::Class());
84  // Input slot #1 is needed for the weights input file:
85  if(useParticleWeights)
86  {
87  DefineInput(1, TList::Class());
88  }
89  // Output slot #0 is reserved
90  // Output slot #1 writes into a TList container
91  DefineOutput(1, TList::Class());
92 }
93 
96 fMultSelection(NULL),
97 fAnalysisUtil(NULL),
98 fEvent(NULL),
99 fMH(NULL),
100 fListHistos(NULL),
101 fHarmonic(0),
102 fNoOfMultipicityBins(0),
103 fMultipicityBinWidth(0),
104 fMinMultiplicity(0),
105 fOppositeChargesPOI(kFALSE),
106 fEvaluateDifferential3pCorrelator(kFALSE),
107 fCorrectForDetectorEffects(kFALSE),
108 fPrintOnTheScreen(kFALSE),
109 fCalculateVsM(kFALSE),
110 fShowBinLabelsVsM(kFALSE),
111 fUseParticleWeights(kFALSE),
112 fUsePhiWeights(kFALSE),
113 fUsePtWeights(kFALSE),
114 fUseEtaWeights(kFALSE),
115 fRejectPileUp(kFALSE),
116 fRejectPileUpTight(kFALSE),
117 fWeightsList(NULL)
118 {
119  // Dummy constructor
120  cout<<"AliAnalysisTaskMixedHarmonics::AliAnalysisTaskMixedHarmonics()"<<endl;
121 }
122 
123 //================================================================================================================
124 
126 {
127 
128  fAnalysisUtil = new AliAnalysisUtils();
129  fAnalysisUtil->SetUseMVPlpSelection(kTRUE);
130  fAnalysisUtil->SetUseOutOfBunchPileUp(kTRUE);
131 
132 
133  // Called at every worker node to initialize
134  cout<<"AliAnalysisTaskMixedHarmonics::UserCreateOutputObjects()"<<endl;
135 
136  // Analyser:
138 
139  // Common:
151  {
152  // Pass the flags to class:
156  // Get data from input slot #1 which is used for weights:
157  if(GetNinputs()==2)
158  {
159  fWeightsList = (TList*)GetInputData(1);
160  }
161  // Pass the list with weights to class:
163  }
164 
165  fMH->Init();
166 
167  if(fMH->GetHistList())
168  {
170  // fListHistos->Print();
171  } else
172  {
173  Printf("ERROR: Could not retrieve histogram list (MH, Task::UserCreateOutputObjects()) !!!!");
174  }
175 
176  PostData(1,fListHistos);
177 
178 } // end of void AliAnalysisTaskMixedHarmonics::UserCreateOutputObjects()
179 
180 //================================================================================================================
181 
183 {
184 
185 
186 
187  AliAODEvent *aodEvent = dynamic_cast<AliAODEvent*>(InputEvent());
188 
189  Bool_t kPileupEvent = kFALSE;
190 
191  //cout<<" Run = "<<aodEvent->GetRunNumber();
192 
193  if(fRejectPileUp){
194  kPileupEvent = CheckEventIsPileUp(aodEvent,fRejectPileUpTight);
195  }
196 
197  if(kPileupEvent) return;
198 
199 
200 
201 
202  // main loop (called for each event)
203  fEvent = dynamic_cast<AliFlowEventSimple*>(GetInputData(0));
204 
205  //cout<<" tracks = "<<fEvent->NumberOfTracks()<<endl;
206 
207  // Mixed Harmonics:
208  if(fEvent)
209  {
210  fMH->Make(fEvent);
211  } else
212  {
213  cout<<"WARNING: No input data (MH, Task::UserExec()) !!!!"<<endl;
214  cout<<endl;
215  }
216 
217  PostData(1,fListHistos);
218 }
219 
220 //================================================================================================================
221 
223 {
224  //accessing the merged output list:
225  fListHistos = (TList*)GetOutputData(1);
226 
228 
229  if(fListHistos)
230  {
232  fMH->Finish();
233  PostData(1,fListHistos);
234  } else
235  {
236  cout<<" WARNING: histogram list pointer is empty (MH, Task::Terminate()) !!!!"<<endl;
237  cout<<endl;
238  }
239 
240 } // end of void AliAnalysisTaskMixedHarmonics::Terminate(Option_t *)
241 
242 
243 
244 
245 
246 
248 
249  Bool_t BisPileup=kFALSE;
250 
251  Double_t centrV0M=300;
252  Double_t centrCL1=300;
253  Double_t centrCL0=300;
254  Double_t centrTRK=300;
255 
256  fMultSelection = (AliMultSelection*) InputEvent()->FindListObject("MultSelection");
257 
258  centrV0M = fMultSelection->GetMultiplicityPercentile("V0M");
259  centrCL1 = fMultSelection->GetMultiplicityPercentile("CL1");
260  centrCL0 = fMultSelection->GetMultiplicityPercentile("CL0");
261  centrTRK = fMultSelection->GetMultiplicityPercentile("TRK");
262 
263  if(fRejectPileUp && InputEvent()) {
264  //------------ pileup for 2015 data ----------
265  BisPileup=kFALSE;
266 
267  //-- pile-up a la Dobrin for LHC15o -----
268  if(PileUpMultiVertex(faod)) {
269  //fPileUpCount->Fill(0.5);
270  BisPileup=kTRUE;
271  }
272  Int_t isPileup = faod->IsPileupFromSPD(3);
273  if(isPileup != 0) {
274  //fPileUpCount->Fill(1.5);
275  BisPileup=kTRUE;
276  }
277  if(((AliAODHeader*)faod->GetHeader())->GetRefMultiplicityComb08() < 0) {
278  //fPileUpCount->Fill(2.5);
279  BisPileup=kTRUE;
280  }
281  if(faod->IsIncompleteDAQ()) {
282  //fPileUpCount->Fill(3.5);
283  BisPileup=kTRUE;
284  }
285  if(fabs(centrV0M-centrCL1)>7.5) {
286  //fPileUpCount->Fill(4.5);
287  BisPileup=kTRUE;
288  }
289 
290  // check vertex consistency
291  const AliAODVertex* vtTrc = faod->GetPrimaryVertex();
292  const AliAODVertex* vtSPD = faod->GetPrimaryVertexSPD();
293 
294  if(vtTrc->GetNContributors() < 2 || vtSPD->GetNContributors()<1) {
295  //fPileUpCount->Fill(5.5);
296  BisPileup=kTRUE;
297  }
298 
299  double covTrc[6], covSPD[6];
300  vtTrc->GetCovarianceMatrix(covTrc);
301  vtSPD->GetCovarianceMatrix(covSPD);
302 
303  double dz = vtTrc->GetZ() - vtSPD->GetZ();
304 
305  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
306  double errTrc = TMath::Sqrt(covTrc[5]);
307  double nsigTot = dz/errTot;
308  double nsigTrc = dz/errTrc;
309 
310  if(TMath::Abs(dz)>0.2 || TMath::Abs(nsigTot)>10 || TMath::Abs(nsigTrc)>20) {
311  //fPileUpCount->Fill(6.5);
312  BisPileup=kTRUE;
313  }
314 
315  //cuts on tracks
316  const Int_t nTracks = faod->GetNumberOfTracks();
317  Int_t multEsd = ((AliAODHeader*)faod->GetHeader())->GetNumberOfESDTracks();
318 
319  //Int_t multTrk = 0;
320  //Int_t multTrkBefC = 0;
321  //Int_t multTrkTOFBefC = 0;
322  Int_t multTPC = 0;
323 
324  for(Int_t it = 0; it < nTracks; it++) {
325  AliAODTrack* aodTrk = (AliAODTrack*)faod->GetTrack(it);
326  if(!aodTrk) {
327  delete aodTrk;
328  continue;
329  }
330  //if(aodTrk->TestFilterBit(32)){
331  // multTrkBefC++;
332  // if(TMath::Abs(aodTrk->GetTOFsignalDz()) <= 10. && aodTrk->GetTOFsignal() >= 12000. && aodTrk->GetTOFsignal() <= 25000.)
333  // multTrkTOFBefC++;
334  // if((TMath::Abs(aodTrk->Eta()) < 0.8) && (aodTrk->GetTPCNcls() >= 70) && (aodTrk->Pt() >= 0.2) && (aodTrk->Pt() < 20.))
335  // multTrk++;
336  //}
337  if(aodTrk->TestFilterBit(128))
338  multTPC++;
339  } // end of for AOD track loop
340 
341  Double_t multTPCn = multTPC;
342  Double_t multEsdn = multEsd;
343  Double_t multESDTPCDif = multEsdn - multTPCn*3.38;
344 
345  if(multESDTPCDif > 15000.){
346  //fPileUpCount->Fill(7.5);
347  BisPileup=kTRUE;
348  }
349  else if(bPileUpTight) {
350  if(multESDTPCDif > 700.) {
351  //fPileUpCount->Fill(8.5);
352  BisPileup=kTRUE;
353  }
354  if(BisPileup==kFALSE) {
355  if(!fMultSelection->GetThisEventIsNotPileup()) BisPileup=kTRUE;
356  if(!fMultSelection->GetThisEventIsNotPileupMV()) BisPileup=kTRUE;
357  if(!fMultSelection->GetThisEventIsNotPileupInMultBins()) BisPileup=kTRUE;
358  if(!fMultSelection->GetThisEventHasNoInconsistentVertices()) BisPileup=kTRUE;
359  if(!fMultSelection->GetThisEventPassesTrackletVsCluster()) BisPileup=kTRUE;
360  if(!fMultSelection->GetThisEventIsNotIncompleteDAQ()) BisPileup=kTRUE;
361  if(!fMultSelection->GetThisEventHasGoodVertex2016()) BisPileup=kTRUE;
362  //if(BisPileup) fPileUpCount->Fill(9.5);
363  }
364  }
365  }
366 
367  return BisPileup;
368 }//-------pile up function ------
369 
370 
371 
373  { // check for multi-vertexer pile-up
374  const int kMinPlpContrib = 5;
375  const double kMaxPlpChi2 = 5.0;
376  const double kMinWDist = 15;
377 
378  const AliVVertex* vtPrm = 0;
379  const AliVVertex* vtPlp = 0;
380 
381  int nPlp = 0;
382 
383  if(!(nPlp=faod->GetNumberOfPileupVerticesTracks()))
384  return kFALSE;
385 
386  vtPrm = faod->GetPrimaryVertex();
387  if(vtPrm == faod->GetPrimaryVertexSPD())
388  return kTRUE; // there are pile-up vertices but no primary
389 
390  //int bcPrim = vtPrm->GetBC();
391 
392  for(int ipl=0;ipl<nPlp;ipl++) {
393  vtPlp = (const AliVVertex*)faod->GetPileupVertexTracks(ipl);
394  if (vtPlp->GetNContributors() < kMinPlpContrib) continue;
395  if (vtPlp->GetChi2perNDF() > kMaxPlpChi2) continue;
396  //int bcPlp = vtPlp->GetBC();
397  //if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2)
398  // return kTRUE; // pile-up from other BC
399 
400  double wDst = GetWDist(vtPrm,vtPlp);
401  if (wDst<kMinWDist) continue;
402 
403  return kTRUE; // pile-up: well separated vertices
404  }
405  return kFALSE;
406 }
407 
408 double AliAnalysisTaskMixedHarmonics::GetWDist(const AliVVertex* v0, const AliVVertex* v1)
409 {
410  // calculate sqrt of weighted distance to other vertex
411  if (!v0 || !v1) {
412  AliDebug(2,"\n\n ::GetWDist => One of vertices is not valid\n\n");
413  return 0;
414  }
415  static TMatrixDSym vVb(3);
416  double dist = -1;
417  double dx = v0->GetX()-v1->GetX();
418  double dy = v0->GetY()-v1->GetY();
419  double dz = v0->GetZ()-v1->GetZ();
420  double cov0[6],cov1[6];
421  v0->GetCovarianceMatrix(cov0);
422  v1->GetCovarianceMatrix(cov1);
423  vVb(0,0) = cov0[0]+cov1[0];
424  vVb(1,1) = cov0[2]+cov1[2];
425  vVb(2,2) = cov0[5]+cov1[5];
426  vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
427  vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
428  vVb.InvertFast();
429  if (!vVb.IsValid()) {
430  AliDebug(2,"Singular Matrix\n");
431  return dist;
432  }
433  dist = vVb(0,0)*dx*dx + vVb(1,1)*dy*dy + vVb(2,2)*dz*dz
434  + 2*vVb(0,1)*dx*dy + 2*vVb(0,2)*dx*dz + 2*vVb(1,2)*dy*dz;
435  return dist>0 ? TMath::Sqrt(dist) : -1;
436 }
437 
438 
439 
440 
441 
442 
443 
444 
445 
446 
447 
448 
449 
450 
451 
452 
453 
454 
455 
virtual void GetOutputHistograms(TList *outputListHistos)
double Double_t
Definition: External.C:58
AliFlowEventSimple * fEvent
Event selection.
AliFlowAnalysisWithMixedHarmonics * fMH
virtual void Make(AliFlowEventSimple *anEvent)
int Int_t
Definition: External.C:63
AliAnalysisUtils * fAnalysisUtil
MultSelection (RUN2 centrality estimator)
virtual void UserExec(Option_t *option)
Bool_t CheckEventIsPileUp(AliAODEvent *faod, Bool_t bPileUpTight=kFALSE)
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
double GetWDist(const AliVVertex *v0, const AliVVertex *v1)
Bool_t PileUpMultiVertex(const AliAODEvent *faod)