TRD/qaRec/AliTRDtrackInfoGen.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: AliTRDtrackInfoGen.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 <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"
  46 #include "AliESDHeader.h"
  47 #include "AliESDtrack.h"
  48 #include "AliMCParticle.h"
  49 #include "AliPID.h"
  50 #include "AliStack.h"
  51 #include "AliTRDtrackV1.h"
  52 #include "AliTrackReference.h"
  53 #include "AliTRDgeometry.h"
  54 #include "TTreeStream.h"
  55
  56 #include <cstdio>
  57 #include <cstring>
  58 #include <iostream>
  59
  60 #include "AliTRDtrackInfoGen.h"
  61 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
  62 #include "AliTRDtrackInfo/AliTRDeventInfo.h"
  63
  64 ClassImp(AliTRDtrackInfoGen)
  65
  66 const Float_t AliTRDtrackInfoGen::xTPC = 290.;
  67 const Float_t AliTRDtrackInfoGen::xTOF = 365.;
  68
  69 //____________________________________________________________________
  70 AliTRDtrackInfoGen::AliTRDtrackInfoGen():
  71   AliTRDrecoTask("InfoGen", "Track List Generator")
  72   ,fESD(0x0)
  73   ,fMC(0x0)
  74   ,fESDfriend(0x0)
  75   ,fTrackInfo(0x0)
  76   ,fEventInfo(0x0)
  77 {
  78   //
  79   // Default constructor
  80   //
  81
  82   DefineInput(0, TChain::Class());
  83   DefineOutput(0, TObjArray::Class());
  84   DefineOutput(1, AliTRDeventInfo::Class());
  85   //DefineOutput(1, TTree::Class());
  86 }
  87
  88 //____________________________________________________________________
  89 AliTRDtrackInfoGen::~AliTRDtrackInfoGen()
  90 {
  91   if(fTrackInfo) delete fTrackInfo; fTrackInfo = 0x0;
  92   if(fEventInfo) delete fEventInfo; fEventInfo = 0x0;
  93   if(fContainer){
  94     fContainer->Delete(); delete fContainer;
  95     fContainer = 0x0;
  96   }
  97 }
  98
  99 //____________________________________________________________________
 100 void AliTRDtrackInfoGen::ConnectInputData(Option_t *)
 101 {
 102   //
 103   // Link the Input Data
 104   //
 105   TTree *tree = dynamic_cast<TChain*>(GetInputData(0));
 106   if(!tree){
 107     printf("ERROR - ESD event not found");
 108   } else {
 109     tree->SetBranchStatus("Tracks", 1);
 110     tree->SetBranchStatus("ESDfriend*",1);
 111   }
 112
 113   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 114   if(!esdH){
 115     printf("ERROR - ESD input handler not found");
 116   } else {
 117     fESD = esdH->GetEvent();
 118     if(!fESD){
 119       printf("ERROR - ESD event not found");
 120     } else {
 121       esdH->SetActiveBranches("ESDfriend*");
 122       fESDfriend = (AliESDfriend *)fESD->FindListObject("AliESDfriend");
 123     }
 124   }
 125   if(HasMCdata()){
 126     AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 127     if(!mcH){
 128       AliError("MC input handler not found");
 129     } else {
 130       fMC = mcH->MCEvent();
 131     }
 132   }
 133 }
 134
 135 //____________________________________________________________________
 136 void AliTRDtrackInfoGen::CreateOutputObjects()
 137 {
 138   //
 139   // Create Output Containers (TObjectArray containing 1D histograms)
 140   //
 141   fTrackInfo = new AliTRDtrackInfo();
 142   fEventInfo = new AliTRDeventInfo();
 143   fContainer = new TObjArray(1000);
 144   fContainer->SetOwner(kTRUE);
 145
 146 /*  OpenFile(1, "RECREATE");
 147   fTree = new TTree("trd", "extract of the TRD detector");
 148   fTree->Branch("info",  &fTrackInfo);
 149   printf("output tree build in %s\n", fTree->GetDirectory()->GetName());*/
 150 }
 151
 152 //____________________________________________________________________
 153 void AliTRDtrackInfoGen::Exec(Option_t *){
 154   //
 155   // Run the Analysis
 156   //
 157   if(!fESD){
 158     puts("Error: ESD not found");
 159     return;
 160   }
 161   if(!fESDfriend){
 162     puts("Error: ESD friend not found");
 163     return;
 164   }
 165   if(HasMCdata() && !fMC){
 166     puts("Error: Monte Carlo Event not available");
 167     return;
 168   }
 169   fContainer->Delete();
 170   fEventInfo->Delete("");
 171   fESD->SetESDfriend(fESDfriend);
 172   new(fEventInfo)AliTRDeventInfo(fESD->GetHeader(), const_cast<AliESDRun *>(fESD->GetESDRun()));
 173
 174   Bool_t *trackMap = 0x0;
 175   AliStack * mStack = 0x0;
 176   if(HasMCdata()){
 177     mStack = fMC->Stack();
 178     if(!mStack){
 179       puts("Error: Cannot get the Monte Carlo Stack");
 180       return;
 181     }
 182     trackMap = new Bool_t[fMC->GetNumberOfTracks()];
 183     memset(trackMap, 0, sizeof(Bool_t) * fMC->GetNumberOfTracks());
 184   }
 185
 186   Int_t nTRD = 0, nTPC = 0, nclsTrklt;
 187   Int_t nTracks = fESD->GetNumberOfTracks();
 188   if(fDebugLevel>=1){
 189     printf("Entry[%3d] Tracks: ESD[%d] MC[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nTracks, HasMCdata() ? mStack->GetNtrack() : 0);
 190   }
 191   AliESDtrack *esdTrack = 0x0;
 192   AliESDfriendTrack *esdFriendTrack = 0x0;
 193   TObject *calObject = 0x0;
 194   AliTRDtrackV1 *track = 0x0;
 195   for(Int_t itrk = 0; itrk < nTracks; itrk++){
 196     esdTrack = fESD->GetTrack(itrk);
 197     if(fDebugLevel>=2) printf("\n%3d ITS[%d] TPC[%d] TRD[%d]\n", itrk, esdTrack->GetNcls(0), esdTrack->GetNcls(1), esdTrack->GetNcls(2));
 198     if(esdTrack->GetNcls(1)) nTPC++;
 199     if(esdTrack->GetNcls(2)) nTRD++;
 200
 201     // look at external track param
 202     const AliExternalTrackParam *op = esdTrack->GetOuterParam();
 203     Double_t xyz[3];
 204     if(op){
 205         op->GetXYZ(xyz);
 206         op->Global2LocalPosition(xyz, op->GetAlpha());
 207         if(fDebugLevel>=2) printf("op @ X[%7.3f]\n", xyz[0]);
 208     }
 209
 210     // read MC info
 211     Int_t fPdg = -1;
 212     Int_t label = -1;
 213     if(HasMCdata()){
 214       label = esdTrack->GetLabel();
 215       if(label < fMC->GetNumberOfTracks()) trackMap[TMath::Abs(label)] = kTRUE; // register the track
 216       //if (TMath::Abs(label) > mStack->GetNtrack()) continue;
 217       AliMCParticle *mcParticle = 0x0;
 218       if(!(mcParticle = fMC->GetTrack(TMath::Abs(label)))){
 219         printf("E - AliTRDtrackInfoGen::Exec() : MC particle missing for ESD label %d\n", label);
 220         continue;
 221       }
 222       fPdg = mcParticle->Particle()->GetPdgCode();
 223       Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
 224       Int_t iref = 0; AliTrackReference *ref = 0x0;
 225       while(iref<nRefs){
 226         ref = mcParticle->GetTrackReference(iref);
 227         if(ref->LocalX() > xTPC) break;
 228         //printf("\ttrackRef[%2d] @ %7.3f\n", iref, ref->LocalX());
 229         iref++;
 230       }
 231
 232       new(fTrackInfo) AliTRDtrackInfo();
 233       fTrackInfo->SetPDG(fPdg);
 234       fTrackInfo->SetPrimary(mcParticle->Particle()->IsPrimary());
 235       Int_t jref = iref;//, kref = 0;
 236       while(jref<nRefs){
 237         ref = mcParticle->GetTrackReference(jref);
 238         if(ref->LocalX() > xTOF) break;
 239         if(fDebugLevel>=3) printf("\ttrackRef[%2d (%2d)] @ %7.3f OK\n", jref-iref, jref, ref->LocalX());
 240         fTrackInfo->AddTrackRef(ref);
 241         jref++;
 242       }
 243       if(fDebugLevel>=2) printf("NtrackRefs[%d(%d)]\n", fTrackInfo->GetNTrackRefs(), nRefs);
 244     } else {
 245       new (fTrackInfo) AliTRDtrackInfo();
 246       fTrackInfo->SetPDG(fPdg);
 247     }
 248
 249     // copy some relevant info to TRD track info
 250     fTrackInfo->SetStatus(esdTrack->GetStatus());
 251     fTrackInfo->SetTrackId(esdTrack->GetID());
 252     Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
 253     fTrackInfo->SetESDpid(p);
 254     fTrackInfo->SetESDpidQuality(esdTrack->GetTRDpidQuality());
 255     fTrackInfo->SetLabel(label);
 256     fTrackInfo->SetNumberOfClustersRefit(esdTrack->GetNcls(2));
 257     // some other Informations which we may wish to store in order to find problematic cases
 258     fTrackInfo->SetKinkIndex(esdTrack->GetKinkIndex(0));
 259     fTrackInfo->SetTPCncls(static_cast<UShort_t>(esdTrack->GetNcls(1)));
 260     nclsTrklt = 0;
 261
 262
 263     // read REC info
 264     esdFriendTrack = fESDfriend->GetTrack(itrk);
 265     if(esdFriendTrack){
 266       Int_t icalib = 0;
 267       while((calObject = esdFriendTrack->GetCalibObject(icalib++))){
 268         if(strcmp(calObject->IsA()->GetName(),"AliTRDtrackV1") != 0) continue; // Look for the TRDtrack
 269         if(!(track = dynamic_cast<AliTRDtrackV1*>(calObject))) break;
 270         nTRD++;
 271         if(fDebugLevel>=3) printf("TRD track OK\n");
 272         fTrackInfo->SetTrack(track);
 273         break;
 274       }
 275       if(fDebugLevel>=2) printf("Ntracklets[%d]\n", fTrackInfo->GetNTracklets());
 276     } else if(fDebugLevel>=2) printf("No ESD friends\n");
 277     if(op) fTrackInfo->SetOuterParam(op);
 278
 279     if(fDebugLevel >= 1){
 280       AliTRDtrackInfo info(*fTrackInfo);
 281       (*fDebugStream) << "trackInfo"
 282       << "TrackInfo.=" << &info
 283       << "\n";
 284       info.Delete("");
 285     }
 286
 287     fContainer->Add(new AliTRDtrackInfo(*fTrackInfo));
 288     fTrackInfo->Delete("");
 289   }
 290   if(fDebugLevel>=2) printf("%3d Tracks: TPC[%d] TRD[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nTPC, nTRD);
 291
 292   // Insert also MC tracks which are passing TRD where the track is not reconstructed
 293   if(HasMCdata()){
 294     if(fDebugLevel > 10){
 295       printf("Output of the MC track map:\n");
 296       for(Int_t itk = 0; itk < fMC->GetNumberOfTracks();  itk++)
 297         printf("trackMap[%d] = %s\n", itk, trackMap[itk] == kTRUE ? "TRUE" : "kFALSE");
 298     }
 299
 300     for(Int_t itk = 0; itk < fMC->GetNumberOfTracks(); itk++){
 301       if(fDebugLevel >=2 ) printf("Number of MC tracks: %d\n", fMC->GetNumberOfTracks());
 302       if(trackMap[itk]) continue;
 303       AliMCParticle *mcParticle = fMC->GetTrack(TMath::Abs(itk));
 304       Int_t fPdg = mcParticle->Particle()->GetPdgCode();
 305       Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
 306       Int_t iref = 0; AliTrackReference *ref = 0x0;
 307       Int_t nRefsTRD = 0;
 308       new(fTrackInfo) AliTRDtrackInfo();
 309       fTrackInfo->SetPDG(fPdg);
 310       while(iref<nRefs){
 311         ref = mcParticle->GetTrackReference(iref);
 312         if(fDebugLevel > 3) printf("\ttrackRef[%2d] @ %7.3f", iref, ref->LocalX());
 313         if(ref->LocalX() > 250. && ref->LocalX() < 370.){
 314           if(fDebugLevel > 3) printf(" OK\n");
 315           fTrackInfo->AddTrackRef(ref);
 316           nRefsTRD++;
 317         }
 318         else
 319           if(fDebugLevel > 3) printf("\n");
 320         iref++;
 321       }
 322       if(!nRefsTRD){
 323         // In this stage we at least require 1 hit inside TRD. What will be done with this tracks is a task for the
 324         // analysis job
 325         fTrackInfo->Delete("");
 326         continue;
 327       }
 328       fTrackInfo->SetPrimary(mcParticle->Particle()->IsPrimary());
 329       fTrackInfo->SetLabel(itk);
 330       if(fDebugLevel >= 1){
 331         AliTRDtrackInfo info(*fTrackInfo);
 332         (*fDebugStream) << "trackInfo"
 333         << "TrackInfo.=" << &info
 334         << "\n";
 335         info.Delete("");
 336       }
 337       if(fDebugLevel > 2)printf("Registering rejected MC track with label %d\n", itk);
 338       fContainer->Add(new AliTRDtrackInfo(*fTrackInfo));
 339       fTrackInfo->Delete("");
 340     }
 341     delete[] trackMap;
 342   }
 343   PostData(0, fContainer);
 344   PostData(1, fEventInfo);
 345 }
 346
 347
 348 //____________________________________________________________________
 349 void AliTRDtrackInfoGen::Terminate(Option_t *)
 350 {
 351   //
 352   // Stays empty because we are only interested in the tree
 353   //
 354   if(fDebugLevel>=1) printf("Terminate:\n");
 355   //TFile *f =((TFile*)gROOT->FindObject("TRD.TrackInfo.root"));
 356   //f->cd(); f->Write(); f->Close();
 357
 358   if(fDebugStream){
 359     delete fDebugStream;
 360     fDebugStream = 0x0;
 361   }
 362 }