-add cut qa class
[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliDielectronV0Cuts.cxx
1 /*************************************************************************
2 * Copyright(c) 1998-2009, 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 //   Cut class providing cuts to V0 candidates                           //
18 //   Selection or deselection of V0 candiates can be done.               //
19 //                                                                       //
20 // Authors:                                                              //
21 //   Julian Book <Julian.Book@cern.ch>                                   //
22 /*
23
24 Class to provide some V0 selection using exactley the same cuts in ESDs as in AODs.
25 It implements the PID cut class AliDielectronPID and the standard AliDielectronVarCuts for
26 the configuration of leg respective pair cuts. These pair cuts are applied on the KFparticle
27 build by the legs.
28
29 Some QA cuts for the tracks are applied before the V0 pair is build. These cuts are:
30 AliDielectronVarCuts dauQAcuts1;
31 dauQAcuts1.AddCut(AliDielectronVarManager::kPt,            0.05, 100.0);
32 dauQAcuts1.AddCut(AliDielectronVarManager::kEta,          -0.9,    0.9);
33 dauQAcuts1.AddCut(AliDielectronVarManager::kNclsTPC,      50.0,  160.0);
34 dauQAcuts1.AddCut(AliDielectronVarManager::kTPCchi2Cl,     0.0,    4.0);
35 AliDielectronTrackCuts dauQAcuts2;
36 dauQAcuts2.SetRequireTPCRefit(kTRUE);
37
38
39
40 Example configuration:
41
42   AliDielectronV0Cuts *gammaV0Cuts = new AliDielectronV0Cuts("IsGamma","IsGamma");
43
44   // which V0 finder you want to use
45   gammaV0Cuts->SetV0finder(kOnTheFly);  // kAll(default), kOffline or kOnTheFly
46
47   // add some pdg codes (they are used then by the KF package and important for gamma conversions)
48   gammaV0Cuts->SetPdgCodes(22,11,11); // mother, daughter1 and 2
49
50   // add default PID cuts (defined in AliDielectronPID)
51   gammaV0Cuts->SetDefaultPID(16);
52
53   // add the pair cuts for V0 candidates
54   gammaV0Cuts->AddCut(AliDielectronVarManager::kCosPointingAngle, TMath::Cos(0.02),   1.0, kFALSE);
55   gammaV0Cuts->AddCut(AliDielectronVarManager::kChi2NDF,                       0.0,  10.0, kFALSE);
56   gammaV0Cuts->AddCut(AliDielectronVarManager::kLegDist,                       0.0,   0.25, kFALSE);
57   gammaV0Cuts->AddCut(AliDielectronVarManager::kR,                             3.0,  90.0, kFALSE);
58   gammaV0Cuts->AddCut(AliDielectronVarManager::kPsiPair,                       0.0,   0.05, kFALSE);
59   gammaV0Cuts->AddCut(AliDielectronVarManager::kM,                             0.0,   0.05, kFALSE);
60   gammaV0Cuts->AddCut(AliDielectronVarManager::kArmPt,                         0.0,   0.05, kFALSE);
61
62   // selection or rejection of V0 tracks
63   gammaV0Cuts->SetExcludeTracks(kTRUE);
64
65   // add the V0cuts directly to the track filter or to some cut group of it
66
67 */
68 //                                                                       //
69 ///////////////////////////////////////////////////////////////////////////
70
71
72 #include "AliDielectronV0Cuts.h"
73 #include "AliDielectronVarManager.h"
74 #include "AliDielectronTrackCuts.h"
75 #include "AliDielectronPID.h"
76 #include "AliESDv0.h"
77 #include "AliAODv0.h"
78
79 ClassImp(AliDielectronV0Cuts)
80
81
82 AliDielectronV0Cuts::AliDielectronV0Cuts() :
83   AliDielectronVarCuts(),
84   fV0TrackArr(0),
85   fExcludeTracks(kTRUE),
86   fV0finder(kAll),
87   fMotherPdg(0),
88   fNegPdg(0),
89   fPosPdg(0),
90   fPID(-1),
91   fOrbit(0),
92   fPeriod(0),
93   fBunchCross(0)
94 {
95   //
96   // Default costructor
97   //
98 }
99
100 //________________________________________________________________________
101 AliDielectronV0Cuts::AliDielectronV0Cuts(const char* name, const char* title) :
102   AliDielectronVarCuts(name,title),
103   fV0TrackArr(0),
104   fExcludeTracks(kTRUE),
105   fV0finder(kAll),
106   fMotherPdg(0),
107   fNegPdg(0),
108   fPosPdg(0),
109   fPID(-1),
110   fOrbit(0),
111   fPeriod(0),
112   fBunchCross(0)
113 {
114   //
115   // Named contructor
116   //
117 }
118
119 //________________________________________________________________________
120 AliDielectronV0Cuts::~AliDielectronV0Cuts()
121 {
122   //
123   // Destructor
124   //
125
126 }
127
128 //________________________________________________________________________
129 void AliDielectronV0Cuts::InitEvent(AliVTrack *trk)
130 {
131   //
132   // Init the V0 candidates
133   //
134
135   // take current event from the track
136   // TODO: this should be simplyfied by AliVTrack::GetEvent() as soon as implemented
137   const AliVEvent *ev=0;
138   if(trk->IsA() == AliAODTrack::Class())
139     ev=static_cast<const AliVEvent*>((static_cast<const AliAODTrack*>(trk))->GetAODEvent());
140   else if(trk->IsA() == AliESDtrack::Class())
141     ev=static_cast<const AliVEvent*>((static_cast<const AliESDtrack*>(trk))->GetESDEvent());
142   else
143     return;
144
145
146   // IsNewEvent
147   if(!ev) return;
148   if(!IsNewEvent(ev)) return;
149   //  printf(" Build V0 candidates according to the applied cuts \n");
150
151   // rest booleans
152   fV0TrackArr.ResetAllBits();
153
154   // basic quality cut, /*at least one*/ both of the V0 daughters has to fullfill
155   AliDielectronVarCuts dauQAcuts1;
156   dauQAcuts1.AddCut(AliDielectronVarManager::kPt,            0.05, 100.0);
157   dauQAcuts1.AddCut(AliDielectronVarManager::kEta,          -0.9,    0.9);
158   dauQAcuts1.AddCut(AliDielectronVarManager::kNclsTPC,      50.0,  160.0);
159   dauQAcuts1.AddCut(AliDielectronVarManager::kTPCchi2Cl,     0.0,    4.0);
160   AliDielectronTrackCuts dauQAcuts2;
161   //  dauQAcuts2.SetRequireITSRefit(kTRUE);
162   dauQAcuts2.SetRequireTPCRefit(kTRUE);
163   AliDielectronPID dauPIDcuts;
164   if(fPID>=0) dauPIDcuts.SetDefaults(fPID);
165
166   Int_t nV0s      = 0;
167   Int_t nV0stored = 0;
168   AliDielectronPair candidate;
169   candidate.SetPdgCode(fMotherPdg);
170
171   // ESD or AOD event
172   if(ev->IsA() == AliESDEvent::Class()) {
173     const AliESDEvent *esdev = static_cast<const AliESDEvent*>(ev);
174
175     //printf("there are %d V0s in the event \n",esdev->GetNumberOfV0s());
176     // loop over V0s
177     for (Int_t iv=0; iv<esdev->GetNumberOfV0s(); ++iv){
178       AliESDv0 *v = esdev->GetV0(iv);
179       if(!v) continue;
180       
181       // check the v0 finder
182       if( v->GetOnFlyStatus() && fV0finder==AliDielectronV0Cuts::kOffline  ) continue;
183       if(!v->GetOnFlyStatus() && fV0finder==AliDielectronV0Cuts::kOnTheFly ) continue;
184
185       // should we make use of AliESDv0Cuts::GetPdgCode() to preselect candiadtes, e.g.:
186       // if(fMotherPdg!=v->GetPdgCode()) continue;
187
188       AliESDtrack *trNeg=esdev->GetTrack(v->GetIndex(0));
189       AliESDtrack *trPos=esdev->GetTrack(v->GetIndex(1));
190       if(!trNeg || !trPos){
191         printf("Error: Couldn't get V0 daughter: %p - %p\n",trNeg,trPos);
192         continue;
193       }
194
195       // protection against LS v0s
196       if(trNeg->Charge() == trPos->Charge()) continue;
197
198       // PID default cuts
199       if(fPID>=0) {
200         if( !dauPIDcuts.IsSelected(trNeg) ) continue;
201         if( !dauPIDcuts.IsSelected(trPos) ) continue;
202       }
203
204       // basic track cuts
205       if( !dauQAcuts2.IsSelected(trNeg) ) continue;
206       if( !dauQAcuts2.IsSelected(trPos) ) continue;
207       if( !dauQAcuts1.IsSelected(trNeg) ) continue;
208       if( !dauQAcuts1.IsSelected(trPos) ) continue;
209
210       if(fMotherPdg==22) candidate.SetGammaTracks(trNeg, 11, trPos, 11);
211       else candidate.SetTracks(trNeg, (trNeg->Charge()<0?fNegPdg:fPosPdg), trPos, (trPos->Charge()<0?fNegPdg:fPosPdg));
212       // eventually take the external trackparam and build the KFparticles by hand (see AliESDv0::GetKFInfo)
213       // the folowing is not needed, because the daughters where used in the v0 vertex fit (I guess)
214       //      AliKFVertex v0vtx = *v;
215       //      candidate.SetProductionVertex(v0vtx);
216
217       if(AliDielectronVarCuts::IsSelected(&candidate)) {
218         nV0s++;
219         //printf(" gamma found for vtx %p dau1id %d dau2id %d \n",v,trNeg->GetID(),trPos->GetID());
220         fV0TrackArr.SetBitNumber(trNeg->GetID());
221         fV0TrackArr.SetBitNumber(trPos->GetID());
222       }
223     }
224
225   }
226   else if(ev->IsA() == AliAODEvent::Class()) {
227     const AliAODEvent *aodEv = static_cast<const AliAODEvent*>(ev);
228     if(!aodEv->GetV0s()) return; // protection for nano AODs
229
230     // loop over vertices
231     for (Int_t ivertex=0; ivertex<aodEv->GetNumberOfV0s(); ++ivertex){
232       AliAODv0 *v=aodEv->GetV0(ivertex);
233       if(!v) continue;
234
235       // check the v0 finder
236       if( v->GetOnFlyStatus() && fV0finder==AliDielectronV0Cuts::kOffline  ) continue;
237       if(!v->GetOnFlyStatus() && fV0finder==AliDielectronV0Cuts::kOnTheFly ) continue;
238
239       AliAODTrack *trNeg=dynamic_cast<AliAODTrack*>(v->GetDaughter(0));
240       AliAODTrack *trPos=dynamic_cast<AliAODTrack*>(v->GetDaughter(1));
241       if(!trNeg || !trPos){
242         printf("Error: Couldn't get V0 daughter: %p - %p\n",trNeg,trPos);
243         continue;
244       }
245       nV0stored++;
246
247       // protection against LS v0s
248       if(trNeg->Charge() == trPos->Charge()) continue;
249
250       // PID default cuts
251       if(fPID>=0) {
252         if( !dauPIDcuts.IsSelected(trNeg) ) continue;
253         if( !dauPIDcuts.IsSelected(trPos) ) continue;
254       }
255
256       // basic track cuts
257       if( !dauQAcuts2.IsSelected(trNeg) ) continue;
258       if( !dauQAcuts2.IsSelected(trPos) ) continue;
259       if( !dauQAcuts1.IsSelected(trNeg) ) continue;
260       if( !dauQAcuts1.IsSelected(trPos) ) continue;
261
262       AliKFVertex v0vtx = *(v->GetSecondaryVtx());
263       if(fMotherPdg==22) candidate.SetGammaTracks(trNeg, 11, trPos, 11);
264       else candidate.SetTracks(trNeg, (trNeg->Charge()<0?fNegPdg:fPosPdg), trPos, (trPos->Charge()<0?fNegPdg:fPosPdg));
265       candidate.SetProductionVertex(v0vtx);
266
267       if(AliDielectronVarCuts::IsSelected(&candidate)) {
268         nV0s++;
269         //printf(" gamma found for vtx %p dau1id %d dau2id %d \n",v,trNeg->GetID(),trPos->GetID());
270         fV0TrackArr.SetBitNumber(trNeg->GetID());
271         fV0TrackArr.SetBitNumber(trPos->GetID());
272       }
273     }
274     //printf("there are %d V0s in the event \n",nV0stored);
275   }
276   else
277     return;
278
279   //  printf(" Number of V0s candiates found %d/%d \n",nV0s,nV0stored);
280
281 }
282 //________________________________________________________________________
283 Bool_t AliDielectronV0Cuts::IsSelected(TObject* track)
284 {
285   //
286   // Make cut decision
287   //
288
289   if(!track) return kFALSE;
290
291   AliVTrack *vtrack = static_cast<AliVTrack*>(track);
292   InitEvent(vtrack);
293   //printf(" track ID %d selected result %d %d \n",vtrack->GetID(),(fV0TrackArr.TestBitNumber(vtrack->GetID())),fExcludeTracks);
294   return ( (fV0TrackArr.TestBitNumber(vtrack->GetID()))^fExcludeTracks );
295 }
296
297 //________________________________________________________________________
298 Bool_t AliDielectronV0Cuts::IsNewEvent(const AliVEvent *ev)
299 {
300   //
301   // check weather we process a new event
302   //
303
304   //  printf(" current ev %d %d %d \n",fBunchCross, fOrbit, fPeriod);
305   //  printf(" new event %p %d %d %d \n",ev, ev->GetBunchCrossNumber(), ev->GetOrbitNumber(), ev->GetPeriodNumber());
306
307   if( fBunchCross == ev->GetBunchCrossNumber() ) {
308     if( fOrbit == ev->GetOrbitNumber() )         {
309       if( fPeriod == ev->GetPeriodNumber() )     {
310         return kFALSE;
311       }
312     }
313   }
314
315   fBunchCross = ev->GetBunchCrossNumber();
316   fOrbit      = ev->GetOrbitNumber();
317   fPeriod     = ev->GetPeriodNumber();
318   return kTRUE;
319 }