AliPhysics  59e0e03 (59e0e03)
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 fFillQAHistograms(kFALSE),
77 fTPCvsGlobalTrkBefore(NULL),
78 fTPCvsGlobalTrkAfter(NULL),
79 fTPCvsESDTrk(NULL),
80 fWeightsList(NULL)
81 {
82  // constructor
83  cout<<"AliAnalysisTaskMixedHarmonics::AliAnalysisTaskMixedHarmonics(const char *name, Bool_t useParticleWeights)"<<endl;
84 
85  // Define input and output slots here
86  // Input slot #0 works with an AliFlowEventSimple
87  DefineInput(0, AliFlowEventSimple::Class());
88  // Input slot #1 is needed for the weights input file:
89  if(useParticleWeights)
90  {
91  DefineInput(1, TList::Class());
92  }
93  // Output slot #0 is reserved
94  // Output slot #1 writes into a TList container
95  DefineOutput(1, TList::Class());
96 }
97 
100 fMultSelection(NULL),
101 fAnalysisUtil(NULL),
102 fEvent(NULL),
103 fMH(NULL),
104 fListHistos(NULL),
105 fHarmonic(0),
106 fNoOfMultipicityBins(0),
107 fMultipicityBinWidth(0),
108 fMinMultiplicity(0),
109 fOppositeChargesPOI(kFALSE),
110 fEvaluateDifferential3pCorrelator(kFALSE),
111 fCorrectForDetectorEffects(kFALSE),
112 fPrintOnTheScreen(kFALSE),
113 fCalculateVsM(kFALSE),
114 fShowBinLabelsVsM(kFALSE),
115 fUseParticleWeights(kFALSE),
116 fUsePhiWeights(kFALSE),
117 fUsePtWeights(kFALSE),
118 fUseEtaWeights(kFALSE),
119 fRejectPileUp(kFALSE),
120 fRejectPileUpTight(kFALSE),
121 fFillQAHistograms(kFALSE),
122 fTPCvsGlobalTrkBefore(NULL),
123 fTPCvsGlobalTrkAfter(NULL),
124 fTPCvsESDTrk(NULL),
125 fWeightsList(NULL)
126 {
127  // Dummy constructor
128  cout<<"AliAnalysisTaskMixedHarmonics::AliAnalysisTaskMixedHarmonics()"<<endl;
129 }
130 
131 //================================================================================================================
132 
134 {
135 
136  fAnalysisUtil = new AliAnalysisUtils();
137  fAnalysisUtil->SetUseMVPlpSelection(kTRUE);
138  fAnalysisUtil->SetUseOutOfBunchPileUp(kTRUE);
139 
140 
141  // Called at every worker node to initialize
142  cout<<"AliAnalysisTaskMixedHarmonics::UserCreateOutputObjects()"<<endl;
143 
144  // Analyser:
146 
147  // Common:
159  {
160  // Pass the flags to class:
164  // Get data from input slot #1 which is used for weights:
165  if(GetNinputs()==2)
166  {
167  fWeightsList = (TList*)GetInputData(1);
168  }
169  // Pass the list with weights to class:
171  }
172 
173  fMH->Init();
174 
175  if(fMH->GetHistList())
176  {
178  // fListHistos->Print();
179  } else
180  {
181  Printf("ERROR: Could not retrieve histogram list (MH, Task::UserCreateOutputObjects()) !!!!");
182  }
183 
184 
185  //Add QA histograms:
186  fTPCvsGlobalTrkBefore = new TH2F("fTPCvsGlobalTrkBefore","Global(Fb32) vs TPC(FB128)",250,0,5000,250,0,5000);
188  fTPCvsGlobalTrkAfter = new TH2F("fTPCvsGlobalTrkAfter","Global(Fb32) vs TPC(FB128)",250,0,5000,250,0,5000);
190 
191 
192  fTPCvsESDTrk = new TH2F("fTPCvsESDTrk","ESDTrk vs TPC(FB128)",1000,0,20000,250,0,5000);
194 
195 
196  PostData(1,fListHistos);
197 
198 } // end of void AliAnalysisTaskMixedHarmonics::UserCreateOutputObjects()
199 
200 //================================================================================================================
201 
203 {
204 
205 
206 
207  AliAODEvent *aodEvent = dynamic_cast<AliAODEvent*>(InputEvent());
208 
209  Bool_t kPileupEvent = kFALSE;
210 
211  //cout<<" Run = "<<aodEvent->GetRunNumber();
212 
214  kPileupEvent = CheckEventIsPileUp(aodEvent);
215  }
216 
217  if(fRejectPileUp && kPileupEvent) return;
218 
219 
220 
221 
222  // main loop (called for each event)
223  fEvent = dynamic_cast<AliFlowEventSimple*>(GetInputData(0));
224 
225  //cout<<" tracks = "<<fEvent->NumberOfTracks()<<endl;
226 
227  // Mixed Harmonics:
228  if(fEvent)
229  {
230  fMH->Make(fEvent);
231  }
232  else
233  {
234  cout<<"WARNING: No input data (MH, Task::UserExec()) !!!!"<<endl;
235  cout<<endl;
236  }
237 
238  PostData(1,fListHistos);
239 
240 }
241 
242 
243 
244 
245 
246 
247 
248 //================================================================================================================
249 
251 {
252  //accessing the merged output list:
253  fListHistos = (TList*)GetOutputData(1);
254 
256 
257  if(fListHistos)
258  {
260  fMH->Finish();
261  PostData(1,fListHistos);
262  } else
263  {
264  cout<<" WARNING: histogram list pointer is empty (MH, Task::Terminate()) !!!!"<<endl;
265  cout<<endl;
266  }
267 
268 } // end of void AliAnalysisTaskMixedHarmonics::Terminate(Option_t *)
269 
270 
271 
272 
273 //----- PileUp removal function -------
274 
276 
277  Bool_t BisPileup=kFALSE;
278 
279  Double_t centrV0M=300;
280  Double_t centrCL1=300;
281  Double_t centrCL0=300;
282  Double_t centrTRK=300;
283 
284  fMultSelection = (AliMultSelection*) InputEvent()->FindListObject("MultSelection");
285 
286  if(!fMultSelection) {
287  printf("\n\n **WARNING** ::UserExec() AliMultSelection object not found.\n\n");
288  exit(1);
289  }
290  centrV0M = fMultSelection->GetMultiplicityPercentile("V0M");
291  centrCL1 = fMultSelection->GetMultiplicityPercentile("CL1");
292  centrCL0 = fMultSelection->GetMultiplicityPercentile("CL0");
293  centrTRK = fMultSelection->GetMultiplicityPercentile("TRK");
294 
295  //-- pile-up a la Dobrin for LHC15o -----
296  if(PileUpMultiVertex(faod)) {
297  //fPileUpCount->Fill(0.5);
298  BisPileup=kTRUE;
299  }
300  Int_t isPileup = faod->IsPileupFromSPD(3);
301  if(isPileup != 0) {
302  //fPileUpCount->Fill(1.5);
303  BisPileup=kTRUE;
304  }
305  if(((AliAODHeader*)faod->GetHeader())->GetRefMultiplicityComb08() < 0) {
306  //fPileUpCount->Fill(2.5);
307  BisPileup=kTRUE;
308  }
309  if(faod->IsIncompleteDAQ()) {
310  //fPileUpCount->Fill(3.5);
311  BisPileup=kTRUE;
312  }
313  if(fabs(centrV0M-centrCL1)> 5.0) {//default: 7.5
314  //fPileUpCount->Fill(4.5);
315  BisPileup=kTRUE;
316  }
317 
318  // check vertex consistency
319  const AliAODVertex* vtTrc = faod->GetPrimaryVertex();
320  const AliAODVertex* vtSPD = faod->GetPrimaryVertexSPD();
321 
322  if(vtTrc->GetNContributors() < 2 || vtSPD->GetNContributors()<1) {
323  //fPileUpCount->Fill(5.5);
324  BisPileup=kTRUE;
325  }
326 
327  double covTrc[6], covSPD[6];
328  vtTrc->GetCovarianceMatrix(covTrc);
329  vtSPD->GetCovarianceMatrix(covSPD);
330 
331  double dz = vtTrc->GetZ() - vtSPD->GetZ();
332 
333  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
334  double errTrc = TMath::Sqrt(covTrc[5]);
335  double nsigTot = dz/errTot;
336  double nsigTrc = dz/errTrc;
337 
338  if(TMath::Abs(dz)>0.2 || TMath::Abs(nsigTot)>10 || TMath::Abs(nsigTrc)>20) {
339  //fPileUpCount->Fill(6.5);
340  BisPileup=kTRUE;
341  }
342 
343  //cuts on tracks
344  //Int_t multTrk = 0;
345  //Int_t multTrkBefC = 0;
346  //Int_t multTrkTOFBefC = 0;
347 
348  Int_t multTPC = 0;
349  Int_t multITSfb96 = 0;
350  Int_t multITSfb32 = 0;
351 
352  Int_t multTPCFE = 0;
353  Int_t multGlobal = 0;
354  Int_t multTPCuncut = 0;
355 
356  Int_t multEsd = ((AliAODHeader*)faod->GetHeader())->GetNumberOfESDTracks();
357 
358  const Int_t nTracks = faod->GetNumberOfTracks();
359 
360  for(Int_t iTracks = 0; iTracks < nTracks; iTracks++) {
361  //AliNanoAODTrack* track = dynamic_cast<AliNanoAODTrack*>(faod->GetTrack(iTracks));
362  AliAODTrack* track = (AliAODTrack*)faod->GetTrack(iTracks);
363  if(!track) continue;
364  //---------- old method -----------
365  if(track->TestFilterBit(128))
366  multTPC++;
367  if(track->TestFilterBit(96))
368  multITSfb96++;
369  //----------------------------------
370  if(track->TestFilterBit(1)) multTPCuncut++;
371  if(track->TestFilterBit(32)) multITSfb32++;
372 
373 
374  if(track->Pt()<0.2 || track->Pt()>5.0 || TMath::Abs(track->Eta())>0.8 || track->GetTPCNcls()<70 || track->GetTPCsignal()<10.0)
375  continue;
376  if(track->TestFilterBit(1) && track->Chi2perNDF()>0.2) multTPCFE++;
377  if(!track->TestFilterBit(16) || track->Chi2perNDF()<0.1) continue;
378 
379  Double_t b[2] = {-99., -99.};
380  Double_t bCov[3] = {-99., -99., -99.};
381 
382  AliAODTrack copy(*track);
383  Double_t magField = faod->GetMagneticField();
384 
385  if(magField!=0){
386  if(track->PropagateToDCA(faod->GetPrimaryVertex(), magField, 100., b, bCov) && TMath::Abs(b[0]) < 0.3 && TMath::Abs(b[1]) < 0.3) multGlobal++;
387  }
388  }
389 
390  Double_t multTPCn = multTPC;
391  Double_t multEsdn = multEsd;
392 
393  //fixed for test:
394  Double_t fPileUpSlopeParm = 3.55;
395  Double_t fPileUpConstParm = 90;
396 
397  Double_t multESDTPCDif = multEsdn - fPileUpSlopeParm*multTPCn;
398  //Double_t multTPCGlobDif = multTPCFE - fPileUpSlopeParm*multGlobal;
399 
400  if(fFillQAHistograms){
401  //fGlobalTracks->Fill(multGlobal);
402  fTPCvsGlobalTrkBefore->Fill(multITSfb32,multTPC);
403  fTPCvsESDTrk->Fill(multEsd,multTPC);
404  }
405 
406  /*if(multESDTPCDif > (fRejectPileUpTight?700.:15000.)) {
407  //fPileUpCount->Fill(7.5);
408  BisPileup=kTRUE;
409  }*/
410  /*if(multESDTPCDif > 15000.) { //default: 15000
411  //fPileUpCount->Fill(7.5);
412  BisPileup=kTRUE;
413  }*/
414 
415 
416 
417  if(fRejectPileUp) {
418  if(multESDTPCDif > fPileUpConstParm) {
419  //fPileUpCount->Fill(7.5);
420  BisPileup=kTRUE;
421  }
422  if(BisPileup==kFALSE) {
423  if(!fMultSelection->GetThisEventIsNotPileup()) BisPileup=kTRUE;
424  if(!fMultSelection->GetThisEventIsNotPileupMV()) BisPileup=kTRUE;
425  if(!fMultSelection->GetThisEventIsNotPileupInMultBins()) BisPileup=kTRUE;
426  if(!fMultSelection->GetThisEventHasNoInconsistentVertices()) BisPileup=kTRUE;
427  if(!fMultSelection->GetThisEventPassesTrackletVsCluster()) BisPileup=kTRUE;
428  if(!fMultSelection->GetThisEventIsNotIncompleteDAQ()) BisPileup=kTRUE;
429  if(!fMultSelection->GetThisEventHasGoodVertex2016()) BisPileup=kTRUE;
430  //if(BisPileup) fPileUpCount->Fill(9.5);
431  }
432  }
433 
434  if(!BisPileup){
435  fTPCvsGlobalTrkAfter->Fill(multITSfb32,multTPC);
436  }
437 
438 
439  return BisPileup;
440 } //-------pile up function ------
441 
442 
443 
444 
445 
446 
447 
448 
450  { // check for multi-vertexer pile-up
451  const int kMinPlpContrib = 5;
452  const double kMaxPlpChi2 = 5.0;
453  const double kMinWDist = 15;
454 
455  const AliVVertex* vtPrm = 0;
456  const AliVVertex* vtPlp = 0;
457 
458  int nPlp = 0;
459 
460  if(!(nPlp=faod->GetNumberOfPileupVerticesTracks()))
461  return kFALSE;
462 
463  vtPrm = faod->GetPrimaryVertex();
464  if(vtPrm == faod->GetPrimaryVertexSPD())
465  return kTRUE; // there are pile-up vertices but no primary
466 
467  //int bcPrim = vtPrm->GetBC();
468 
469  for(int ipl=0;ipl<nPlp;ipl++) {
470  vtPlp = (const AliVVertex*)faod->GetPileupVertexTracks(ipl);
471  if (vtPlp->GetNContributors() < kMinPlpContrib) continue;
472  if (vtPlp->GetChi2perNDF() > kMaxPlpChi2) continue;
473  //int bcPlp = vtPlp->GetBC();
474  //if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2)
475  // return kTRUE; // pile-up from other BC
476 
477  double wDst = GetWDist(vtPrm,vtPlp);
478  if (wDst<kMinWDist) continue;
479 
480  return kTRUE; // pile-up: well separated vertices
481  }
482  return kFALSE;
483 }
484 
485 double AliAnalysisTaskMixedHarmonics::GetWDist(const AliVVertex* v0, const AliVVertex* v1)
486 {
487  // calculate sqrt of weighted distance to other vertex
488  if (!v0 || !v1) {
489  AliDebug(2,"\n\n ::GetWDist => One of vertices is not valid\n\n");
490  return 0;
491  }
492  static TMatrixDSym vVb(3);
493  double dist = -1;
494  double dx = v0->GetX()-v1->GetX();
495  double dy = v0->GetY()-v1->GetY();
496  double dz = v0->GetZ()-v1->GetZ();
497  double cov0[6],cov1[6];
498  v0->GetCovarianceMatrix(cov0);
499  v1->GetCovarianceMatrix(cov1);
500  vVb(0,0) = cov0[0]+cov1[0];
501  vVb(1,1) = cov0[2]+cov1[2];
502  vVb(2,2) = cov0[5]+cov1[5];
503  vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
504  vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
505  vVb.InvertFast();
506  if (!vVb.IsValid()) {
507  AliDebug(2,"Singular Matrix\n");
508  return dist;
509  }
510  dist = vVb(0,0)*dx*dx + vVb(1,1)*dy*dy + vVb(2,2)*dz*dz
511  + 2*vVb(0,1)*dx*dy + 2*vVb(0,2)*dx*dz + 2*vVb(1,2)*dy*dz;
512  return dist>0 ? TMath::Sqrt(dist) : -1;
513 }
514 
515 
516 
517 
518 
519 
520 
521 
522 
523 
524 
525 
526 
527 
528 
529 
530 
531 
532 
virtual void GetOutputHistograms(TList *outputListHistos)
double Double_t
Definition: External.C:58
Definition: External.C:236
AliFlowEventSimple * fEvent
Event selection.
AliFlowAnalysisWithMixedHarmonics * fMH
virtual void Make(AliFlowEventSimple *anEvent)
TH2F * fTPCvsESDTrk
Global vs TPC tracks for QA.
int Int_t
Definition: External.C:63
AliAnalysisUtils * fAnalysisUtil
MultSelection (RUN2 centrality estimator)
virtual void UserExec(Option_t *option)
TList * fWeightsList
Global vs TPC tracks for QA.
const char Option_t
Definition: External.C:48
TH2F * fTPCvsGlobalTrkAfter
Global vs TPC tracks for QA.
bool Bool_t
Definition: External.C:53
double GetWDist(const AliVVertex *v0, const AliVVertex *v1)
Bool_t PileUpMultiVertex(const AliAODEvent *faod)