adding track cut filter. Only add used ME tracks to ME
[u/mrichter/AliRoot.git] / PWGGA / GammaConv / ConvCorrelations / AliConversionTrackCuts.cxx
1 /**************************************************************************
2  * This file is property of and copyright by the ALICE HLT Project        *
3  * ALICE Experiment at CERN, All rights reserved.                         *
4  *                                                                        *
5  * Primary Author: Svein Lindal <slindal@fys.uio.no>                      *
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 /// @file   AliConversionTrackCuts.cxx
17 /// @author Svein Lindal
18 /// @brief  Base class for analysation of conversion particle - track correlations
19
20
21 #include "AliConversionTrackCuts.h"
22 //#include "AliAODTrack.h"
23 #include "AliAODEvent.h"
24 #include <TFormula.h>
25
26
27
28 using namespace std;
29 ClassImp(AliConversionTrackCuts)
30
31 //________________________________________________________________________
32 AliConversionTrackCuts::AliConversionTrackCuts() : 
33 AliAnalysisCuts(),
34   fFlagsOn(0x0),
35   fFlagsOff(0x0),
36   fRejectKinkDaughters(kTRUE),
37   fDCARfixed(kTRUE),
38   fDCARptFormula(""),
39   fDCARmax(1E20),
40   fDCAZfixed(kTRUE),
41   fDCAZptFormula(""),
42   fDCAZmax(1E20),
43   fSPDminNClusters(0),
44   fITSminNClusters(0),
45   fITSmaxChi2(1E20),
46   fTPCminNClusters(0),
47   fTPCmaxChi2(1E20),
48   fAODTestFilterBit(-1)
49 {
50   //Constructor
51 }
52 //________________________________________________________________________
53 AliConversionTrackCuts::AliConversionTrackCuts(TString name, TString title = "title") : 
54   AliAnalysisCuts(name, title),
55   fFlagsOn(0x0),
56   fFlagsOff(0x0),
57   fRejectKinkDaughters(kTRUE),
58   fDCARfixed(kTRUE),
59   fDCARptFormula(""),
60   fDCARmax(1E20),
61   fDCAZfixed(kTRUE),
62   fDCAZptFormula(""),
63   fDCAZmax(1E20),
64   fSPDminNClusters(0),
65   fITSminNClusters(0),
66   fITSmaxChi2(1E20),
67   fTPCminNClusters(0),
68   fTPCmaxChi2(1E20),
69   fAODTestFilterBit(-1)
70
71 {
72   //Constructor
73 }
74
75
76 //________________________________________________________________________________
77 // AliConversionTrackCuts::~AliConversionTrackCuts() {
78 //   ///destructor
79 // }
80
81 Bool_t AliConversionTrackCuts::AcceptTrack(AliAODTrack * track, AliAODEvent * aodEvent) {
82   // Check an AOD track.
83 // This is done doing directly all checks, since there is not
84 // an equivalend checker for AOD tracks
85 //
86
87    // try to retrieve the reference AOD event
88    // AliAODEvent *aodEvent = 0x0;
89    // if (fEvent) aodEvent = fEvent->GetRefAOD();
90    // if (!aodEvent) {
91    //    AliError("AOD reference event is not initialized!");
92    //    return kFALSE;
93    // }
94
95    // step #0: check SPD and ITS clusters
96    Int_t nSPD = 0;
97    nSPD  = TESTBIT(track->GetITSClusterMap(), 0);
98    nSPD += TESTBIT(track->GetITSClusterMap(), 1);
99    if (nSPD < fSPDminNClusters) {
100       AliDebug(AliLog::kDebug + 2, "Not enough SPD clusters in this track. Rejected");
101       return kFALSE;
102    }
103
104    // step #1: check number of clusters in TPC
105    if (track->GetTPCNcls() < fTPCminNClusters) {
106       AliDebug(AliLog::kDebug + 2, "Too few TPC clusters. Rejected");
107       return kFALSE;
108    }
109    if (track->GetITSNcls() < fITSminNClusters) {
110       AliDebug(AliLog::kDebug + 2, "Too few ITS clusters. Rejected");
111       return kFALSE;
112    }
113
114    // step #2: check chi square
115    if (track->Chi2perNDF() > fTPCmaxChi2) {
116       AliDebug(AliLog::kDebug + 2, "Bad chi2. Rejected");
117       return kFALSE;
118    }
119    if (track->Chi2perNDF() > fITSmaxChi2) {
120       AliDebug(AliLog::kDebug + 2, "Bad chi2. Rejected");
121       return kFALSE;
122    }
123
124    // step #3: reject kink daughters
125    AliAODVertex *vertex = track->GetProdVertex();
126    if (vertex && fRejectKinkDaughters) {
127       if (vertex->GetType() == AliAODVertex::kKink) {
128          AliDebug(AliLog::kDebug + 2, "Kink daughter. Rejected");
129          return kFALSE;
130       }
131    }
132
133    // step #4: DCA cut (transverse)
134    Double_t b[2], cov[3];
135    vertex = aodEvent->GetPrimaryVertex();
136    if (!vertex) {
137       AliDebug(AliLog::kDebug + 2, "NULL vertex");
138       return kFALSE;
139    }
140    if (!track->PropagateToDCA(vertex, aodEvent->GetMagneticField(), kVeryBig, b, cov)) {
141       AliDebug(AliLog::kDebug + 2, "Failed propagation to vertex");
142       return kFALSE;
143    }
144    // if the DCA cut is not fixed, compute current value
145    if (!fDCARfixed) {
146       static TString str(fDCARptFormula);
147       str.ReplaceAll("pt", "x");
148       static const TFormula dcaXY(Form("%s_dcaXY", GetName()), str.Data());
149       fDCARmax = dcaXY.Eval(track->Pt());
150    }
151    // check the cut
152    if (TMath::Abs(b[0]) > fDCARmax) {
153       AliDebug(AliLog::kDebug + 2, "Too large transverse DCA");
154       return kFALSE;
155    }
156
157    // step #5: DCA cut (longitudinal)
158    // the DCA has already been computed above
159    // if the DCA cut is not fixed, compute current value
160    if (!fDCAZfixed) {
161       static TString str(fDCAZptFormula);
162       str.ReplaceAll("pt", "x");
163       static const TFormula dcaZ(Form("%s_dcaXY", GetName()), str.Data());
164       fDCAZmax = dcaZ.Eval(track->Pt());
165    }
166    // check the cut
167    if (TMath::Abs(b[1]) > fDCAZmax) {
168       AliDebug(AliLog::kDebug + 2, "Too large longitudinal DCA");
169       return kFALSE;
170    }
171
172    // step #6: check eta/pt range
173    if (track->Eta() < fEta[0] || track->Eta() > fEta[1]) {
174       AliDebug(AliLog::kDebug + 2, "Outside ETA acceptance");
175       return kFALSE;
176    }
177    if (track->Pt() < fPt[0] || track->Pt() > fPt[1]) {
178       AliDebug(AliLog::kDebug + 2, "Outside PT acceptance");
179       return kFALSE;
180    }
181
182    // if we are here, all cuts were passed and no exit point was got
183    return kTRUE;
184 }
185
186 //_________________________________________________________________________________________________
187 void AliConversionTrackCuts::Print(const Option_t *) const
188 {
189 //
190 // Print information on this cut
191 //
192
193    AliInfo(Form("Cut name                : %s", GetName()));
194    AliInfo(Form("Required flags (off, on): %lx %lx", fFlagsOn, fFlagsOff));
195    AliInfo(Form("Ranges in eta, pt       : %.2f - %.2f, %.2f - %.2f", fEta[0], fEta[1], fPt[0], fPt[1]));
196    AliInfo(Form("Kink daughters are      : %s", (fRejectKinkDaughters ? "rejected" : "accepted")));
197    AliInfo(Form("TPC requirements        : min. cluster = %d, max chi2 = %f", fTPCminNClusters, fTPCmaxChi2));
198    AliInfo(Form("ITS requirements        : min. cluster = %d (all), %d (SPD), max chi2 = %f", fITSminNClusters, fSPDminNClusters, fITSmaxChi2));
199
200    if (fDCARfixed) {
201          AliInfo(Form("DCA r cut               : fixed to %f cm", fDCARmax));
202    } else {
203          AliInfo(Form("DCA r cut formula       : %s", fDCARptFormula.Data()));
204    }
205    
206    if (fDCAZfixed) {
207          AliInfo(Form("DCA z cut               : fixed to %f cm", fDCAZmax));
208    } else {
209          AliInfo(Form("DCA z cut formula       : %s", fDCAZptFormula.Data()));
210    }
211    
212    
213 }
214
215
216
217
218
219
220
221  
222 //________________________________________________________________________________
223 // void AliConversionTrackCuts::CreateHistograms() {
224 //   CreateBaseHistograms();
225 // }
226
227
228 // ///________________________________________________________________________________
229 // void AliConversionTrackCuts::SetUpDefaultBins() {
230 //   //Set up default bins
231 //   Double_t ptbins[19] = {0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0, 6.0, 8.0, 10.0, 12.5, 15, 20, 25, 30, 50, 100};
232 //   fAxisdEta.Set(160, -1.6, 1.6);
233 //   fAxisdEta.SetNameTitle("dEta", "delta eta");
234
235 //   fAxisdPhi.Set(64, -TMath::PiOver2(), 3*TMath::PiOver2());
236 //   fAxisdPhi.SetNameTitle("dPhi", "delta Phi");
237
238 //   fAxistPt.Set(18, ptbins);
239 //   fAxistPt.SetNameTitle("tPt", "trigger Pt");
240
241 //   fAxiscPt.Set(18, ptbins);
242 //   fAxiscPt.SetNameTitle("cPt", "track Pt");
243
244 //   fAxisIso.Set(3, -0.5, 2.5);
245 //   fAxisIso.SetNameTitle("iso", "isolation");
246
247 //   fAxesList.AddAt(&fAxisdEta, 0);
248 //   fAxesList.AddAt(&fAxisdPhi, 1);
249 //   fAxesList.AddAt(&fAxistPt, 2);
250 //   fAxesList.AddAt(&fAxiscPt, 3);
251 //   fAxesList.AddAt(&fAxisIso, 4);
252
253 //   fAxisMEEta.Set(160, -0.8, 0.8);
254 //   fAxisMEEta.SetNameTitle("eta", "eta");
255   
256 //   fAxisMEPhi.Set(64, 0, TMath::TwoPi());
257 //   fAxisMEPhi.SetNameTitle("phi", "phi");
258
259 //   fTrackAxisList.AddAt(&fAxisMEEta, 0);
260 //   fTrackAxisList.AddAt(&fAxisMEPhi, 1);
261 //   fTrackAxisList.AddAt(&fAxistPt, 2);
262 //   fTrackAxisList.AddAt(&fAxiscPt, 3);
263 //   fTrackAxisList.AddAt(&fAxisIso, 4);
264
265 //   fTrigAxisList.AddAt(&fAxisMEEta, 0);
266 //   fTrigAxisList.AddAt(&fAxisMEPhi, 1);
267 //   fTrigAxisList.AddAt(&fAxistPt, 2);
268 //   fTrigAxisList.AddAt(&fAxisIso, 3);
269
270 //   for(int iIso = 0; iIso < 2; iIso++) {
271 //     fHNTriggers[iIso] = NULL;
272 //   }
273 // }
274
275
276 // //________________________________________________________________________
277 // void AliConversionTrackCuts::CreateBaseHistograms() {
278 //   //Create histograms add, to outputlis
279
280 //   cout << "Creating histograms for "<< GetName() << endl;
281
282 //   fHistograms = new TList();
283 //   fHistograms->SetOwner(kTRUE);
284 //   fHistograms->SetName(fName);
285
286 //   for(int iIso = 0; iIso < 2; iIso++) {
287
288 //     fHNTriggers[iIso] = new TH1F(Form("%s_%s_fNTriggers", fName.Data(), (iIso==0)?"nonIso":"isolated"), 
289 //                                                               Form("%s_%s_fNTriggers", fName.Data(), (iIso==0)?"nonIso":"isolated"), 
290 //                                                               fAxistPt.GetNbins(), fAxistPt.GetXbins()->GetArray());
291 //     fHNTriggers[iIso]->Sumw2();
292 //     fHistograms->Add(fHNTriggers[iIso]);
293   
294 //   }
295
296 //   fCorrSparse = CreateSparse(GetName(), GetTitle(), &fAxesList);
297 //   fHistograms->Add(fCorrSparse);
298
299 //   fTrackSparse = CreateSparse(Form("%s_%s", GetName(), "METrack"), Form("%s %s", GetTitle(), "ME Tracks"), &fTrackAxisList);
300 //   fHistograms->Add(fTrackSparse);
301
302 //   fTrigSparse = CreateSparse(Form("%s_%s", GetName(), "METrig"), Form("%s %s", GetTitle(), "ME Triggers"), &fTrigAxisList);
303 //   fHistograms->Add(fTrigSparse);
304
305 // }
306
307 // ///________________________________________________________________________
308 // THnSparseF * AliConversionTrackCuts::CreateSparse(TString nameString, TString titleString, TList * axesList) {
309 //   //Create sparse
310 //   const Int_t dim = axesList->GetSize();
311   
312 //   cout << nameString << " " << titleString << " " <<   "    dimesion: " << dim << endl;
313
314 //   TAxis * axes[dim];
315 //   Int_t   bins[dim];
316 //   Double_t min[dim];
317 //   Double_t max[dim];
318
319 //   for(Int_t i = 0; i<dim; i++) {
320 //      TAxis * axis = dynamic_cast<TAxis*>(axesList->At(i));
321 //      if(axis) axes[i] = axis;
322 //      else {
323 //        cout << "AliAnalysisTaskdPhi::CreateSparse: Error error, all the axes are not present in axis list" << endl;
324 //        return NULL;
325 //      }
326 //   }
327
328 //   for(Int_t i = 0; i<dim; i++) {
329 //      cout << axes[i]->GetTitle() << endl;
330 //      bins[i] = axes[i]->GetNbins(); 
331 //      min[i] = axes[i]->GetBinLowEdge(1);
332 //      max[i] = axes[i]->GetBinUpEdge(axes[i]->GetNbins());
333 //   }
334
335 //   THnSparseF * sparse = new THnSparseF(Form("%s", nameString.Data()), 
336 //                                                                         Form("%s", titleString.Data()), 
337 //                                                                         dim, bins, min, max);
338   
339 //   for(Int_t i = 0; i<dim; i++) {
340 //      sparse->GetAxis(i)->SetNameTitle(axes[i]->GetName(), axes[i]->GetTitle() );
341 //      if(axes[i]->GetXbins()->GetSize() > 0) {
342 //        sparse->SetBinEdges(i, axes[i]->GetXbins()->GetArray() );
343 //      }
344 //   }
345 //   return sparse;
346 // }
347
348
349 // ///____________________________________________________________________________
350 // // void AliConversionTrackCuts::FillTriggerCounters(Float_t tPt, Bool_t isolated){ 
351 // //   //Fill histogram with trigger counters
352
353 // //   fHNTriggers[0]->Fill(tPt);
354   
355 // //   if(isolated) {
356 // //     fHNTriggers[isolated]->Fill(tPt);
357     
358 // //   }
359 // // }
360
361 // // ///_____________________________________________________________________________
362 // // void AliConversionTrackCuts::FillHistograms(Float_t tPt, Float_t cPt, Float_t dPhi, Float_t dEta, Bool_t isolated) {
363 // //   //Fill histograms
364
365 // //   if(dEta) { ;}
366 // //   //fHdPhi[0]->Fill(tPt, cPt, dPhi);
367 // //   if(isolated) {
368 // //     //fHdPhi[isolated]->Fill(tPt, cPt, dPhi);
369 // //   }
370 // // }
371
372 // //_______________________________________________________________________________
373
374 // void AliConversionTrackCuts::PrintStatistics()  { 
375 //   //Print some statistics between each file
376 //   for(Int_t i = 1; i <= fHNTriggers[0]->GetNbinsX(); i++) {
377 //     Int_t nTrig = (Int_t) fHNTriggers[0]->GetBinContent(i+1);
378 //     cout << "triggers: " << nTrig << endl;
379
380 //   }
381 // }
382
383
384 // //_______________________________________________________________________________
385 // void AliConversionTrackCuts::FillTriggerCounters(const AliAODConversionParticle * particle, Bool_t leading) {
386 //   fHNTriggers[leading]->Fill(particle->Pt());
387 // }
388
389
390 // //________________________________________________________________
391 // void AliConversionTrackCuts::CorrelateWithTracks(AliAODConversionParticle * particle, TObjArray * tracks, Int_t const tIDs[4], Int_t isolated = 0) {
392 //   //Correlate particle with tracks
393
394 //   FillTriggerCounters(particle, isolated);
395
396 //   Int_t nDim = fAxesList.GetSize();
397 //   Double_t dphivalues[nDim];
398 //   Double_t trackValues[nDim];
399 //   Double_t trigValues[nDim - 1];
400
401 //   for(int ij = 0; ij < tracks->GetEntriesFast(); ij++) {
402 //      AliVTrack * track = static_cast<AliVTrack*>(tracks->UncheckedAt(ij));
403 //      Int_t tid = track->GetID();
404
405 //      if((tid > 0) && (tid == tIDs[0] || tid == tIDs[1] || tid == tIDs[2] || tid == tIDs[3]) ) {
406 //        continue;
407 //      }
408         
409 //      dphivalues[0] = particle->Eta() - track->Eta();
410 //      dphivalues[1] = GetDPhi(particle->Phi() - track->Phi());
411 //      dphivalues[2] = particle->Pt();
412 //      dphivalues[3] = track->Pt();
413 //      dphivalues[4] = isolated;
414
415 //      trackValues[0] = track->Eta();
416 //      trackValues[1] = track->Phi();
417 //      trackValues[2] = particle->Pt();
418 //      trackValues[3] = track->Pt();
419 //      trackValues[4] = isolated;
420
421 //      trigValues[0] = particle->Eta();
422 //      trigValues[1] = particle->Phi();
423 //      trigValues[2] = particle->Pt();
424 //      trigValues[4] = isolated;
425
426 //      if(nDim > 4) {
427 //        dphivalues[5] = particle->M();
428 //        trackValues[5] = particle->M();
429 //        trigValues[4] = particle->M();
430 //      }
431         
432 //      fCorrSparse->Fill(dphivalues);
433 //      fTrackSparse->Fill(trackValues);
434 //      fTrigSparse->Fill(trigValues);
435 //   }
436 // }
437