]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/qaRec/AliTRDinfoGen.cxx
restore functionality of "mc" and "friends" flags
[u/mrichter/AliRoot.git] / TRD / qaRec / AliTRDinfoGen.cxx
CommitLineData
814ecea4 1/**************************************************************************
d2381af5 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**************************************************************************/
814ecea4 15
873458ab 16/* $Id: AliTRDinfoGen.cxx 27496 2008-07-22 08:35:45Z cblume $ */
814ecea4 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 <TObject.h>
30#include <TH1F.h>
31#include <TFile.h>
32#include <TTree.h>
33#include <TROOT.h>
34#include <TChain.h>
35#include <TParticle.h>
36
37#include "AliLog.h"
38#include "AliAnalysisManager.h"
39#include "AliESDEvent.h"
40#include "AliMCEvent.h"
41#include "AliESDInputHandler.h"
42#include "AliMCEventHandler.h"
43
44#include "AliESDfriend.h"
45#include "AliESDfriendTrack.h"
a24151d1 46#include "AliESDHeader.h"
814ecea4 47#include "AliESDtrack.h"
48#include "AliMCParticle.h"
c46a7947 49#include "AliPID.h"
814ecea4 50#include "AliStack.h"
51#include "AliTRDtrackV1.h"
52#include "AliTrackReference.h"
53#include "AliTRDgeometry.h"
5c6b33ff 54#include "AliTRDcluster.h"
55#include "AliTRDseedV1.h"
814ecea4 56#include "TTreeStream.h"
57
58#include <cstdio>
d068cde5 59#include <climits>
814ecea4 60#include <cstring>
afefec95 61#include <iostream>
814ecea4 62
873458ab 63#include "AliTRDinfoGen.h"
64#include "info/AliTRDtrackInfo.h"
65#include "info/AliTRDeventInfo.h"
98233fc0 66#include "info/AliTRDv0Info.h"
814ecea4 67
873458ab 68ClassImp(AliTRDinfoGen)
814ecea4 69
873458ab 70const Float_t AliTRDinfoGen::xTPC = 290.;
71const Float_t AliTRDinfoGen::xTOF = 365.;
814ecea4 72
73//____________________________________________________________________
873458ab 74AliTRDinfoGen::AliTRDinfoGen():
6da3eee3 75 AliTRDrecoTask("infoGen", "MC-REC TRD-track list generator")
814ecea4 76 ,fESD(0x0)
77 ,fMC(0x0)
78 ,fESDfriend(0x0)
79 ,fTrackInfo(0x0)
a24151d1 80 ,fEventInfo(0x0)
98233fc0 81 ,fV0container(0x0)
82 ,fV0Info(0x0)
814ecea4 83{
84 //
85 // Default constructor
86 //
87
88 DefineInput(0, TChain::Class());
89 DefineOutput(0, TObjArray::Class());
a24151d1 90 DefineOutput(1, AliTRDeventInfo::Class());
98233fc0 91 DefineOutput(2, TObjArray::Class());
814ecea4 92}
93
94//____________________________________________________________________
873458ab 95AliTRDinfoGen::~AliTRDinfoGen()
ed383798 96{
94f34623 97 if(fTrackInfo) delete fTrackInfo; fTrackInfo = 0x0;
98 if(fEventInfo) delete fEventInfo; fEventInfo = 0x0;
98233fc0 99 if(fV0Info) delete fV0Info; fV0Info = 0x0;
94f34623 100 if(fContainer){
101 fContainer->Delete(); delete fContainer;
102 fContainer = 0x0;
103 }
98233fc0 104 if(fV0container){
105 fV0container->Delete(); delete fV0container;
106 fV0container = 0x0;
107 }
ed383798 108}
109
110//____________________________________________________________________
873458ab 111void AliTRDinfoGen::ConnectInputData(Option_t *)
814ecea4 112{
113 //
114 // Link the Input Data
115 //
d2381af5 116 TTree *tree = dynamic_cast<TChain*>(GetInputData(0));
117 if(!tree){
814ecea4 118 printf("ERROR - ESD event not found");
d2381af5 119 } else {
120 tree->SetBranchStatus("Tracks", 1);
121 tree->SetBranchStatus("ESDfriend*",1);
814ecea4 122 }
123
124 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
125 if(!esdH){
d2381af5 126 printf("ERROR - ESD input handler not found");
814ecea4 127 } else {
128 fESD = esdH->GetEvent();
129 if(!fESD){
d2381af5 130 printf("ERROR - ESD event not found");
814ecea4 131 } else {
132 esdH->SetActiveBranches("ESDfriend*");
133 fESDfriend = (AliESDfriend *)fESD->FindListObject("AliESDfriend");
814ecea4 134 }
135 }
4b8f8a35 136 if(HasMCdata()){
d2381af5 137 AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
138 if(!mcH){
139 AliError("MC input handler not found");
140 } else {
141 fMC = mcH->MCEvent();
142 }
814ecea4 143 }
144}
145
146//____________________________________________________________________
873458ab 147void AliTRDinfoGen::CreateOutputObjects()
814ecea4 148{
149 //
150 // Create Output Containers (TObjectArray containing 1D histograms)
151 //
152 fTrackInfo = new AliTRDtrackInfo();
a24151d1 153 fEventInfo = new AliTRDeventInfo();
98233fc0 154 fV0Info = new AliTRDv0Info();
3d86166d 155 fContainer = new TObjArray(1000);
94f34623 156 fContainer->SetOwner(kTRUE);
98233fc0 157 fV0container = new TObjArray(50);
158 fV0container->SetOwner(kTRUE);
814ecea4 159
814ecea4 160}
161
162//____________________________________________________________________
873458ab 163void AliTRDinfoGen::Exec(Option_t *){
d2381af5 164 //
165 // Run the Analysis
166 //
167 if(!fESD){
b94feda9 168 AliError("Failed retrieving ESD event");
d2381af5 169 return;
170 }
171 if(!fESDfriend){
b94feda9 172 AliError("Failed retrieving ESD friend event");
d2381af5 173 return;
174 }
175 if(HasMCdata() && !fMC){
b94feda9 176 AliError("Failed retrieving MC event");
d2381af5 177 return;
178 }
3d86166d 179 fContainer->Delete();
98233fc0 180 fV0container->Delete();
afefec95 181 fEventInfo->Delete("");
d2381af5 182 fESD->SetESDfriend(fESDfriend);
a24151d1 183 new(fEventInfo)AliTRDeventInfo(fESD->GetHeader(), const_cast<AliESDRun *>(fESD->GetESDRun()));
d2381af5 184
185 Bool_t *trackMap = 0x0;
186 AliStack * mStack = 0x0;
2d65e748 187 Int_t nTracksMC = HasMCdata() ? fMC->GetNumberOfTracks() : 0, nTracksESD = fESD->GetNumberOfTracks();
d2381af5 188 if(HasMCdata()){
189 mStack = fMC->Stack();
190 if(!mStack){
b94feda9 191 AliError("Failed retrieving MC Stack");
d2381af5 192 return;
193 }
2d65e748 194 trackMap = new Bool_t[nTracksMC];
195 memset(trackMap, 0, sizeof(Bool_t) * nTracksMC);
d2381af5 196 }
197
198 Int_t nTRD = 0, nTPC = 0, nclsTrklt;
d2381af5 199 if(fDebugLevel>=1){
2d65e748 200 printf("Entry[%3d] Tracks: ESD[%d] MC[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nTracksESD, nTracksMC);
d2381af5 201 }
202 AliESDtrack *esdTrack = 0x0;
203 AliESDfriendTrack *esdFriendTrack = 0x0;
204 TObject *calObject = 0x0;
205 AliTRDtrackV1 *track = 0x0;
5c6b33ff 206 AliTRDseedV1 *tracklet = 0x0;
207 AliTRDcluster *cl = 0x0;
2d65e748 208 for(Int_t itrk = 0; itrk < nTracksESD; itrk++){
d2381af5 209 esdTrack = fESD->GetTrack(itrk);
210 if(fDebugLevel>=2) printf("\n%3d ITS[%d] TPC[%d] TRD[%d]\n", itrk, esdTrack->GetNcls(0), esdTrack->GetNcls(1), esdTrack->GetNcls(2));
211 if(esdTrack->GetNcls(1)) nTPC++;
212 if(esdTrack->GetNcls(2)) nTRD++;
814ecea4 213
303f052c 214 // look at external track param
d2381af5 215 const AliExternalTrackParam *op = esdTrack->GetOuterParam();
216 Double_t xyz[3];
217 if(op){
218 op->GetXYZ(xyz);
219 op->Global2LocalPosition(xyz, op->GetAlpha());
220 if(fDebugLevel>=2) printf("op @ X[%7.3f]\n", xyz[0]);
221 }
814ecea4 222
d2381af5 223 // read MC info
224 Int_t fPdg = -1;
d068cde5 225 Int_t label = -1; UInt_t alab=UINT_MAX;
d2381af5 226 if(HasMCdata()){
b94feda9 227 label = esdTrack->GetLabel();
228 alab = TMath::Abs(label);
229 // register the track
230 if(alab < UInt_t(nTracksMC)){
231 trackMap[alab] = kTRUE;
232 } else {
82b33eb2 233 AliError(Form("MC label[%d] outside scope for Ev[%d] Trk[%d].", label, (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), itrk));
b94feda9 234 continue;
235 }
5d6dc395 236 AliMCParticle *mcParticle = 0x0;
7aad0c47 237 if(!(mcParticle = (AliMCParticle*) fMC->GetTrack(alab))){
82b33eb2 238 AliError(Form("MC particle label[%d] missing for Ev[%d] Trk[%d].", label, (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), itrk));
5d6dc395 239 continue;
240 }
d2381af5 241 fPdg = mcParticle->Particle()->GetPdgCode();
242 Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
243 Int_t iref = 0; AliTrackReference *ref = 0x0;
244 while(iref<nRefs){
245 ref = mcParticle->GetTrackReference(iref);
303f052c 246 if(ref->LocalX() > xTPC) break;
d2381af5 247 iref++;
248 }
814ecea4 249
93e41bce 250 new(fTrackInfo) AliTRDtrackInfo();
251 fTrackInfo->SetPDG(fPdg);
d2381af5 252 fTrackInfo->SetPrimary(mcParticle->Particle()->IsPrimary());
253 Int_t jref = iref;//, kref = 0;
254 while(jref<nRefs){
255 ref = mcParticle->GetTrackReference(jref);
303f052c 256 if(ref->LocalX() > xTOF) break;
d2381af5 257 if(fDebugLevel>=3) printf("\ttrackRef[%2d (%2d)] @ %7.3f OK\n", jref-iref, jref, ref->LocalX());
258 fTrackInfo->AddTrackRef(ref);
259 jref++;
260 }
d2381af5 261 if(fDebugLevel>=2) printf("NtrackRefs[%d(%d)]\n", fTrackInfo->GetNTrackRefs(), nRefs);
262 } else {
93e41bce 263 new (fTrackInfo) AliTRDtrackInfo();
264 fTrackInfo->SetPDG(fPdg);
d2381af5 265 }
814ecea4 266
d2381af5 267 // copy some relevant info to TRD track info
268 fTrackInfo->SetStatus(esdTrack->GetStatus());
269 fTrackInfo->SetTrackId(esdTrack->GetID());
c46a7947 270 Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
271 fTrackInfo->SetESDpid(p);
ed15ef4f 272 fTrackInfo->SetESDpidQuality(esdTrack->GetTRDntrackletsPID());
d2381af5 273 fTrackInfo->SetLabel(label);
274 fTrackInfo->SetNumberOfClustersRefit(esdTrack->GetNcls(2));
52a79f8d 275 // some other Informations which we may wish to store in order to find problematic cases
276 fTrackInfo->SetKinkIndex(esdTrack->GetKinkIndex(0));
277 fTrackInfo->SetTPCncls(static_cast<UShort_t>(esdTrack->GetNcls(1)));
d2381af5 278 nclsTrklt = 0;
279
814ecea4 280
d2381af5 281 // read REC info
282 esdFriendTrack = fESDfriend->GetTrack(itrk);
283 if(esdFriendTrack){
284 Int_t icalib = 0;
285 while((calObject = esdFriendTrack->GetCalibObject(icalib++))){
286 if(strcmp(calObject->IsA()->GetName(),"AliTRDtrackV1") != 0) continue; // Look for the TRDtrack
b718144c 287 if(!(track = dynamic_cast<AliTRDtrackV1*>(calObject))) break;
288 nTRD++;
d2381af5 289 if(fDebugLevel>=3) printf("TRD track OK\n");
5c6b33ff 290 // Set the clusters to unused
291 for(Int_t ipl = 0; ipl < AliTRDgeometry::kNlayer; ipl++){
292 if(!(tracklet = track->GetTracklet(ipl))) continue;
293 tracklet->ResetClusterIter();
294 while((cl = tracklet->NextCluster())) cl->Use(0);
295 }
c46a7947 296 fTrackInfo->SetTrack(track);
d2381af5 297 break;
298 }
299 if(fDebugLevel>=2) printf("Ntracklets[%d]\n", fTrackInfo->GetNTracklets());
300 } else if(fDebugLevel>=2) printf("No ESD friends\n");
301 if(op) fTrackInfo->SetOuterParam(op);
814ecea4 302
b718144c 303 if(fDebugLevel >= 1){
304 AliTRDtrackInfo info(*fTrackInfo);
d2381af5 305 (*fDebugStream) << "trackInfo"
b718144c 306 << "TrackInfo.=" << &info
d2381af5 307 << "\n";
b718144c 308 info.Delete("");
d2381af5 309 }
310
3d86166d 311 fContainer->Add(new AliTRDtrackInfo(*fTrackInfo));
d2381af5 312 fTrackInfo->Delete("");
313 }
b718144c 314 if(fDebugLevel>=2) printf("%3d Tracks: TPC[%d] TRD[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nTPC, nTRD);
4b8f8a35 315
98233fc0 316 AliESDv0 *v0 = 0x0;
317 for(Int_t iv0=0; iv0<fESD->GetNumberOfV0s(); iv0++){
318 if(!(v0 = fESD->GetV0(iv0))) continue;
319 fV0container->Add(new AliTRDv0Info(v0));
320 }
321
d2381af5 322 // Insert also MC tracks which are passing TRD where the track is not reconstructed
323 if(HasMCdata()){
324 if(fDebugLevel > 10){
325 printf("Output of the MC track map:\n");
2d65e748 326 for(Int_t itk = 0; itk < nTracksMC; itk++)
d2381af5 327 printf("trackMap[%d] = %s\n", itk, trackMap[itk] == kTRUE ? "TRUE" : "kFALSE");
328 }
329
2d65e748 330 for(Int_t itk = 0; itk < nTracksMC; itk++){
d2381af5 331 if(trackMap[itk]) continue;
7aad0c47 332 AliMCParticle *mcParticle = (AliMCParticle*) fMC->GetTrack(TMath::Abs(itk));
d2381af5 333 Int_t fPdg = mcParticle->Particle()->GetPdgCode();
334 Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
335 Int_t iref = 0; AliTrackReference *ref = 0x0;
336 Int_t nRefsTRD = 0;
93e41bce 337 new(fTrackInfo) AliTRDtrackInfo();
338 fTrackInfo->SetPDG(fPdg);
d2381af5 339 while(iref<nRefs){
340 ref = mcParticle->GetTrackReference(iref);
341 if(fDebugLevel > 3) printf("\ttrackRef[%2d] @ %7.3f", iref, ref->LocalX());
342 if(ref->LocalX() > 250. && ref->LocalX() < 370.){
343 if(fDebugLevel > 3) printf(" OK\n");
344 fTrackInfo->AddTrackRef(ref);
345 nRefsTRD++;
346 }
347 else
348 if(fDebugLevel > 3) printf("\n");
349 iref++;
350 }
351 if(!nRefsTRD){
352 // In this stage we at least require 1 hit inside TRD. What will be done with this tracks is a task for the
353 // analysis job
354 fTrackInfo->Delete("");
355 continue;
356 }
357 fTrackInfo->SetPrimary(mcParticle->Particle()->IsPrimary());
358 fTrackInfo->SetLabel(itk);
359 if(fDebugLevel >= 1){
b718144c 360 AliTRDtrackInfo info(*fTrackInfo);
d2381af5 361 (*fDebugStream) << "trackInfo"
b718144c 362 << "TrackInfo.=" << &info
d2381af5 363 << "\n";
b718144c 364 info.Delete("");
d2381af5 365 }
366 if(fDebugLevel > 2)printf("Registering rejected MC track with label %d\n", itk);
3d86166d 367 fContainer->Add(new AliTRDtrackInfo(*fTrackInfo));
d2381af5 368 fTrackInfo->Delete("");
369 }
370 delete[] trackMap;
371 }
3d86166d 372 PostData(0, fContainer);
a24151d1 373 PostData(1, fEventInfo);
98233fc0 374 PostData(2, fV0container);
814ecea4 375}
376
814ecea4 377
378//____________________________________________________________________
873458ab 379void AliTRDinfoGen::Terminate(Option_t *)
814ecea4 380{
381 //
382 // Stays empty because we are only interested in the tree
383 //
3d86166d 384 if(fDebugStream){
385 delete fDebugStream;
386 fDebugStream = 0x0;
387 }
814ecea4 388}