AliPhysics  5eaf189 (5eaf189)
AliAnalysisTaskCheckEvSel.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2018, 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 /* $Id: $ */
17 
18 //*************************************************************************
19 // Class AliAnalysisTaskCheckEvSel
20 // Task to check event selection with D2H code
21 // Author: Francesco Prino
22 //*************************************************************************
23 
24 
25 #include <TList.h>
26 #include <TH1F.h>
27 #include <TH2F.h>
28 
29 #include "AliAnalysisManager.h"
30 #include "AliInputEventHandler.h"
31 #include "AliPIDResponse.h"
32 #include "AliAODHandler.h"
33 #include "AliAODEvent.h"
34 #include "AliVVertex.h"
35 #include "AliAODTrack.h"
36 #include "AliVertexingHFUtils.h"
38 #include "AliRDHFCutsD0toKpi.h"
40 
42 
43 
44 //________________________________________________________________________
47  fOutput(0x0),
48  fReadMC(kFALSE),
49  fSystem(0),
50  fCutOnzVertexSPD(0),
51  fAODProtection(1),
52  fHistNEvents(0x0),
53  fHistNEventsVsCent(0x0),
54  fHistNEventsVsCL1(0x0),
55  fHistNEventsVsWhyRej(0x0),
56  fHistNTrackletsBeforePileup(0x0),
57  fHistNTrackletsAfterPileup(0x0),
58  fHistNCL1BeforePileup(0x0),
59  fHistNCL1AfterPileup(0x0),
60  fHistCentrality(0x0),
61  fHistNTracksTPCoutVsV0Cent(0x0),
62  fHistNTracksFB4VsV0Cent(0x0),
63  fHistNTracksBC0VsV0Cent(0x0),
64  fHistNTrackletsVsV0Cent(0x0),
65  fHistNTracksTPCoutVsNTracklets(0x0),
66  fHistNTracksFB4VsNTracklets(0x0),
67  fHistNTracksBC0VsNTracksFB4(0x0),
68  fHistZVertexSPDBeforeCuts(0x0),
69  fHistZVertexSPDBeforeSPDCut(0x0),
70  fHistZVertexSPDAfterCuts(0x0),
71  fHistZVertexSPDBadTrackVert(0x0),
72  fCounter(0),
73  fAnalysisCuts(0x0)
74 {
75  // default constructor
76  fAnalysisCuts=new AliRDHFCutsD0toKpi("EvSelCuts");
77  fAnalysisCuts->SetTriggerMask(AliVEvent::kMB | AliVEvent::kINT7);
78  fAnalysisCuts->SetTriggerClass("");
79  fAnalysisCuts->SetOptPileup(AliRDHFCuts::kNoPileupSelection);
80 }
81 //________________________________________________________________________
83  AliAnalysisTaskSE("EvSelTask"),
84  fOutput(0x0),
85  fReadMC(readMC),
86  fSystem(system),
87  fCutOnzVertexSPD(0),
88  fAODProtection(1),
89  fHistNEvents(0x0),
90  fHistNEventsVsCent(0x0),
91  fHistNEventsVsCL1(0x0),
92  fHistNEventsVsWhyRej(0x0),
93  fHistNTrackletsBeforePileup(0x0),
94  fHistNTrackletsAfterPileup(0x0),
95  fHistNCL1BeforePileup(0x0),
96  fHistNCL1AfterPileup(0x0),
97  fHistCentrality(0x0),
98  fHistNTracksTPCoutVsV0Cent(0x0),
99  fHistNTracksFB4VsV0Cent(0x0),
100  fHistNTracksBC0VsV0Cent(0x0),
101  fHistNTrackletsVsV0Cent(0x0),
102  fHistNTracksTPCoutVsNTracklets(0x0),
103  fHistNTracksFB4VsNTracklets(0x0),
104  fHistNTracksBC0VsNTracksFB4(0x0),
105  fHistZVertexSPDBeforeCuts(0x0),
106  fHistZVertexSPDBeforeSPDCut(0x0),
107  fHistZVertexSPDAfterCuts(0x0),
108  fHistZVertexSPDBadTrackVert(0x0),
109  fCounter(0),
110  fAnalysisCuts(cuts)
111 {
112  // default constructor
113 
114  DefineOutput(1,TList::Class()); //My private output
115  DefineOutput(2,AliNormalizationCounter::Class());
116 }
117 
118 //________________________________________________________________________
120 {
121  //
122  // Destructor
123  //
124  if(fOutput && !fOutput->IsOwner()){
125  delete fHistNEvents;
126  delete fHistNEventsVsCent;
127  delete fHistNEventsVsCL1;
128  delete fHistNEventsVsWhyRej;
131  delete fHistNCL1BeforePileup;
132  delete fHistNCL1AfterPileup;
133  delete fHistCentrality;
143  delete fHistZVertexSPDAfterCuts;
145  }
146  delete fOutput;
147  delete fCounter;
148 }
149 
150 //________________________________________________________________________
152 {
153  // Create the output container
154  //
155  if(fDebug > 1) printf("AnalysisTaskCheckEvSel::UserCreateOutputObjects() \n");
156 
157  Double_t maxMult=200.;
158  if(fSystem==1) maxMult=5000.;
159  else if(fSystem==2) maxMult=500.;
160 
161  fOutput = new TList();
162  fOutput->SetOwner();
163  fOutput->SetName("OutputHistos");
164 
165  fHistNEvents = new TH1F("hNEvents", "number of events ",18,-0.5,17.5);
166  ConfigureEvSelAxis(fHistNEvents->GetXaxis());
167  fHistNEvents->SetMinimum(0);
168  fOutput->Add(fHistNEvents);
169 
170  fHistNEventsVsCent = new TH2F("hNEventsVsCent", " ; ; Centrality ",18,-0.5,17.5,101,0.,101.);
173 
174  fHistNEventsVsCL1 = new TH2F("hNEventsVsCL1", " ; ; N_{CL1}",18,-0.5,17.5,200,-0.5,2*maxMult-0.5);
177 
178  fHistNEventsVsWhyRej = new TH2F("hNEventsVsWhyRej", " ; ; WhyRej ",18,-0.5,17.5,11,-0.5,10.5);
181 
182  fHistNTrackletsBeforePileup = new TH1F("hNTrackletsBeforePileup"," ; N_{tracklets}",200,-0.5,maxMult-0.5);
183  fHistNTrackletsAfterPileup = new TH1F("hNTrackletsAfterPileup"," ; N_{tracklets}",200,-0.5,maxMult-0.5);
184  fHistNCL1BeforePileup = new TH1F("hNCL1BeforePileup"," ; N_{CL1}",200,-0.5,maxMult-0.5);
185  fHistNCL1AfterPileup = new TH1F("hNCL1AfterPileup"," ; N_{CL1}",200,-0.5,maxMult-0.5);
190 
191  fHistCentrality = new TH1F("hCentrality"," ; Centrality ; ",105,0.,105.);
192  fOutput->Add(fHistCentrality);
193 
194  fHistNTracksTPCoutVsV0Cent = new TH2F("hNTracksTPCoutVsV0Cent"," ; Centrality ; N_{tracks, TPCout}",105,0.,105.,200,-0.5,2*maxMult-0.5);
195  fHistNTracksFB4VsV0Cent = new TH2F("hNTracksFB4VsV0Cent"," ; Centrality ; N_{tracks, FiltBit4}",105,0.,105.,200,-0.5,maxMult-0.5);
196  fHistNTracksBC0VsV0Cent = new TH2F("hNTracksBC0VsV0Cent"," ; Centrality ; N_{tracks, TOFBC=0}",105,0.,105.,200,-0.5,maxMult-0.5);
197  fHistNTrackletsVsV0Cent = new TH2F("hNTrackletsVsV0Cent"," ; Centrality ; N_{tracklets}",105,0.,105.,200,-0.5,maxMult-0.5);
198  fHistNTracksTPCoutVsNTracklets = new TH2F("hNTracksTPCoutVsNTracklets"," ; N_{tracklets} ; N_{tracks, TPCout}",200,-0.5,maxMult-0.5,200,-0.5,2*maxMult-0.5);
199  fHistNTracksFB4VsNTracklets = new TH2F("hNTracksFB4VsNTracklets"," ; N_{tracklets} ; N_{tracks, FiltBit4}",200,-0.5,maxMult-0.5,200,-0.5,maxMult-0.5);
200  fHistNTracksBC0VsNTracksFB4 = new TH2F("hNTracksBC0VsNTracksFB4"," ; N_{tracks, FiltBit4}; N_{tracks, TOFBC=0}",200,-0.5,maxMult-0.5,200,-0.5,maxMult-0.5);
208 
209  fHistZVertexSPDBeforeCuts = new TH2F("hZVertexSPDBeforeCuts"," ; z_{SPDvertex} ; z_{TRKvertex}",400,-20.,20.,400,-20.,20.);
211  fHistZVertexSPDBeforeSPDCut = new TH2F("hZVertexSPDBeforeSPDCut"," ; z_{SPDvertex} ; z_{TRKvertex}",400,-20.,20.,400,-20.,20.);
213  fHistZVertexSPDAfterCuts = new TH2F("hZVertexSPDAfterCuts"," ; z_{SPDvertex} ; z_{TRKvertex}",400,-20.,20.,400,-20.,20.);
215  fHistZVertexSPDBadTrackVert =new TH2F("hZVertexSPDBadTrackVert"," ; z_{SPDvertex} ; z_{TRKvertex}",400,-20.,20.,400,-20.,20.);
217 
218  TString normName="NormalizationCounter";
219  AliAnalysisDataContainer *cont = GetOutputSlot(2)->GetContainer();
220  if(cont)normName=(TString)cont->GetName();
221  fCounter = new AliNormalizationCounter(normName.Data());
222  fCounter->Init();
223 
224  PostData(1,fOutput);
225 }
226 
227 //________________________________________________________________________
231  ax->SetBinLabel(1,"nEvents read");
232  ax->SetBinLabel(2,"Rejected due to mismatch in trees");
233  ax->SetBinLabel(3,"nEvents with good AOD");
234  ax->SetBinLabel(4,"Rejected due to trigger");
235  ax->SetBinLabel(5,"Rejected due to phys sel");
236  ax->SetBinLabel(6,"Rejected due to bad centrality estimator");
237  ax->SetBinLabel(7,"Rejected due to centrality flattening");
238  ax->SetBinLabel(8,"Rejected due to centrality out of range");
239  ax->SetBinLabel(9,"Rejected due to not reco vertex");
240  ax->SetBinLabel(10,"Rejected for contr vertex");
241  ax->SetBinLabel(11,"Rejected for bad track vertex");
242  ax->SetBinLabel(12,"Rejected for vertex out of accept");
243  ax->SetBinLabel(13,"Rejected for pileup events");
244  ax->SetBinLabel(14,"Passing Phys sel + trigger");
245  ax->SetBinLabel(15,"Passing Phys sel + trigger + Pileup");
246  ax->SetBinLabel(16,"Passing Phys sel + trigger + Pileup + zVertex");
247  ax->SetBinLabel(17,"Passing Phys sel + trigger + Pileup + zVertex + Centrality");
248  ax->SetBinLabel(18,"Passing IsEventSelected");
249  ax->SetNdivisions(1,kFALSE);
250 }
251 //________________________________________________________________________
253  //Build the 3-track combinatorics (+-+ and -+-) for D+->Kpipi decays
254 
255  if(fDebug > 1) printf("AnalysisTaskCheckEvSel::UserExec() \n");
256  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
257 
258  fHistNEvents->Fill(0); // count event
259 
260  if(fAODProtection>=0){
261  // Protection against different number of events in the AOD and deltaAOD
262  // In case of discrepancy the event is rejected.
263  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
264  if (matchingAODdeltaAODlevel<0 || (matchingAODdeltaAODlevel==0 && fAODProtection==1)) {
265  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
266  fHistNEvents->Fill(1);
267  return;
268  }
269  }
270 
271  if(!aod && AODEvent() && IsStandardAOD()) {
272  // In case there is an AOD handler writing a standard AOD, use the AOD
273  // event in memory rather than the input (ESD) event.
274  aod = dynamic_cast<AliAODEvent*> (AODEvent());
275  }
276  if(!aod){
277  printf("AliAnalysisTaskCheckEvSel::UserExec: AOD not found!\n");
278  return;
279  }
280 
281  // fix for temporary bug in ESDfilter
282  // the AODs with null vertex pointer didn't pass the PhysSel
283  if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
284 
285 
286  // Post the data already here
287 
288  Bool_t isEvSel=fAnalysisCuts->IsEventSelected(aod);
291  const AliVVertex *vertex = aod->GetPrimaryVertex();
292  const AliVVertex *vertexSPD = aod->GetPrimaryVertexSPD();
294  Double_t ncl1 = aod->GetNumberOfITSClusters(1);
295  Double_t zvSPD=vertexSPD->GetZ();
296  Double_t zvTRK=vertex->GetZ();
297 
307 
308  fHistZVertexSPDBeforeCuts->Fill(zvSPD,zvTRK);
313  vertexSPD->GetNContributors()>=1){
314  fHistZVertexSPDBeforeSPDCut->Fill(zvSPD,zvTRK);
315  Double_t dz = vertexSPD->GetZ()-vertex->GetZ();
316  Bool_t okSpdTrk=kTRUE;
317  if(fCutOnzVertexSPD==2 && TMath::Abs(dz)>0.5) okSpdTrk=kFALSE;
318  if(okSpdTrk && fCutOnzVertexSPD==3){
319  double covTrc[6],covSPD[6];
320  vertex->GetCovarianceMatrix(covTrc);
321  vertexSPD->GetCovarianceMatrix(covSPD);
322  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
323  double errTrc = TMath::Sqrt(covTrc[5]);
324  double nsigTot = TMath::Abs(dz)/errTot, nsigTrc = TMath::Abs(dz)/errTrc;
325  if (TMath::Abs(dz)>0.2 || nsigTot>10 || nsigTrc>20) okSpdTrk=kFALSE;
326  }
327  if(!okSpdTrk) fHistZVertexSPDBadTrackVert->Fill(zvSPD,zvTRK);
328  }
329  if(isEvSel) fHistZVertexSPDAfterCuts->Fill(zvSPD,zvTRK);
330 
331  fHistNEvents->Fill(2);
332  fHistNEventsVsCent->Fill(2,centr);
333  fHistNEventsVsCL1->Fill(2,ncl1);
334 
335  Int_t binToFill=-1;
337  if(fAnalysisCuts->IsEventRejectedDueToTrigger()) binToFill=3;
339  }else{
341  if(wrej==3) binToFill=5;
342  else if(wrej==4) binToFill=6;
343  else binToFill=7;
344  }else{
349  }else{
351  else if(fAnalysisCuts->IsEventRejectedDueToPileup()) binToFill=12;
352  }
353  }
354  }
355  if(binToFill>0){
356  fHistNEvents->Fill(binToFill);
357  fHistNEventsVsCent->Fill(binToFill,centr);
358  fHistNEventsVsCL1->Fill(binToFill,ncl1);
359  }
361  fHistNEvents->Fill(13);
362  fHistNEventsVsCent->Fill(13,centr);
363  fHistNEventsVsCL1->Fill(13,ncl1);
364  fHistNEventsVsWhyRej->Fill(13,wrej);
366  fHistNEvents->Fill(14);
367  fHistNEventsVsCent->Fill(14,centr);
368  fHistNEventsVsCL1->Fill(14,ncl1);
369  fHistNEventsVsWhyRej->Fill(14,wrej);
371  fHistNEvents->Fill(15);
372  fHistNEventsVsCent->Fill(15,centr);
373  fHistNEventsVsCL1->Fill(15,ncl1);
374  fHistNEventsVsWhyRej->Fill(15,wrej);
376  fHistNEvents->Fill(16);
377  fHistNEventsVsCent->Fill(16,centr);
378  fHistNEventsVsCL1->Fill(16,ncl1);
379  fHistNEventsVsWhyRej->Fill(16,wrej);
380  }
381  }
382  }
383  }
384  if(isEvSel || (!isEvSel && wrej==1)){
385  fHistNTrackletsBeforePileup->Fill(ntrkl);
386  fHistNCL1BeforePileup->Fill(ncl1);
387  }
388  if(isEvSel){
389  fHistNEvents->Fill(17);
390  fHistNEventsVsCent->Fill(17,centr);
391  fHistNEventsVsCL1->Fill(17,ncl1);
392  fHistNEventsVsWhyRej->Fill(17,wrej);
393  fHistNTrackletsAfterPileup->Fill(ntrkl);
394  fHistNCL1AfterPileup->Fill(ncl1);
395  }
396 
398  // events not passing the centrality selection can be removed immediately. For the others we must count the generated D mesons
400 
401  if(isEvSel){
402 
403  Int_t ntracksTPCout=0;
404  Int_t ntracksFB4=0;
405  Int_t ntracksBC0=0;
406 
407  Double_t magField = aod->GetMagneticField();
408 
409  for(Int_t iTr=0; iTr<aod->GetNumberOfTracks(); iTr++){
410  AliAODTrack* track=dynamic_cast<AliAODTrack*>(aod->GetTrack(iTr));
411 
412  if(track->GetStatus() & AliESDtrack::kTPCout) ntracksTPCout++;
413  if(track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)){
414  ntracksFB4++;
415  Int_t tofBC=track->GetTOFBunchCrossing(magField);
416  if(tofBC==0) ntracksBC0++;
417  }
418  }
419  fHistCentrality->Fill(centr);
420  fHistNTracksTPCoutVsV0Cent->Fill(centr,ntracksTPCout);
421  fHistNTracksFB4VsV0Cent->Fill(centr,ntracksFB4);
422  fHistNTracksBC0VsV0Cent->Fill(centr,ntracksBC0);
423  fHistNTrackletsVsV0Cent->Fill(centr,ntrkl);
424  fHistNTracksTPCoutVsNTracklets->Fill(ntracksTPCout,ntrkl);
425  fHistNTracksFB4VsNTracklets->Fill(ntracksFB4,ntrkl);
426  fHistNTracksBC0VsNTracksFB4->Fill(ntracksBC0,ntracksFB4);
427  }
428 
429  PostData(1,fOutput);
430  PostData(2,fCounter);
431 
432  return;
433 }
434 
435 //_________________________________________________________________
437 {
438  // Terminate analysis
439  //
440  if(fDebug > 1) printf("AliAnalysisTaskCheckEvSel: Terminate() \n");
441  fOutput = dynamic_cast<TList*> (GetOutputData(1));
442  if (!fOutput) {
443  printf("ERROR: fOutput not available\n");
444  return;
445  }
446  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
447  if(fHistNEvents){
448  printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
449  }else{
450  printf("ERROR: fHistNEvents not available\n");
451  return;
452  }
453  return;
454 }
455 
Bool_t IsEventRejectedDueToCentrality() const
Definition: AliRDHFCuts.h:330
Bool_t IsEventRejectedDueToZVertexOutsideFiducialRegion() const
Definition: AliRDHFCuts.h:324
Bool_t IsEventRejectedDueToNotRecoVertex() const
Definition: AliRDHFCuts.h:318
double Double_t
Definition: External.C:58
TH2F * fHistNTracksTPCoutVsV0Cent
! Centrality-multiplicity correl
Definition: External.C:236
TH1F * fHistCentrality
! hist. of centrality distribution
AliNormalizationCounter * fCounter
!Counter for normalization
Int_t IsEventSelectedInCentrality(AliVEvent *event)
TH2F * fHistNTracksBC0VsNTracksFB4
! Centrality-multiplicity correl
static Int_t CheckMatchingAODdeltaAODevents()
Bool_t IsEventRejectedDueToVertexContributors() const
Definition: AliRDHFCuts.h:321
Int_t GetWhyRejection() const
Definition: AliRDHFCuts.h:313
TH2F * fHistNTrackletsVsV0Cent
! Centrality-multiplicity correl
TH2F * fHistNEventsVsCent
! hist. for No. of events
virtual void Terminate(Option_t *option)
TList * fOutput
! list send on output slot 0
TH2F * fHistNTracksTPCoutVsNTracklets
! Centrality-multiplicity correl
static Int_t GetNumberOfTrackletsInEtaRange(AliAODEvent *ev, Double_t mineta, Double_t maxeta)
TH2F * fHistZVertexSPDBadTrackVert
! z-vertex distr.
TH1F * fHistNEvents
! hist. for No. of events
TH2F * fHistNEventsVsWhyRej
! hist. for No. of events
int Int_t
Definition: External.C:63
TH2F * fHistZVertexSPDAfterCuts
! z-vertex distr.
Bool_t IsEventRejectedDueToBadTrackVertex() const
Definition: AliRDHFCuts.h:339
TH2F * fHistNTracksFB4VsV0Cent
! Centrality-multiplicity correl
virtual void UserExec(Option_t *option)
Float_t GetCentrality(AliAODEvent *aodEvent)
Definition: AliRDHFCuts.h:257
Bool_t IsEventRejectedDueToPileup() const
Definition: AliRDHFCuts.h:327
TH1F * fHistNCL1BeforePileup
! hist. for No. of tracklets
TH2F * fHistNTracksBC0VsV0Cent
! Centrality-multiplicity correl
Bool_t IsEventRejectedDuePhysicsSelection() const
Definition: AliRDHFCuts.h:345
Int_t fAODProtection
cut on zSPD vertex to remove outliers in centrality vs. tracklets (0=no cut, 1= cut at 12 cm...
TH2F * fHistNTracksFB4VsNTracklets
! Centrality-multiplicity correl
Bool_t IsEventSelected(AliVEvent *event)
void StoreEvent(AliVEvent *, AliRDHFCuts *, Bool_t mc=kFALSE, Int_t multiplicity=-9999, Double_t spherocity=-99.)
TH1F * fHistNTrackletsBeforePileup
! hist. for No. of tracklets
TH2F * fHistZVertexSPDBeforeSPDCut
! z-vertex distr.
TH2F * fHistNEventsVsCL1
! hist. for No. of events
TH1F * fHistNTrackletsAfterPileup
! hist. for No. of tracklets
const char Option_t
Definition: External.C:48
Bool_t IsEventRejectedDueToTrigger() const
Definition: AliRDHFCuts.h:315
bool Bool_t
Definition: External.C:53
TH1F * fHistNCL1AfterPileup
! hist. for No. of tracklets
Int_t GetUseCentrality() const
Definition: AliRDHFCuts.h:279
Int_t fCutOnzVertexSPD
0=pp, 1=Pb-Pb, 2=p-Pb
Bool_t IsEventRejectedDueToCentralityFlattening() const
Definition: AliRDHFCuts.h:333
Double_t maxMult
TH2F * fHistZVertexSPDBeforeCuts
! z-vertex distr.