ad68d85c2235d23f5e0b550e4f2be0432c5a348d
[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
229     // loop over vertices
230     for (Int_t ivertex=0; ivertex<aodEv->GetNumberOfV0s(); ++ivertex){
231       AliAODv0 *v=aodEv->GetV0(ivertex);
232       if(!v) continue;
233
234       // check the v0 finder
235       if( v->GetOnFlyStatus() && fV0finder==AliDielectronV0Cuts::kOffline  ) continue;
236       if(!v->GetOnFlyStatus() && fV0finder==AliDielectronV0Cuts::kOnTheFly ) continue;
237
238       AliAODTrack *trNeg=dynamic_cast<AliAODTrack*>(v->GetDaughter(0));
239       AliAODTrack *trPos=dynamic_cast<AliAODTrack*>(v->GetDaughter(1));
240       if(!trNeg || !trPos){
241         printf("Error: Couldn't get V0 daughter: %p - %p\n",trNeg,trPos);
242         continue;
243       }
244       nV0stored++;
245
246       // protection against LS v0s
247       if(trNeg->Charge() == trPos->Charge()) continue;
248
249       // PID default cuts
250       if(fPID>=0) {
251         if( !dauPIDcuts.IsSelected(trNeg) ) continue;
252         if( !dauPIDcuts.IsSelected(trPos) ) continue;
253       }
254
255       // basic track cuts
256       if( !dauQAcuts2.IsSelected(trNeg) ) continue;
257       if( !dauQAcuts2.IsSelected(trPos) ) continue;
258       if( !dauQAcuts1.IsSelected(trNeg) ) continue;
259       if( !dauQAcuts1.IsSelected(trPos) ) continue;
260
261       AliKFVertex v0vtx = *(v->GetSecondaryVtx());
262       if(fMotherPdg==22) candidate.SetGammaTracks(trNeg, 11, trPos, 11);
263       else candidate.SetTracks(trNeg, (trNeg->Charge()<0?fNegPdg:fPosPdg), trPos, (trPos->Charge()<0?fNegPdg:fPosPdg));
264       candidate.SetProductionVertex(v0vtx);
265
266       if(AliDielectronVarCuts::IsSelected(&candidate)) {
267         nV0s++;
268         //printf(" gamma found for vtx %p dau1id %d dau2id %d \n",v,trNeg->GetID(),trPos->GetID());
269         fV0TrackArr.SetBitNumber(trNeg->GetID());
270         fV0TrackArr.SetBitNumber(trPos->GetID());
271       }
272     }
273     //printf("there are %d V0s in the event \n",nV0stored);
274   }
275   else
276     return;
277
278   //  printf(" Number of V0s candiates found %d/%d \n",nV0s,nV0stored);
279
280 }
281 //________________________________________________________________________
282 Bool_t AliDielectronV0Cuts::IsSelected(TObject* track)
283 {
284   //
285   // Make cut decision
286   //
287
288   if(!track) return kFALSE;
289
290   AliVTrack *vtrack = static_cast<AliVTrack*>(track);
291   InitEvent(vtrack);
292   //printf(" track ID %d selected result %d %d \n",vtrack->GetID(),(fV0TrackArr.TestBitNumber(vtrack->GetID())),fExcludeTracks);
293   return ( (fV0TrackArr.TestBitNumber(vtrack->GetID()))^fExcludeTracks );
294 }
295
296 //________________________________________________________________________
297 Bool_t AliDielectronV0Cuts::IsNewEvent(const AliVEvent *ev)
298 {
299   //
300   // check weather we process a new event
301   //
302
303   //  printf(" current ev %d %d %d \n",fBunchCross, fOrbit, fPeriod);
304   //  printf(" new event %p %d %d %d \n",ev, ev->GetBunchCrossNumber(), ev->GetOrbitNumber(), ev->GetPeriodNumber());
305
306   if( fBunchCross == ev->GetBunchCrossNumber() ) {
307     if( fOrbit == ev->GetOrbitNumber() )         {
308       if( fPeriod == ev->GetPeriodNumber() )     {
309         return kFALSE;
310       }
311     }
312   }
313
314   fBunchCross = ev->GetBunchCrossNumber();
315   fOrbit      = ev->GetOrbitNumber();
316   fPeriod     = ev->GetPeriodNumber();
317   return kTRUE;
318 }