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