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