]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWG/FLOW/Tasks/AliFlowEventCuts.cxx
Moving/split PWG2/FLOW to PWGCF/FLOW, PWG/FLOW/Base, PWG/FLOW/Tasks, PWG/Glauber
[u/mrichter/AliRoot.git] / PWG / FLOW / Tasks / AliFlowEventCuts.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, 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 // AliFlowEventCuts:
19 // An event cut class for the flow framework
20 //
21 // origin: Mikolaj Krzewicki (mikolaj.krzewicki@cern.ch)
22
23 #include <limits.h>
24 #include <float.h>
25 #include <TList.h>
26 #include <TH1F.h>
27 #include <TH2F.h>
28 #include <TBrowser.h>
29 #include "TMath.h"
30 #include "TNamed.h"
31 #include "AliVVertex.h"
32 #include "AliVEvent.h"
33 #include "AliESDEvent.h"
34 #include "AliCentrality.h"
35 #include "AliESDVZERO.h"
36 #include "AliMultiplicity.h"
37 #include "AliMCEvent.h"
38 #include "AliFlowEventCuts.h"
39 #include "AliFlowTrackCuts.h"
40 #include "AliTriggerAnalysis.h"
41
42 ClassImp(AliFlowEventCuts)
43
44 //-----------------------------------------------------------------------
45 AliFlowEventCuts::AliFlowEventCuts():
46   TNamed(),
47   fQA(NULL),
48   fCutNumberOfTracks(kFALSE),
49   fNumberOfTracksMax(INT_MAX),
50   fNumberOfTracksMin(INT_MIN),
51   fCutRefMult(kFALSE),
52   fRefMultMethod(kTPConly),
53   fUseAliESDtrackCutsRefMult(kFALSE),
54   fRefMultMethodAliESDtrackCuts(AliESDtrackCuts::kTrackletsITSTPC),
55   fRefMultMax(INT_MAX),
56   fRefMultMin(INT_MIN),
57   fRefMultCuts(NULL),
58   fMeanPtCuts(NULL),
59   fStandardTPCcuts(NULL),
60   fStandardGlobalCuts(NULL),
61   fCutPrimaryVertexX(kFALSE),
62   fPrimaryVertexXmax(INT_MAX),
63   fPrimaryVertexXmin(INT_MIN),
64   fCutPrimaryVertexY(kFALSE),
65   fPrimaryVertexYmax(INT_MAX),
66   fPrimaryVertexYmin(INT_MIN),
67   fCutPrimaryVertexZ(kFALSE),
68   fPrimaryVertexZmax(INT_MAX),
69   fPrimaryVertexZmin(INT_MIN),
70   fCutNContributors(kFALSE),
71   fNContributorsMax(INT_MAX),
72   fNContributorsMin(INT_MIN),
73   fCutMeanPt(kFALSE),
74   fMeanPtMax(-DBL_MAX),
75   fMeanPtMin(DBL_MAX),
76   fCutSPDvertexerAnomaly(kFALSE),
77   fCutTPCmultiplicityOutliers(kFALSE),
78   fCutCentralityPercentile(kFALSE),
79   fUseCentralityUnchecked(kFALSE),
80   fCentralityPercentileMethod(kTPConly),
81   fCentralityPercentileMax(100.),
82   fCentralityPercentileMin(0.),
83   fCutZDCtiming(kFALSE),
84   fTrigAna()
85 {
86   //constructor 
87 }
88
89 //-----------------------------------------------------------------------
90 AliFlowEventCuts::AliFlowEventCuts(const char* name, const char* title):
91   TNamed(name, title),
92   fQA(NULL),
93   fCutNumberOfTracks(kFALSE),
94   fNumberOfTracksMax(INT_MAX),
95   fNumberOfTracksMin(INT_MIN),
96   fCutRefMult(kFALSE),
97   fRefMultMethod(kTPConly),
98   fUseAliESDtrackCutsRefMult(kFALSE),
99   fRefMultMethodAliESDtrackCuts(AliESDtrackCuts::kTrackletsITSTPC),
100   fRefMultMax(INT_MAX),
101   fRefMultMin(INT_MIN),
102   fRefMultCuts(NULL),
103   fMeanPtCuts(NULL),
104   fStandardTPCcuts(AliFlowTrackCuts::GetStandardTPCStandaloneTrackCuts2010()),
105   fStandardGlobalCuts(AliFlowTrackCuts::GetStandardGlobalTrackCuts2010()),
106   fCutPrimaryVertexX(kFALSE),
107   fPrimaryVertexXmax(INT_MAX),
108   fPrimaryVertexXmin(INT_MIN),
109   fCutPrimaryVertexY(kFALSE),
110   fPrimaryVertexYmax(INT_MAX),
111   fPrimaryVertexYmin(INT_MIN),
112   fCutPrimaryVertexZ(kFALSE),
113   fPrimaryVertexZmax(INT_MAX),
114   fPrimaryVertexZmin(INT_MIN),
115   fCutNContributors(kFALSE),
116   fNContributorsMax(INT_MAX),
117   fNContributorsMin(INT_MIN),
118   fCutMeanPt(kFALSE),
119   fMeanPtMax(-DBL_MAX),
120   fMeanPtMin(DBL_MAX),
121   fCutSPDvertexerAnomaly(kFALSE),
122   fCutTPCmultiplicityOutliers(kFALSE),
123   fCutCentralityPercentile(kFALSE),
124   fUseCentralityUnchecked(kFALSE),
125   fCentralityPercentileMethod(kTPConly),
126   fCentralityPercentileMax(100.),
127   fCentralityPercentileMin(0.),
128   fCutZDCtiming(kFALSE),
129   fTrigAna()
130 {
131   //constructor 
132 }
133
134 ////-----------------------------------------------------------------------
135 AliFlowEventCuts::AliFlowEventCuts(const AliFlowEventCuts& that):
136   TNamed(that),
137   fQA(NULL),
138   fCutNumberOfTracks(that.fCutNumberOfTracks),
139   fNumberOfTracksMax(that.fNumberOfTracksMax),
140   fNumberOfTracksMin(that.fNumberOfTracksMin),
141   fCutRefMult(that.fCutRefMult),
142   fRefMultMethod(that.fRefMultMethod),
143   fUseAliESDtrackCutsRefMult(that.fUseAliESDtrackCutsRefMult),
144   fRefMultMethodAliESDtrackCuts(that.fRefMultMethodAliESDtrackCuts),
145   fRefMultMax(that.fRefMultMax),
146   fRefMultMin(that.fRefMultMin),
147   fRefMultCuts(NULL),
148   fMeanPtCuts(NULL),
149   fStandardTPCcuts(NULL),
150   fStandardGlobalCuts(NULL),
151   fCutPrimaryVertexX(that.fCutPrimaryVertexX),
152   fPrimaryVertexXmax(that.fPrimaryVertexXmax),
153   fPrimaryVertexXmin(that.fPrimaryVertexXmin),
154   fCutPrimaryVertexY(that.fCutPrimaryVertexX),
155   fPrimaryVertexYmax(that.fPrimaryVertexYmax),
156   fPrimaryVertexYmin(that.fPrimaryVertexYmin),
157   fCutPrimaryVertexZ(that.fCutPrimaryVertexX),
158   fPrimaryVertexZmax(that.fPrimaryVertexZmax),
159   fPrimaryVertexZmin(that.fPrimaryVertexZmin),
160   fCutNContributors(that.fCutNContributors),
161   fNContributorsMax(that.fNContributorsMax),
162   fNContributorsMin(that.fNContributorsMin),
163   fCutMeanPt(that.fCutMeanPt),
164   fMeanPtMax(that.fMeanPtMax),
165   fMeanPtMin(that.fMeanPtMin),
166   fCutSPDvertexerAnomaly(that.fCutSPDvertexerAnomaly),
167   fCutTPCmultiplicityOutliers(that.fCutTPCmultiplicityOutliers),
168   fCutCentralityPercentile(that.fCutCentralityPercentile),
169   fUseCentralityUnchecked(that.fUseCentralityUnchecked),
170   fCentralityPercentileMethod(that.fCentralityPercentileMethod),
171   fCentralityPercentileMax(that.fCentralityPercentileMax),
172   fCentralityPercentileMin(that.fCentralityPercentileMin),
173   fCutZDCtiming(that.fCutZDCtiming),
174   fTrigAna()
175 {
176   if (that.fQA) DefineHistograms();
177   //copy constructor 
178   if (that.fRefMultCuts)
179     fRefMultCuts = new AliFlowTrackCuts(*(that.fRefMultCuts));
180   if (that.fMeanPtCuts)
181     fMeanPtCuts = new AliFlowTrackCuts(*(that.fMeanPtCuts));
182   fStandardTPCcuts = AliFlowTrackCuts::GetStandardTPCStandaloneTrackCuts2010();
183   fStandardGlobalCuts = AliFlowTrackCuts::GetStandardGlobalTrackCuts2010();
184 }
185
186 ////-----------------------------------------------------------------------
187 AliFlowEventCuts::~AliFlowEventCuts()
188 {
189   //dtor
190   delete fMeanPtCuts;
191   delete fRefMultCuts;
192   delete fStandardGlobalCuts;
193   delete fStandardTPCcuts;
194   if (fQA) { fQA->SetOwner(); fQA->Delete(); delete fQA; }
195 }
196
197 ////-----------------------------------------------------------------------
198 AliFlowEventCuts& AliFlowEventCuts::operator=(const AliFlowEventCuts& that)
199 {
200   //assignment
201   if (this==&that) return *this;
202
203   if (that.fQA)
204   {
205     if (fQA)
206     {
207       fQA->Delete();
208       delete fQA;
209     }
210     fQA = static_cast<TList*>(that.fQA->Clone());
211   }
212   else
213   {
214     fQA->Delete();
215     delete fQA;
216     fQA=NULL;
217   }
218
219   fCutNumberOfTracks=that.fCutNumberOfTracks;
220   fNumberOfTracksMax=that.fNumberOfTracksMax;
221   fNumberOfTracksMin=that.fNumberOfTracksMin;
222   fCutRefMult=that.fCutRefMult;
223   fRefMultMethod=that.fRefMultMethod;
224   fUseAliESDtrackCutsRefMult=that.fUseAliESDtrackCutsRefMult;
225   fRefMultMethodAliESDtrackCuts=that.fRefMultMethodAliESDtrackCuts;
226   fRefMultMax=that.fRefMultMax;
227   fRefMultMin=that.fRefMultMin;
228   if (that.fRefMultCuts) *fRefMultCuts=*(that.fRefMultCuts);
229   if (that.fMeanPtCuts) *fMeanPtCuts=*(that.fMeanPtCuts);
230   fStandardTPCcuts = AliFlowTrackCuts::GetStandardTPCStandaloneTrackCuts2010();
231   fStandardGlobalCuts = AliFlowTrackCuts::GetStandardGlobalTrackCuts2010();
232   fCutPrimaryVertexX=that.fCutPrimaryVertexX;
233   fPrimaryVertexXmax=that.fPrimaryVertexXmax;
234   fPrimaryVertexXmin=that.fPrimaryVertexXmin;
235   fCutPrimaryVertexY=that.fCutPrimaryVertexY;
236   fPrimaryVertexYmax=that.fPrimaryVertexYmax;
237   fPrimaryVertexYmin=that.fPrimaryVertexYmin;
238   fCutPrimaryVertexZ=that.fCutPrimaryVertexZ;
239   fPrimaryVertexZmax=that.fPrimaryVertexZmax;
240   fPrimaryVertexZmin=that.fPrimaryVertexZmin;
241   fCutNContributors=that.fCutNContributors;
242   fNContributorsMax=that.fNContributorsMax;
243   fNContributorsMin=that.fNContributorsMin;
244   fCutMeanPt=that.fCutMeanPt;
245   fMeanPtMax=that.fMeanPtMax;
246   fMeanPtMin=that.fMeanPtMin;
247   fCutSPDvertexerAnomaly=that.fCutSPDvertexerAnomaly;
248   fCutTPCmultiplicityOutliers=that.fCutTPCmultiplicityOutliers;
249   fCutCentralityPercentile=that.fCutCentralityPercentile;
250   fUseCentralityUnchecked=that.fUseCentralityUnchecked;
251   fCentralityPercentileMethod=that.fCentralityPercentileMethod;
252   fCentralityPercentileMax=that.fCentralityPercentileMax;
253   fCentralityPercentileMin=that.fCentralityPercentileMin;
254   fCutZDCtiming=that.fCutZDCtiming;
255   return *this;
256 }
257
258 //----------------------------------------------------------------------- 
259 Bool_t AliFlowEventCuts::IsSelected(TObject* obj)
260 {
261   //check cuts
262   AliVEvent* vevent = dynamic_cast<AliVEvent*>(obj);
263   if (vevent) return PassesCuts(vevent);
264   return kFALSE;  //when passed wrong type of object
265 }
266 //----------------------------------------------------------------------- 
267 Bool_t AliFlowEventCuts::PassesCuts(AliVEvent *event)
268 {
269   ///check if event passes cuts
270   const AliVVertex* pvtx=event->GetPrimaryVertex();
271   Double_t pvtxx = pvtx->GetX();
272   Double_t pvtxy = pvtx->GetY();
273   Double_t pvtxz = pvtx->GetZ();
274   Int_t ncontrib = pvtx->GetNContributors();
275   Bool_t pass=kTRUE;
276   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*>(event);
277   Int_t multTPC = 0;
278   Int_t multGlobal = 0; 
279   if (fQA)
280   {
281     multTPC = fStandardTPCcuts->Count(event);
282     multGlobal = fStandardGlobalCuts->Count(event);
283     QAbefore(0)->Fill(pvtxz);
284     QAbefore(1)->Fill(multGlobal,multTPC);
285   }
286   if (fCutTPCmultiplicityOutliers)
287   {
288     //this is pretty slow as we check the event track by track twice
289     //this cut will work for 2010 PbPb data and is dependent on
290     //TPC and ITS reco efficiency (e.g. geometry, calibration etc)
291     if (!fQA)
292     {
293       multTPC = fStandardTPCcuts->Count(event);
294       multGlobal = fStandardGlobalCuts->Count(event);
295     }
296     if (multTPC > ( 23+1.216*multGlobal)) {pass=kFALSE;}
297     if (multTPC < (-20+1.087*multGlobal)) {pass=kFALSE;}
298   }
299   if (fCutNContributors)
300   {
301     if (ncontrib < fNContributorsMin || ncontrib >= fNContributorsMax) pass=kFALSE;
302   }
303   if (fCutPrimaryVertexX)
304   {
305     if (pvtxx < fPrimaryVertexXmin || pvtxx >= fPrimaryVertexXmax) pass=kFALSE;
306   }
307   if (fCutPrimaryVertexY)
308   {
309     if (pvtxy < fPrimaryVertexYmin || pvtxy >= fPrimaryVertexYmax) pass=kFALSE;
310   }
311   if (fCutPrimaryVertexZ)
312   {
313     if (pvtxz < fPrimaryVertexZmin || pvtxz >= fPrimaryVertexZmax)
314       pass=kFALSE;
315   }
316   if (fCutCentralityPercentile&&esdevent)
317   {
318     AliCentrality* centr = esdevent->GetCentrality();
319     if (fUseCentralityUnchecked)
320     {
321       if (!centr->IsEventInCentralityClassUnchecked( fCentralityPercentileMin,
322                                                      fCentralityPercentileMax,
323                                                      CentrMethName(fCentralityPercentileMethod) ))
324       {
325         pass=kFALSE;
326       }
327     }
328     else
329     {
330       if (!centr->IsEventInCentralityClass( fCentralityPercentileMin,
331                                             fCentralityPercentileMax,
332                                             CentrMethName(fCentralityPercentileMethod) ))
333       {
334         pass=kFALSE;
335       }
336     }
337   }
338   if (fCutSPDvertexerAnomaly&&esdevent)
339   {
340     const AliESDVertex* sdpvertex = esdevent->GetPrimaryVertexSPD();
341     if (sdpvertex->GetNContributors()<1) pass=kFALSE;
342     if (sdpvertex->GetDispersion()>0.04) pass=kFALSE;
343     if (sdpvertex->GetZRes()>0.25) pass=kFALSE;
344     const AliESDVertex* tpcvertex = esdevent->GetPrimaryVertexTPC();
345     if (tpcvertex->GetNContributors()<1) pass=kFALSE;
346     const AliMultiplicity* tracklets = esdevent->GetMultiplicity();
347     if (tpcvertex->GetNContributors()<(-10.0+0.25*tracklets->GetNumberOfITSClusters(0)))
348     {
349       pass=kFALSE;
350     }
351   }
352   if (fCutZDCtiming&&esdevent)
353   {
354     if (!fTrigAna.ZDCTimeTrigger(esdevent))
355     {
356       pass=kFALSE;
357     }
358   }
359   if(fCutNumberOfTracks) {if ( event->GetNumberOfTracks() < fNumberOfTracksMin ||
360                                event->GetNumberOfTracks() >= fNumberOfTracksMax ) pass=kFALSE;}
361   if(fCutRefMult&&esdevent)
362   {
363     //reference multiplicity still to be defined
364     Double_t refMult = RefMult(event);
365     if (refMult < fRefMultMin || refMult >= fRefMultMax )
366     {
367       pass=kFALSE;
368     }
369   }
370   if (fCutMeanPt)
371   {
372     Float_t meanpt=0.0;
373     Int_t ntracks=event->GetNumberOfTracks();
374     Int_t nselected=0;
375     for (Int_t i=0; i<ntracks; i++)
376     {
377       AliVParticle* track = event->GetTrack(i);
378       if (!track) continue;
379       Bool_t localpass=kTRUE;
380       if (fMeanPtCuts) localpass=fMeanPtCuts->IsSelected(track);
381       if (localpass) 
382       {
383         meanpt += track->Pt();
384         nselected++;
385       }
386     }
387     if (nselected) meanpt=meanpt/nselected;
388     if (meanpt<fMeanPtMin || meanpt >= fMeanPtMax) pass=kFALSE;
389   }
390   if (fQA&&pass) 
391   {
392     QAafter(1)->Fill(multGlobal,multTPC);
393     QAafter(0)->Fill(pvtxz);
394   }
395   return pass;
396 }
397
398 //----------------------------------------------------------------------- 
399 const char* AliFlowEventCuts::CentrMethName(refMultMethod method) const
400 {
401   //get the string for refmultmethod, for use with AliCentrality in
402   //the cut on centrality percentile
403   switch (method)
404   {
405     case kSPDtracklets:
406       return "TKL";
407     case kSPD1clusters:
408       return "CL1";
409     case kTPConly:
410       return "TRK";
411     case kV0:
412       return "V0M";
413     default:
414       return "";
415   }
416 }
417 //----------------------------------------------------------------------- 
418 AliFlowEventCuts* AliFlowEventCuts::StandardCuts()
419 {
420   //make a set of standard event cuts, caller becomes owner
421   AliFlowEventCuts* cuts = new AliFlowEventCuts();
422   return cuts;
423 }
424
425 //----------------------------------------------------------------------- 
426 Int_t AliFlowEventCuts::RefMult(AliVEvent* event)
427 {
428   //calculate the reference multiplicity, if all fails return 0
429   AliESDVZERO* vzero = NULL;
430   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*>(event);
431
432   if (fUseAliESDtrackCutsRefMult && esdevent)
433   {
434     //use the standard ALICE reference multiplicity with the default eta range
435     return AliESDtrackCuts::GetReferenceMultiplicity(esdevent, fRefMultMethodAliESDtrackCuts);
436   }
437
438   if (fRefMultMethod==kTPConly && !fRefMultCuts)
439   {
440     fRefMultCuts = AliFlowTrackCuts::GetStandardTPCStandaloneTrackCuts();
441     fRefMultCuts->SetEtaRange(-0.8,0.8);
442     fRefMultCuts->SetPtMin(0.15);
443   }
444   else if (fRefMultMethod==kSPDtracklets && !fRefMultCuts)
445   {
446     fRefMultCuts = new AliFlowTrackCuts("tracklet refmult cuts");
447     fRefMultCuts->SetParamType(AliFlowTrackCuts::kSPDtracklet);
448     fRefMultCuts->SetEtaRange(-0.8,0.8);
449   }
450   else if (fRefMultMethod==kV0)
451   {
452     if (!esdevent) return 0;
453     vzero=esdevent->GetVZEROData();
454     if (!vzero) return 0;
455     return TMath::Nint(vzero->GetMTotV0A()+vzero->GetMTotV0C());
456   }
457   else if (fRefMultMethod==kSPD1clusters)
458   {
459     if (!esdevent) return 0;
460     const AliMultiplicity* mult = esdevent->GetMultiplicity();
461     if (!mult) return 0;
462     return mult->GetNumberOfITSClusters(1);
463   }
464
465   Int_t refmult=0;
466   fRefMultCuts->SetEvent(event);
467   for (Int_t i=0; i<fRefMultCuts->GetNumberOfInputObjects(); i++)
468   {
469     if (fRefMultCuts->IsSelected(fRefMultCuts->GetInputObject(i),i))
470       refmult++;
471   }
472   return refmult;
473 }
474 //_____________________________________________________________________________
475 void AliFlowEventCuts::DefineHistograms()
476 {
477   //define QA histos
478   if (fQA) return;
479
480   Bool_t adddirstatus = TH1::AddDirectoryStatus();
481   TH1::AddDirectory(kFALSE);
482   fQA = new TList(); fQA->SetOwner();
483   fQA->SetName(Form("%s QA",GetName()));
484   TList* before = new TList(); before->SetOwner();
485   before->SetName("before");
486   TList* after = new TList(); after->SetOwner();
487   after->SetName("after");
488   fQA->Add(before);
489   fQA->Add(after);
490   before->Add(new TH1F("zvertex",";z;event cout",500,-15.,15.)); //0
491   after->Add(new TH1F("zvertex",";z;event cout",500,-15.,15.)); //0
492   before->Add(new TH2F("fTPCvsGlobalMult","TPC only vs Global track multiplicity;global;TPC only",500,0,2500,500,0,3500));//1
493   after->Add(new TH2F("fTPCvsGlobalMult","TPC only vs Global track multiplicity;global;TPC only",500,0,2500,500,0,3500));//1
494   TH1::AddDirectory(adddirstatus);
495 }
496
497 //---------------------------------------------------------------//
498 void AliFlowEventCuts::Browse(TBrowser* b)
499 {
500   //some browsing capabilities
501   if (fQA) b->Add(fQA);
502 }
503
504 //---------------------------------------------------------------//
505 Long64_t AliFlowEventCuts::Merge(TCollection* list)
506 {
507   //merge
508   Int_t number=0;
509   AliFlowEventCuts* obj;
510   if (!list) return 0;
511   if (list->GetEntries()<1) return 0;
512   TIter next(list);
513   while ( (obj = dynamic_cast<AliFlowEventCuts*>(next())) )
514   {
515     if (obj==this) continue;
516     TList listwrapper;
517     listwrapper.Add(obj->GetQA());
518     fQA->Merge(&listwrapper);
519     number++;
520   }
521   return number;
522 }
523