set reco param on an event by event basis
[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDefficiency.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: AliTRDefficiency.cxx 27496 2008-07-22 08:35:45Z cblume $ */
17
18 ////////////////////////////////////////////////////////////////////////////
19 //                                                                        //
20 //  Reconstruction QA                                                     //
21 //                                                                        //
22 //  Authors:                                                              //
23 //    Markus Fasel <M.Fasel@gsi.de>                                       //
24 //                                                                        //
25 ////////////////////////////////////////////////////////////////////////////
26
27 #include <TClonesArray.h>
28 #include <TObjArray.h>
29 #include <TProfile.h>
30 #include <TPad.h>
31 #include <TLegend.h>
32 #include "TTreeStream.h"
33
34 #include "AliPID.h"
35 #include "AliESDtrack.h"
36 #include "AliTrackReference.h"
37 #include "AliExternalTrackParam.h"
38 #include "AliTracker.h"
39 #include "AliAnalysisManager.h"
40
41 #include "Cal/AliTRDCalPID.h"
42 #include "AliTRDefficiency.h"
43 #include "info/AliTRDtrackInfo.h"
44
45 ClassImp(AliTRDefficiency)
46
47 //____________________________________________________________________
48 AliTRDefficiency::AliTRDefficiency()
49   :AliTRDrecoTask()
50   ,fMissed(NULL)
51 {
52   //
53   // Default constructor
54   //
55   SetNameTitle("TRDefficiency", "TRD barrel tracking efficiency checker");
56 }
57
58 //____________________________________________________________________
59 AliTRDefficiency::AliTRDefficiency(char* name)
60   :AliTRDrecoTask(name, "TRD barrel tracking efficiency checker")
61   ,fMissed(NULL)
62 {
63   //
64   // Default constructor
65   //
66 }
67
68 //____________________________________________________________________
69 AliTRDefficiency::~AliTRDefficiency()
70 {
71   // Destructor
72   if(fMissed){
73     fMissed->Delete();
74     delete fMissed;
75   }
76 }
77
78 //____________________________________________________________________
79 void  AliTRDefficiency::UserCreateOutputObjects()
80 {
81   //
82   // Create output objects
83   //
84
85   const Int_t nbins = AliTRDCalPID::kNMom;
86   Float_t xbins[nbins+1] = {.5, .7, .9, 1.3, 1.7, 2.4, 3.5, 4.5, 5.5, 7., 9., 11.};
87
88   TH1 *h = NULL;
89   fContainer = new TObjArray(); fContainer->SetOwner();
90   for(Int_t is=0; is<AliPID::kSPECIES; is++){
91     fContainer->Add(h = new TProfile(Form("h%s", AliTRDCalPID::GetPartSymb(is)), AliPID::ParticleShortName(is), nbins, xbins));
92     h->SetLineColor(AliTRDCalPID::GetPartColor(is));
93     h->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
94     h->SetMarkerStyle(24);
95   }
96   fContainer->Add(h = new TProfile("h", "", nbins, xbins));
97   h->SetMarkerStyle(7);
98   PostData(1, fContainer);
99
100
101 //____________________________________________________________________
102 void AliTRDefficiency::UserExec(Option_t *)
103 {
104   //
105   // Do it
106   //
107
108   Int_t labelsacc[10000]; 
109   memset(labelsacc, 0, sizeof(Int_t) * 10000);
110         
111   fTracks = dynamic_cast<TObjArray *>(GetInputData(1));
112   if(!fTracks) return;
113   if(!fTracks->GetEntriesFast()) return;
114   else AliDebug(2, Form("Tracks[%d] for %s", fTracks->GetEntriesFast(), GetName()));
115   if(!fMissed){ 
116     fMissed = new TClonesArray("AliTRDtrackInfo", 10);
117     fMissed->SetOwner();
118   }
119
120   Float_t mom;
121   Int_t selection[10000], nselect = 0;
122   ULong_t status; Int_t pidx;
123   Int_t nTRD = 0, nTPC = 0, nMiss = 0;
124   AliTRDtrackInfo     *track = NULL;
125   AliTrackReference     *ref = NULL;
126   AliExternalTrackParam *esd = NULL;
127   for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
128     track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
129
130                 if(!track->HasESDtrack()) continue;
131     status = track->GetStatus();
132
133     // missing TPC propagation - interesting for SA
134     if(!(status&AliESDtrack::kTPCout)) continue;
135
136     // missing MC info.
137     if(HasMCdata() && track->GetNTrackRefs() <= 1) continue;
138    
139     nTPC++;
140     selection[nselect++]=itrk;
141     ref  = track->GetTrackRef(0);
142     esd  = track->GetESDinfo()->GetOuterParam();
143     mom  = ref ? ref->P(): esd->P();
144     pidx = AliTRDCalPID::GetPartIndex(track->GetPDG());
145     pidx = TMath::Max(pidx, 0);
146     AliDebug(4, Form("PID: %d", pidx));
147
148     //Int_t n = track->GetNumberOfClusters(); 
149     // where are this tracklets ???
150     //if(ncls0 > ncls1) printf("%3d ESD[%3d] TRD[%3d|%3d]\n", itrk, ncls0, ncls1, n);
151     if(track->GetNumberOfClustersRefit()){ 
152       ((TProfile*)fContainer->At(pidx))->Fill(mom, 1.);
153                         labelsacc[nTRD] = track->GetLabel();
154       nTRD++;
155       continue;
156     }
157
158
159
160     Float_t xmed, xleng;
161     Int_t iref = 1; Bool_t found = kFALSE;
162     while((ref = track->GetTrackRef(iref))){
163       xmed = .5*(ref->LocalX() + track->GetTrackRef(iref-1)->LocalX());
164       xleng= (ref->LocalX() - track->GetTrackRef(iref-1)->LocalX());
165       if(TMath::Abs(xmed - 298.5) < .5 &&
166         TMath::Abs(xleng - 3.7) < .1){ 
167         found = kTRUE;
168         break;
169       }
170       iref++;
171     }
172     if(!found){ 
173       nTPC--;
174       // track missing first layer. Maybe interesting for SA.
175       continue;
176     }
177     nselect--;
178     new ((*fMissed)[nMiss]) AliTRDtrackInfo(*track);
179     nMiss++;
180   }
181   AliDebug(2, Form("%3d Tracks: ESD[%3d] TPC[%3d] TRD[%3d | %5.2f%%] Off[%d]", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), fTracks->GetEntriesFast(), nTPC, nTRD, nTPC ? 1.E2*nTRD/float(nTPC) : 0., fMissed->GetEntriesFast()));
182
183
184   // Find double tracks
185   Float_t threshold = 10.;
186   AliTrackReference *refMiss = NULL;
187   AliExternalTrackParam *op = NULL;
188   AliTRDtrackInfo       *tt = NULL;
189   for(Int_t imiss=0; imiss<nMiss; imiss++){
190     //printf("Searching missing %d ...\n", imiss);
191
192     // get outer param of missed
193     tt = (AliTRDtrackInfo*)fMissed->UncheckedAt(imiss);
194     op = tt->GetESDinfo()->GetOuterParam();
195     Double_t alpha = op->GetAlpha(), cosa = TMath::Cos(alpha), sina = TMath::Sin(alpha);
196
197     Double_t xyz[3], x0, y0, z0, x, y, z, dx, dy, dz, d;
198
199     Bool_t bFOUND = kFALSE;
200     for(Int_t iselect=0; iselect<nselect; iselect++){
201       track = (AliTRDtrackInfo*)fTracks->UncheckedAt(selection[iselect]);
202
203       // check first MC ... if available
204       d = 0;
205       for(Int_t iref=0; iref<track->GetNTrackRefs(); iref++){
206         if(!(ref = track->GetTrackRef(iref))) continue;
207         if((refMiss = tt->GetTrackRef(iref))){
208           dy = ref->LocalY() - refMiss->LocalY();
209           dz = ref->Z() - refMiss->Z();
210         } else {
211           // compare missOP with refTrackRef in LTC
212           x0 = ref->LocalX();
213           op->GetYAt(x0, AliTracker::GetBz(), y0);
214           op->GetZAt(x0, AliTracker::GetBz(), z0);
215           dy = y0 - ref->LocalY();
216           dz = z0 - ref->Z();
217         }
218         d += (dy*dy + dz*dz);
219       }
220       //printf("\td[%d] = %f N[%d]\n", selection[iselect], d, track->GetNTrackRefs());
221       if((track->GetNTrackRefs())){ 
222         d /= track->GetNTrackRefs();
223         if(d < threshold){
224           //printf("\t\tFound %2d in ref[%3d] : d[%f]\n", imiss, selection[iselect], d/track->GetNTrackRefs());
225           bFOUND = kTRUE; break;
226         }
227       }
228
229       // process outer param ... always available
230       // compare missOP with OP in GTC
231       esd = track->GetESDinfo()->GetOuterParam();
232       esd->GetXYZ(xyz);
233       x0 = esd->GetX();
234       op->GetYAt(x0, AliTracker::GetBz(), y0);
235       op->GetZAt(x0, AliTracker::GetBz(), z0);
236       x = x0*cosa - y0*sina;
237       y = x0*sina + y0*cosa;
238       z = z0;
239       dx=xyz[0]-x;
240       dy=xyz[1]-y;
241       dz=xyz[2]-z;
242       d = dx*dx+dy*dy+dz*dz;
243       //printf("\td[%d] = %f op\n", selection[iselect], d);
244       if(d < threshold){
245         //printf("\t\tFound %2d in op[%3d]  : d[%f] dx[%5.2f] dy[%5.2f] dz[%5.2f]\n", imiss, selection[iselect], d, dx, dy, dz);
246         bFOUND = kTRUE; break;
247       }
248     }
249     if(bFOUND) nTPC--;
250     else{ 
251       ref = tt->GetTrackRef(0);
252       mom = ref ? ref->P(): op->P();
253       pidx = AliTRDCalPID::GetPartIndex(tt->GetPDG());
254       pidx = TMath::Max(pidx, 0);
255       ((TProfile*)fContainer->At(pidx))->Fill(mom, 0.);
256       AliDebug(2, Form("  NOT bFOUND Id[%d] Mom[%f]\n", tt->GetTrackId(), mom));
257     }
258   }
259
260   AliDebug(2, Form("%3d Tracks: ESD[%3d] TPC[%3d] TRD[%3d | %5.2f%%] Off[%d]", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), fTracks->GetEntriesFast(), nTPC, nTRD, nTPC ? 1.E2*nTRD/float(nTPC) : 0., fMissed->GetEntriesFast()));
261
262   //fMissed->Delete();
263         // check for double countings
264         Int_t indices[10000]; memset(indices, 0, sizeof(Int_t) * 10000);
265         TMath::Sort(nTRD, labelsacc, indices);
266         if(DebugLevel() > 2){
267         for(Int_t itk = 0; itk < nTRD - 1; itk++)
268                 if(labelsacc[indices[itk]] ==labelsacc[indices[itk + 1]]) printf("Double counted MC track: %d\n", labelsacc[indices[itk]]);
269         }
270 }
271
272
273 //____________________________________________________________________
274 Bool_t AliTRDefficiency::GetRefFigure(Int_t ifig)
275 {
276 // Steer reference figures
277
278   if(!gPad){
279     AliWarning("Please provide a canvas to draw results.");
280     return kFALSE;
281   }
282   gPad->SetLogx();
283
284   TLegend *leg(NULL);
285   Bool_t bFIRST(kTRUE);
286   TProfile *h(NULL);
287   switch(ifig){
288   case 0:
289     h = (TProfile*)fContainer->At(AliPID::kSPECIES);
290     for(Int_t is=0; is<AliPID::kSPECIES; is++){
291       h->Add((TProfile*)fContainer->At(is));
292     }
293     h->SetMarkerStyle(24);
294     h->SetMarkerColor(kBlack);
295     h->SetLineColor(kBlack);
296     h->SetTitle("TRD Efficiency integrated");
297     h->SetXTitle("p [GeV/c]");
298     h->GetXaxis()->SetMoreLogLabels();
299     h->SetYTitle("Efficiency");
300     h->GetYaxis()->CenterTitle();
301     h->Draw("e1");
302     break;
303   case 1:
304     bFIRST = kTRUE;
305     for(Int_t is=0; is<AliPID::kSPECIES; is++){
306       if(!(h = (TProfile*)fContainer->At(is))) continue;
307       h->SetMarkerStyle(24);
308       if(bFIRST){
309         h->Draw("e1");
310         h->SetXTitle("p [GeV/c]");
311         h->GetXaxis()->SetMoreLogLabels();
312         h->SetYTitle("Efficiency");
313         h->GetYaxis()->CenterTitle();
314         h->GetYaxis()->SetRangeUser(0.8, 1.05);
315         leg=new TLegend(.7, .2, .98, .6);
316         leg->SetHeader("Species");
317         leg->SetBorderSize(0);
318         leg->SetFillStyle(0);
319         leg->AddEntry(h, h->GetTitle(), "pl");
320       } else {
321         leg->AddEntry(h, h->GetTitle(), "pl");
322         h->Draw("same e1");
323       }
324       bFIRST = kFALSE;
325     }
326     if(leg) leg->Draw();
327     break;
328   }
329   return kTRUE;
330 }
331
332
333 //____________________________________________________________________
334 Bool_t AliTRDefficiency::PostProcess()
335 {
336   fNRefFigures = HasMCdata() ? 2 : 1; 
337   return kTRUE;
338 }