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