updated macros for Sigma* analysis (M- Venaruzzo)
[u/mrichter/AliRoot.git] / PWGLF / RESONANCES / AliRsnCutV0.cxx
1 //
2 // Class AliRsnCutV0
3 //
4 // General implementation of a single cut strategy, which can be:
5 // - a value contained in a given interval  [--> IsBetween()   ]
6 // - a value equal to a given reference     [--> MatchesValue()]
7 //
8 // In all cases, the reference value(s) is (are) given as data members
9 // and each kind of cut requires a given value type (Int, UInt, Double),
10 // but the cut check procedure is then automatized and chosen thanks to
11 // an enumeration of the implemented cut types.
12 // At the end, the user (or any other point which uses this object) has
13 // to use the method IsSelected() to check if this cut has been passed.
14 //
15 // authors: Massimo Venaruzzo (massimo.venaruzzo@ts.infn.it)
16 // modified: Enrico Fragiacomo (enrico.fragiacomo@ts.infn.it)
17 //
18
19 #include <Riostream.h>
20 #include <TFormula.h>
21 #include <TBits.h>
22
23 #include "AliLog.h"
24 #include "AliESDtrackCuts.h"
25
26 #include "AliRsnEvent.h"
27 #include "AliRsnDaughter.h"
28 #include "AliRsnCutV0.h"
29
30 ClassImp(AliRsnCutV0)
31
32 //_________________________________________________________________________________________________
33 //AliRsnCutV0::AliRsnCutV0(const char *name, Int_t hypothesis) :
34 AliRsnCutV0::AliRsnCutV0(const char *name, Int_t hypothesis, AliPID::EParticleType pid, AliPID::EParticleType pid2) :
35    AliRsnCut(name, AliRsnTarget::kDaughter),
36    fHypothesis(0),
37    fMass(0.0),
38    fTolerance(0.01),
39    fMaxDCAVertex(0.3),
40    fMinCosPointAngle(0.95),
41    fMaxDaughtersDCA(0.5),
42    fMinTPCcluster(70),
43    fMaxRapidity(0.8),
44    fPID(pid),
45    fPID2(pid2),
46    fPIDCut1(0),
47    fPIDCut2(0),
48    fPIDCut3(0),
49    fESDtrackCuts(0x0),
50    fCutQuality(Form("%sDaughtersQuality", name)),
51    fAODTestFilterBit(5)
52 {
53 //
54 // Default constructor.
55 // Initializes all cuts in such a way that all of them are disabled.
56 //
57
58    SetHypothesis(hypothesis);
59 }
60
61 //_________________________________________________________________________________________________
62 AliRsnCutV0::AliRsnCutV0(const AliRsnCutV0 &copy) :
63    AliRsnCut(copy),
64    fHypothesis(copy.fHypothesis),
65    fMass(copy.fMass),
66    fTolerance(copy.fTolerance),
67    fMaxDCAVertex(copy.fMaxDCAVertex),
68    fMinCosPointAngle(copy.fMinCosPointAngle),
69    fMaxDaughtersDCA(copy.fMaxDaughtersDCA),
70    fMinTPCcluster(copy.fMinTPCcluster),
71    fMaxRapidity(copy.fMaxRapidity),
72    fPID(copy.fPID),
73    fPID2(copy.fPID2),
74    fPIDCut1(copy.fPIDCut1),
75    fPIDCut2(copy.fPIDCut2),
76    fPIDCut3(copy.fPIDCut3),
77    fESDtrackCuts(copy.fESDtrackCuts),
78    fCutQuality(copy.fCutQuality),
79    fAODTestFilterBit(copy.fAODTestFilterBit)
80 {
81 //
82 // Copy constructor.
83 // Just copy all data member values.:IsSelected: Object is not a V0 (RESONANCES/AliRsnCutV0.cxx:149)
84
85 //
86    fCutQuality.SetPtRange(0.15, 1E+20);
87    fCutQuality.SetEtaRange(-0.8, 0.8);
88    fCutQuality.SetDCARmax(0.05);
89    //fCutQuality.SetDCARPtFormula("0.0182+0.0350/pt^1.01");
90    fCutQuality.SetDCAZmax(2.0);
91    fCutQuality.SetSPDminNClusters(1);
92    fCutQuality.SetITSminNClusters(0);
93    fCutQuality.SetITSmaxChi2(1E+20);
94    fCutQuality.SetTPCminNClusters(fMinTPCcluster);
95    fCutQuality.SetTPCmaxChi2(4.0);
96    fCutQuality.SetRejectKinkDaughters();
97    fCutQuality.SetAODTestFilterBit(5);
98
99 }
100
101 //_________________________________________________________________________________________________
102 AliRsnCutV0 &AliRsnCutV0::operator=(const AliRsnCutV0 &copy)
103 {
104 //
105 // Assignment operator.
106 // Just copy all data member values.
107 //
108    if (this == &copy)
109      return *this;
110    fHypothesis = copy.fHypothesis;
111    fMass = copy.fMass;
112    fTolerance = copy.fTolerance;
113    fMaxDCAVertex = copy.fMaxDCAVertex;
114    fMinCosPointAngle = copy.fMinCosPointAngle;
115    fMaxDaughtersDCA = copy.fMaxDaughtersDCA;
116    fMinTPCcluster = copy.fMinTPCcluster;
117    fMaxRapidity = copy.fMaxRapidity;
118    fPID = copy.fPID;
119    fPID2 = copy.fPID2;
120    fPIDCut1 = copy.fPIDCut1;
121    fPIDCut2 = copy.fPIDCut2;
122    fPIDCut3 = copy.fPIDCut3;
123    fESDtrackCuts = copy.fESDtrackCuts;
124    fCutQuality = copy.fCutQuality;
125    fAODTestFilterBit = copy.fAODTestFilterBit;
126
127    return (*this);
128 }
129
130 //_________________________________________________________________________________________________
131 Bool_t AliRsnCutV0::IsSelected(TObject *object)
132 {
133 //:IsSelected: Object is not a V0 (RESONANCES/AliRsnCutV0.cxx:149)
134
135 // Cut checker.
136 // Checks the type of object being evaluated
137 // and then calls the appropriate sub-function (for ESD or AOD)
138 //
139
140    // coherence check
141    if (!TargetOK(object)) return kFALSE;
142
143    // check cast
144    AliESDv0 *v0esd = fDaughter->Ref2ESDv0();
145    AliAODv0 *v0aod = fDaughter->Ref2AODv0();
146    //cout << fDaughter->GetRef()->ClassName() << ' ' << v0esd << ' ' << v0aod << endl;
147
148    // operate depending on cast:IsSelected: Object is not a V0 (RESONANCES/AliRsnCutV0.cxx:149)
149
150    if (v0esd) {
151       return CheckESD(v0esd);
152    } else if (v0aod) {
153       return CheckAOD(v0aod);
154    } else {
155       AliDebugClass(1, "Object is not a V0");
156       return kFALSE;
157    }
158 }
159
160 //_________________________________________________________________________________________________
161 Bool_t AliRsnCutV0::CheckESD(AliESDv0 *v0)
162 {
163 //
164 // Check an ESD V0.
165 // This is done using the default track checker for ESD.
166 // It is declared static, not to recreate it every time.
167 //
168
169    AliDebugClass(1, "Check ESD");
170    if (v0->GetOnFlyStatus()) {
171       AliDebugClass(1, "Rejecting V0 in 'on fly' status");
172       return kFALSE; // if kTRUE, then this V0 is recontructed
173    }
174
175    // retrieve pointer to owner event
176    AliESDEvent *lESDEvent = fEvent->GetRefESD();
177    Double_t xPrimaryVertex = lESDEvent->GetPrimaryVertex()->GetX();
178    Double_t yPrimaryVertex = lESDEvent->GetPrimaryVertex()->GetY();
179    Double_t zPrimaryVertex = lESDEvent->GetPrimaryVertex()->GetZ();
180    AliDebugClass(2, Form("Primary vertex: %f %f %f", xPrimaryVertex, yPrimaryVertex, zPrimaryVertex));
181
182    // retrieve the V0 daughters
183    UInt_t lIdxPos      = (UInt_t) TMath::Abs(v0->GetPindex());
184    UInt_t lIdxNeg      = (UInt_t) TMath::Abs(v0->GetNindex());
185    AliESDtrack *pTrack = lESDEvent->GetTrack(lIdxPos);
186    AliESDtrack *nTrack = lESDEvent->GetTrack(lIdxNeg);
187
188    // check quality cuts
189    if (fESDtrackCuts) {
190       AliDebugClass(2, "Checking quality cuts");
191       if (!fESDtrackCuts->IsSelected(pTrack)) {
192          AliDebugClass(2, "Positive daughter failed quality cuts");
193          return kFALSE;
194       }
195       if (!fESDtrackCuts->IsSelected(nTrack)) {
196          AliDebugClass(2, "Negative daughter failed quality cuts");
197          return kFALSE;
198       }
199    }
200
201    // filter like-sign V0
202    if ( TMath::Abs( ((pTrack->GetSign()) - (nTrack->GetSign())) ) < 0.1) {
203       AliDebugClass(2, "Failed like-sign V0 check");
204       return kFALSE;
205    }
206
207
208    // check compatibility with expected species hypothesis
209    v0->ChangeMassHypothesis(fHypothesis);
210    if ((TMath::Abs(v0->GetEffMass() - fMass)) > fTolerance) {
211       AliDebugClass(2, "V0 is not in the expected inv mass range");
212       return kFALSE;
213    }
214
215    // topological checks
216    if (TMath::Abs(v0->GetD(xPrimaryVertex, yPrimaryVertex, zPrimaryVertex)) > fMaxDCAVertex) {
217       AliDebugClass(2, "Failed check on DCA to primary vertes");
218       return kFALSE;
219    }
220    if (TMath::Abs(v0->GetV0CosineOfPointingAngle()) < fMinCosPointAngle) {
221       AliDebugClass(2, "Failed check on cosine of pointing angle");
222       return kFALSE;
223    }
224    if (TMath::Abs(v0->GetDcaV0Daughters()) > fMaxDaughtersDCA) {
225       AliDebugClass(2, "Failed check on DCA between daughters");
226       return kFALSE;
227    }
228    if (TMath::Abs(v0->Y(fHypothesis)) > fMaxRapidity) {
229       AliDebugClass(2, "Failed check on V0 rapidity");
230       return kFALSE;
231    }
232
233
234    // check PID on proton or antiproton from V0
235
236    // check initialization of PID object
237    AliPIDResponse *pid = fEvent->GetPIDResponse();
238    if (!pid) {
239       AliFatal("NULL PID response");
240       return kFALSE;
241    }
242
243    // check if TOF is matched
244    // and computes all values used in the PID cut
245    //Bool_t   isTOFpos  = MatchTOF(ptrack);
246    //Bool_t   isTOFneg  = MatchTOF(ntrack);
247    Double_t pospTPC   = pTrack->GetTPCmomentum();
248    Double_t negpTPC   = nTrack->GetTPCmomentum();
249    //Double_t posp      = pTrack->P();
250    //Double_t negp      = nTrack->P();
251    Double_t posnsTPC   = TMath::Abs(pid->NumberOfSigmasTPC(pTrack, fPID));
252    Double_t posnsTPC2  = TMath::Abs(pid->NumberOfSigmasTPC(pTrack, fPID2));
253    //Double_t posnsTOF  = TMath::Abs(pid->NumberOfSigmasTOF(ptrack, fPID));
254    Double_t negnsTPC   = TMath::Abs(pid->NumberOfSigmasTPC(nTrack, fPID));
255    Double_t negnsTPC2  = TMath::Abs(pid->NumberOfSigmasTPC(nTrack, fPID2));
256    //Double_t negnsTOF  = TMath::Abs(pid->NumberOfSigmasTOF(ntrack, fPID));
257    Double_t maxTPC = 1E20;
258    Double_t maxTPC2 = 1E20;
259    //Double_t maxTOF = 1E20;
260
261    // applies the cut differently depending on the PID and the momentum
262
263    if(fHypothesis==kLambda0) {
264       //if (isTOFpos) {
265       // TPC: 5sigma cut for all
266       //if (posnsTPC > 5.0) return kFALSE;
267       // TOF: 3sigma
268       // maxTOF = 3.0;
269       //return (posnsTOF <= maxTOF);
270       //} else {
271       // TPC:
272       // below 600 MeV: 4sigma
273       // above 600 MeV: 3sigma
274
275       if (pospTPC <= 0.6 && fPID==AliPID::kProton)
276          maxTPC = fPIDCut1;
277       else if (pospTPC > 0.6 && fPID==AliPID::kProton)
278          maxTPC = fPIDCut2;
279       else
280          return kFALSE;
281
282       maxTPC2 = fPIDCut3;
283
284       if (! ((posnsTPC <= maxTPC) && (negnsTPC2 <= maxTPC2)) ) {
285          AliDebugClass(2, "Failed check on V0 PID");
286          return kFALSE;
287       }
288    }
289
290    //}
291
292    if(fHypothesis==kLambda0Bar) {
293       //if (isTOFneg) {
294       // TPC: 5sigma cut for all
295       //if (negnsTPC > 5.0) return kFALSE;
296       // TOF: 3sigma
297       // maxTOF = 3.0;
298       //return (negnsTOF <= maxTOF);
299       //} else {
300       // TPC:
301       // below 600 MeV: 4sigma
302       // above 600 MeV: 3sigma
303
304       if (negpTPC <= 0.6 && fPID==AliPID::kProton)
305          maxTPC = fPIDCut1;
306       else if (negpTPC > 0.6 && fPID==AliPID::kProton)
307          maxTPC = fPIDCut2;
308       else
309          return kFALSE;
310
311       maxTPC2 = fPIDCut3;
312
313       if(! ((negnsTPC <= maxTPC) && (posnsTPC2 <= maxTPC2)) ) {
314          AliDebugClass(2, "Failed check on V0 PID");
315          return kFALSE;
316       }
317    }
318    //}
319
320
321    // if we reach this point, all checks were successful
322    AliDebugClass(2, "Good V0 (hallelujah)");
323    return kTRUE;
324 }
325
326 //_________________________________________________________________________________________________
327 Bool_t AliRsnCutV0::CheckAOD(AliAODv0 *v0)
328 {
329 //
330 // Check an AOD V0.
331 // This is done doing directly all checks, since there is not
332 // an equivalend checker for AOD tracks
333 //
334
335    AliDebugClass(2, "Check AOD");
336    if (v0->GetOnFlyStatus()) {
337       AliDebugClass(2, "Rejecting V0 in 'on fly' status");
338       return kFALSE; // if kTRUE, then this V0 is recontructed
339    }
340
341
342
343    // retrieve pointer to owner event
344    AliAODEvent *lAODEvent = fEvent->GetRefAOD();
345    Double_t xPrimaryVertex = lAODEvent->GetPrimaryVertex()->GetX();
346    Double_t yPrimaryVertex = lAODEvent->GetPrimaryVertex()->GetY();
347    Double_t zPrimaryVertex = lAODEvent->GetPrimaryVertex()->GetZ();
348    AliDebugClass(2, Form("Primary vertex: %f %f %f", xPrimaryVertex, yPrimaryVertex, zPrimaryVertex));
349
350    // retrieve the V0 daughters
351    AliAODTrack *pTrack = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0));
352    AliAODTrack *nTrack = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1));
353
354
355    // check quality cuts
356    UInt_t  filtermapP = 9999;
357    UInt_t  filtermapN = 9999;
358    filtermapP = pTrack->GetFilterMap();
359    filtermapN = nTrack->GetFilterMap();
360
361
362    if ( !pTrack->TestFilterBit(fAODTestFilterBit)   ) {
363       AliDebugClass(2, Form("Positive daughter failed quality cuts filtermapP=%d",filtermapP));
364       return kFALSE;
365    }
366    if ( !nTrack->TestFilterBit(fAODTestFilterBit)   ) {
367       AliDebugClass(2, Form("Negative daughter failed quality cuts filtermapN=%d",filtermapN));
368       return kFALSE;
369    }
370
371
372
373    /*AliDebugClass(1, Form("fESDtrackCuts=%p",fESDtrackCuts));
374    if (fESDtrackCuts) { // use fESDtrackCuts to retrieve cuts values
375
376
377       AliDebugClass(2, "Checking quality cuts");
378
379
380
381       const Bool_t  CutAcceptKinkDaughters  = fESDtrackCuts->GetAcceptKinkDaughters(); // 0 = kFalse
382       const Float_t CutMaxChi2PerClusterTPC = fESDtrackCuts->GetMaxChi2PerClusterTPC();
383       const Int_t   CutMinNClusterTPC       = fESDtrackCuts->GetMinNClusterTPC();
384       const Bool_t  CutRequireTPCRefit      = fESDtrackCuts->GetRequireTPCRefit();
385       //AliDebugClass(2, Form("accept kink=%d maxchi2=%f minnclSTPC=%d requireTPCrefit=%d", CutAcceptKinkDaughters,
386       //        CutMaxChi2PerClusterTPC, CutMinNClusterTPC, CutRequireTPCRefit));
387       //AliDebugClass(2, Form("pTrack->TPCNcls=%d", pTrack->GetTPCNcls() ));
388       //AliDebugClass(2, Form("nTrack->TPCNcls=%d", nTrack->GetTPCNcls() ));
389
390
391       if(pTrack->GetTPCNcls() < CutMinNClusterTPC) {AliDebugClass(2, "Positive daughter not MinNclsTPC"); return kFALSE;}
392       if(nTrack->GetTPCNcls() < CutMinNClusterTPC) {AliDebugClass(2, "Negative daughter not MinNclsTPC"); return kFALSE;}
393
394    }*/
395
396    // filter like-sign V0
397    if ( TMath::Abs( ((pTrack->Charge()) - (nTrack->Charge())) ) < 0.1) {
398       AliDebugClass(2, "Failed like-sign V0 check");
399       return kFALSE;
400    }
401
402    // check compatibility with expected species hypothesis
403    Double_t mass = 0.0;
404    if(fHypothesis==kLambda0) {
405       mass = v0->MassLambda();
406    }
407    else if (fHypothesis==kLambda0Bar) {
408       mass = v0->MassAntiLambda();
409    }
410    if ((TMath::Abs(mass - fMass)) > fTolerance) {
411       AliDebugClass(2, Form("V0 is not in the expected inv mass range  Mass: %d %f %f", fHypothesis, fMass, mass));
412       return kFALSE;
413    }
414    AliDebugClass(2, Form("Mass: %d %f %f", fHypothesis, fMass, mass));
415
416
417    // topological checks
418    if (TMath::Abs(v0->DcaV0ToPrimVertex()) > fMaxDCAVertex) {
419       AliDebugClass(2, Form("Failed check on DCA to primary vertes dca=%f maxdca=%f",TMath::Abs(v0->DcaV0ToPrimVertex()),fMaxDCAVertex));
420       return kFALSE;
421    }
422    if (TMath::Abs( TMath::Cos(v0->OpenAngleV0()) ) < fMinCosPointAngle) {
423       AliDebugClass(2, "Failed check on cosine of pointing angle");
424       return kFALSE;
425    }
426    if (TMath::Abs(v0->DcaV0Daughters()) > fMaxDaughtersDCA) {
427       AliDebugClass(2, "Failed check on DCA between daughters");
428       return kFALSE;
429    }
430    if (TMath::Abs(v0->RapLambda()) > fMaxRapidity) {
431       AliDebugClass(2, "Failed check on V0 rapidity");
432       return kFALSE;
433    }
434
435
436
437
438
439    // check initialization of PID object
440    AliPIDResponse *pid = fEvent->GetPIDResponse();
441    if (!pid) {
442       AliFatal("NULL PID response");
443       return kFALSE;
444    }
445
446    // check if TOF is matched
447    // and computes all values used in the PID cut
448    //Bool_t   isTOFpos  = MatchTOF(ptrack);
449    //Bool_t   isTOFneg  = MatchTOF(ntrack);
450    Double_t pospTPC   = pTrack->GetTPCmomentum();
451    Double_t negpTPC   = nTrack->GetTPCmomentum();
452    //Double_t posp      = pTrack->P();
453    //Double_t negp      = nTrack->P();
454    Double_t posnsTPC   = TMath::Abs(pid->NumberOfSigmasTPC(pTrack, fPID));
455    Double_t posnsTPC2  = TMath::Abs(pid->NumberOfSigmasTPC(pTrack, fPID2));
456    //Double_t posnsTOF  = TMath::Abs(pid->NumberOfSigmasTOF(ptrack, fPID));
457    Double_t negnsTPC   = TMath::Abs(pid->NumberOfSigmasTPC(nTrack, fPID));
458    Double_t negnsTPC2  = TMath::Abs(pid->NumberOfSigmasTPC(nTrack, fPID2));
459    //Double_t negnsTOF  = TMath::Abs(pid->NumberOfSigmasTOF(ntrack, fPID));
460    Double_t maxTPC = 1E20;
461    Double_t maxTPC2 = 1E20;
462    //Double_t maxTOF = 1E20;
463
464    // applies the cut differently depending on the PID and the momentum
465
466    if(fHypothesis==kLambda0) {
467
468
469       //if (isTOFpos) {
470       // TPC: 5sigma cut for all
471       //if (posnsTPC > 5.0) return kFALSE;
472       // TOF: 3sigma
473       // maxTOF = 3.0;
474       //return (posnsTOF <= maxTOF);
475       //} else {
476       // TPC:
477       // below 600 MeV: 4sigma
478       // above 600 MeV: 3sigma
479
480       if (pospTPC <= 0.6 && fPID==AliPID::kProton)
481          maxTPC = fPIDCut1; // EF safer value to run on MC
482       else if (pospTPC > 0.6 && fPID==AliPID::kProton)
483          maxTPC = fPIDCut2; // EF safer value to run on MC
484       else
485          return kFALSE;
486
487       maxTPC2 = fPIDCut3; // EF safer value to run on MC
488
489       if (! ((posnsTPC <= maxTPC) && (negnsTPC2 <= maxTPC2)) ) {
490          AliDebugClass(2, "Failed check on V0 PID");
491          return kFALSE;
492       }
493    }
494
495    if(fHypothesis==kLambda0Bar) {
496
497       //if (isTOFneg) {
498       // TPC: 5sigma cut for all
499       //if (negnsTPC > 5.0) return kFALSE;
500       // TOF: 3sigma
501       // maxTOF = 3.0;
502       //return (negnsTOF <= maxTOF);
503       //} else {
504       // TPC:
505       // below 600 MeV: 4sigma
506       // above 600 MeV: 3sigma
507
508       if (negpTPC <= 0.6 && fPID==AliPID::kProton)
509          maxTPC = fPIDCut1; // EF safer value to run on MC
510       else if (negpTPC > 0.6 && fPID==AliPID::kProton)
511          maxTPC = fPIDCut2; // EF safer value to run on MC
512       else
513          return kFALSE;
514
515       maxTPC2 = fPIDCut3; // EF safer value to run on MC
516
517       if(! ((negnsTPC <= maxTPC) && (posnsTPC2 <= maxTPC2)) ) {
518          AliDebugClass(2, "Failed check on V0 PID");
519          return kFALSE;
520       }
521    }
522
523
524
525    // if we reach this point, all checks were successful
526    AliDebugClass(1, "Good AOD V0 (hallelujah)");
527    AliDebugClass(1, Form("Mass: %d %f %f %d %d", fHypothesis, fMass, mass, filtermapP, filtermapN));
528    return kTRUE;
529
530 }
531
532 //_________________________________________________________________________________________________
533 void AliRsnCutV0::Print(const Option_t *) const
534 {
535 //
536 // Print information on this cut
537 //
538 }