AliPhysics  d565ceb (d565ceb)
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(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  } else
232  {
233  cout<<"WARNING: No input data (MH, Task::UserExec()) !!!!"<<endl;
234  cout<<endl;
235  }
236 
237  PostData(1,fListHistos);
238 }
239 
240 //================================================================================================================
241 
243 {
244  //accessing the merged output list:
245  fListHistos = (TList*)GetOutputData(1);
246 
248 
249  if(fListHistos)
250  {
252  fMH->Finish();
253  PostData(1,fListHistos);
254  } else
255  {
256  cout<<" WARNING: histogram list pointer is empty (MH, Task::Terminate()) !!!!"<<endl;
257  cout<<endl;
258  }
259 
260 } // end of void AliAnalysisTaskMixedHarmonics::Terminate(Option_t *)
261 
262 
263 
264 
265 //----- PileUp removal function -------
266 
268 
269  Bool_t BisPileup=kFALSE;
270 
271  Double_t centrV0M=300;
272  Double_t centrCL1=300;
273  Double_t centrCL0=300;
274  Double_t centrTRK=300;
275 
276  fMultSelection = (AliMultSelection*) InputEvent()->FindListObject("MultSelection");
277 
278  if(!fMultSelection) {
279  printf("\n\n **WARNING** ::UserExec() AliMultSelection object not found.\n\n");
280  exit(1);
281  }
282  centrV0M = fMultSelection->GetMultiplicityPercentile("V0M");
283  centrCL1 = fMultSelection->GetMultiplicityPercentile("CL1");
284  centrCL0 = fMultSelection->GetMultiplicityPercentile("CL0");
285  centrTRK = fMultSelection->GetMultiplicityPercentile("TRK");
286 
287  //-- pile-up a la Dobrin for LHC15o -----
288  if(PileUpMultiVertex(faod)) {
289  //fPileUpCount->Fill(0.5);
290  BisPileup=kTRUE;
291  }
292  Int_t isPileup = faod->IsPileupFromSPD(3);
293  if(isPileup != 0) {
294  //fPileUpCount->Fill(1.5);
295  BisPileup=kTRUE;
296  }
297  if(((AliAODHeader*)faod->GetHeader())->GetRefMultiplicityComb08() < 0) {
298  //fPileUpCount->Fill(2.5);
299  BisPileup=kTRUE;
300  }
301  if(faod->IsIncompleteDAQ()) {
302  //fPileUpCount->Fill(3.5);
303  BisPileup=kTRUE;
304  }
305  if(fabs(centrV0M-centrCL1)> 5.0) {//default: 7.5
306  //fPileUpCount->Fill(4.5);
307  BisPileup=kTRUE;
308  }
309 
310  // check vertex consistency
311  const AliAODVertex* vtTrc = faod->GetPrimaryVertex();
312  const AliAODVertex* vtSPD = faod->GetPrimaryVertexSPD();
313 
314  if(vtTrc->GetNContributors() < 2 || vtSPD->GetNContributors()<1) {
315  //fPileUpCount->Fill(5.5);
316  BisPileup=kTRUE;
317  }
318 
319  double covTrc[6], covSPD[6];
320  vtTrc->GetCovarianceMatrix(covTrc);
321  vtSPD->GetCovarianceMatrix(covSPD);
322 
323  double dz = vtTrc->GetZ() - vtSPD->GetZ();
324 
325  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
326  double errTrc = TMath::Sqrt(covTrc[5]);
327  double nsigTot = dz/errTot;
328  double nsigTrc = dz/errTrc;
329 
330  if(TMath::Abs(dz)>0.2 || TMath::Abs(nsigTot)>10 || TMath::Abs(nsigTrc)>20) {
331  //fPileUpCount->Fill(6.5);
332  BisPileup=kTRUE;
333  }
334 
335  //cuts on tracks
336  //Int_t multTrk = 0;
337  //Int_t multTrkBefC = 0;
338  //Int_t multTrkTOFBefC = 0;
339 
340  Int_t multTPC = 0;
341  Int_t multITSfb96 = 0;
342  Int_t multITSfb32 = 0;
343 
344  Int_t multTPCFE = 0;
345  Int_t multGlobal = 0;
346  Int_t multTPCuncut = 0;
347 
348  Int_t multEsd = ((AliAODHeader*)faod->GetHeader())->GetNumberOfESDTracks();
349 
350  const Int_t nTracks = faod->GetNumberOfTracks();
351 
352  for(Int_t iTracks = 0; iTracks < nTracks; iTracks++) {
353  //AliNanoAODTrack* track = dynamic_cast<AliNanoAODTrack*>(faod->GetTrack(iTracks));
354  AliAODTrack* track = (AliAODTrack*)faod->GetTrack(iTracks);
355  if(!track) continue;
356  //---------- old method -----------
357  if(track->TestFilterBit(128))
358  multTPC++;
359  if(track->TestFilterBit(96))
360  multITSfb96++;
361  //----------------------------------
362  if(track->TestFilterBit(1)) multTPCuncut++;
363  if(track->TestFilterBit(32)) multITSfb32++;
364 
365 
366  if(track->Pt()<0.2 || track->Pt()>5.0 || TMath::Abs(track->Eta())>0.8 || track->GetTPCNcls()<70 || track->GetTPCsignal()<10.0)
367  continue;
368  if(track->TestFilterBit(1) && track->Chi2perNDF()>0.2) multTPCFE++;
369  if(!track->TestFilterBit(16) || track->Chi2perNDF()<0.1) continue;
370 
371  Double_t b[2] = {-99., -99.};
372  Double_t bCov[3] = {-99., -99., -99.};
373 
374  AliAODTrack copy(*track);
375  Double_t magField = faod->GetMagneticField();
376 
377  if(magField!=0){
378  if(track->PropagateToDCA(faod->GetPrimaryVertex(), magField, 100., b, bCov) && TMath::Abs(b[0]) < 0.3 && TMath::Abs(b[1]) < 0.3) multGlobal++;
379  }
380  }
381 
382 
383  /*
384  for(Int_t it = 0; it < nTracks; it++) {
385  AliAODTrack* aodTrk = (AliAODTrack*)faod->GetTrack(it);
386  if(!aodTrk) {
387  delete aodTrk;
388  continue;
389  }
390  //if(aodTrk->TestFilterBit(32)){
391  // multTrkBefC++;
392  // if(TMath::Abs(aodTrk->GetTOFsignalDz()) <= 10. && aodTrk->GetTOFsignal() >= 12000. && aodTrk->GetTOFsignal() <= 25000.)
393  // multTrkTOFBefC++;
394  // if((TMath::Abs(aodTrk->Eta()) < 0.8) && (aodTrk->GetTPCNcls() >= 70) && (aodTrk->Pt() >= 0.2) && (aodTrk->Pt() < 20.))
395  // multTrk++;
396  //}
397  if(aodTrk->TestFilterBit(128))
398  multTPC++;
399  if(aodTrk->TestFilterBit(96))
400  multITSfb96++;
401  } // end of for AOD track loop
402 
403  */
404 
405  Double_t multTPCn = multTPC;
406  Double_t multEsdn = multEsd;
407 
408  //fixed for test:
409  Double_t fPileUpSlopeParm = 3.51;
410  Double_t fPileUpConstParm = 50;
411 
412  Double_t multESDTPCDif = multEsdn - fPileUpSlopeParm*multTPCn;
413  //Double_t multTPCGlobDif = multTPCFE - fPileUpSlopeParm*multGlobal;
414 
415  if(fFillQAHistograms){
416  //fGlobalTracks->Fill(multGlobal);
417  fTPCvsGlobalTrkBefore->Fill(multITSfb32,multTPC);
418  fTPCvsESDTrk->Fill(multEsd,multTPC);
419  }
420 
421  /*if(multESDTPCDif > (fRejectPileUpTight?700.:15000.)) {
422  //fPileUpCount->Fill(7.5);
423  BisPileup=kTRUE;
424  }*/
425  /*if(multESDTPCDif > 15000.) { //default: 15000
426  //fPileUpCount->Fill(7.5);
427  BisPileup=kTRUE;
428  }*/
429 
430 
431 
432  if(fRejectPileUp) {
433  if(multESDTPCDif > fPileUpConstParm) {
434  //fPileUpCount->Fill(7.5);
435  BisPileup=kTRUE;
436  }
437  if(BisPileup==kFALSE) {
438  if(!fMultSelection->GetThisEventIsNotPileup()) BisPileup=kTRUE;
439  if(!fMultSelection->GetThisEventIsNotPileupMV()) BisPileup=kTRUE;
440  if(!fMultSelection->GetThisEventIsNotPileupInMultBins()) BisPileup=kTRUE;
441  if(!fMultSelection->GetThisEventHasNoInconsistentVertices()) BisPileup=kTRUE;
442  if(!fMultSelection->GetThisEventPassesTrackletVsCluster()) BisPileup=kTRUE;
443  if(!fMultSelection->GetThisEventIsNotIncompleteDAQ()) BisPileup=kTRUE;
444  if(!fMultSelection->GetThisEventHasGoodVertex2016()) BisPileup=kTRUE;
445  //if(BisPileup) fPileUpCount->Fill(9.5);
447  fTPCvsGlobalTrkAfter->Fill(multITSfb32,multTPC);
448  }
449  }
450 
451 
452  return BisPileup;
453 } //-------pile up function ------
454 
455 
456 
457 
458 
459 
460 
461 
463  { // check for multi-vertexer pile-up
464  const int kMinPlpContrib = 5;
465  const double kMaxPlpChi2 = 5.0;
466  const double kMinWDist = 15;
467 
468  const AliVVertex* vtPrm = 0;
469  const AliVVertex* vtPlp = 0;
470 
471  int nPlp = 0;
472 
473  if(!(nPlp=faod->GetNumberOfPileupVerticesTracks()))
474  return kFALSE;
475 
476  vtPrm = faod->GetPrimaryVertex();
477  if(vtPrm == faod->GetPrimaryVertexSPD())
478  return kTRUE; // there are pile-up vertices but no primary
479 
480  //int bcPrim = vtPrm->GetBC();
481 
482  for(int ipl=0;ipl<nPlp;ipl++) {
483  vtPlp = (const AliVVertex*)faod->GetPileupVertexTracks(ipl);
484  if (vtPlp->GetNContributors() < kMinPlpContrib) continue;
485  if (vtPlp->GetChi2perNDF() > kMaxPlpChi2) continue;
486  //int bcPlp = vtPlp->GetBC();
487  //if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2)
488  // return kTRUE; // pile-up from other BC
489 
490  double wDst = GetWDist(vtPrm,vtPlp);
491  if (wDst<kMinWDist) continue;
492 
493  return kTRUE; // pile-up: well separated vertices
494  }
495  return kFALSE;
496 }
497 
498 double AliAnalysisTaskMixedHarmonics::GetWDist(const AliVVertex* v0, const AliVVertex* v1)
499 {
500  // calculate sqrt of weighted distance to other vertex
501  if (!v0 || !v1) {
502  AliDebug(2,"\n\n ::GetWDist => One of vertices is not valid\n\n");
503  return 0;
504  }
505  static TMatrixDSym vVb(3);
506  double dist = -1;
507  double dx = v0->GetX()-v1->GetX();
508  double dy = v0->GetY()-v1->GetY();
509  double dz = v0->GetZ()-v1->GetZ();
510  double cov0[6],cov1[6];
511  v0->GetCovarianceMatrix(cov0);
512  v1->GetCovarianceMatrix(cov1);
513  vVb(0,0) = cov0[0]+cov1[0];
514  vVb(1,1) = cov0[2]+cov1[2];
515  vVb(2,2) = cov0[5]+cov1[5];
516  vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
517  vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
518  vVb.InvertFast();
519  if (!vVb.IsValid()) {
520  AliDebug(2,"Singular Matrix\n");
521  return dist;
522  }
523  dist = vVb(0,0)*dx*dx + vVb(1,1)*dy*dy + vVb(2,2)*dz*dz
524  + 2*vVb(0,1)*dx*dy + 2*vVb(0,2)*dx*dz + 2*vVb(1,2)*dy*dz;
525  return dist>0 ? TMath::Sqrt(dist) : -1;
526 }
527 
528 
529 
530 
531 
532 
533 
534 
535 
536 
537 
538 
539 
540 
541 
542 
543 
544 
545 
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)