TRD/qaRec/AliTRDtrackingEfficiency.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: AliTRDtrackingEfficiency.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 "AliESDtrack.h"
  33 #include "AliTrackReference.h"
  34 #include "AliExternalTrackParam.h"
  35 #include "AliTracker.h"
  36 #include "AliMagFMaps.h"
  37 #include "AliAnalysisManager.h"
  38
  39 #include "AliTRDtrackingEfficiency.h"
  40 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
  41
  42 ClassImp(AliTRDtrackingEfficiency)
  43
  44 //____________________________________________________________________
  45 AliTRDtrackingEfficiency::AliTRDtrackingEfficiency()
  46   :AliTRDrecoTask("TrackingEff", "Barrel Tracking Effiency")
  47   ,fMissed(0x0)
  48 {
  49   //
  50   // Default constructor
  51   //
  52 }
  53
  54 //____________________________________________________________________
  55 AliTRDtrackingEfficiency::~AliTRDtrackingEfficiency()
  56 {
  57   if(fMissed){
  58     fMissed->Delete();
  59     delete fMissed;
  60   }
  61 }
  62
  63 //____________________________________________________________________
  64 void  AliTRDtrackingEfficiency::CreateOutputObjects()
  65 {
  66   //
  67   // Create output objects
  68   //
  69
  70   OpenFile(0, "RECREATE");
  71   const Int_t nbins = 11;
  72   Float_t xbins[nbins+1] = {.5, .7, .9, 1.3, 1.7, 2.4, 3.5, 4.5, 5.5, 7., 9., 11.};
  73   fContainer = new TObjArray();
  74   fContainer->Add(new TProfile("h", "", nbins, xbins));
  75 }
  76
  77 //____________________________________________________________________
  78 void AliTRDtrackingEfficiency::Exec(Option_t *)
  79 {
  80   //
  81   // Do it
  82   //
  83
  84   if(!AliTracker::GetFieldMap()){
  85     AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., AliMagFMaps::k5kG);
  86     AliTracker::SetFieldMap(field, kTRUE);
  87   }
  88   TProfile *h = (TProfile*)fContainer->At(0);
  89         Int_t labelsacc[10000]; memset(labelsacc, 0, sizeof(Int_t) * 10000);
  90
  91   if(!fMissed){
  92     fMissed = new TClonesArray("AliTRDtrackInfo", 10);
  93     fMissed->SetOwner();
  94   }
  95
  96   Float_t mom;
  97   Int_t selection[10000], nselect = 0;
  98   ULong_t status;
  99   Int_t nTRD = 0, nTPC = 0, nMiss = 0;
 100   AliTRDtrackInfo     *track = 0x0;
 101   AliTrackReference     *ref = 0x0;
 102   AliExternalTrackParam *esd = 0x0;
 103   for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
 104     track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
 105                 if(!track->HasESDtrack()) continue;
 106     status = track->GetStatus();
 107     if(!(status&AliESDtrack::kTPCout)){
 108       // missing TPC propagation - interesting for SA
 109       continue;
 110     }
 111
 112     nTPC++;
 113     selection[nselect++]=itrk;
 114
 115     //Int_t n = track->GetNumberOfClusters();
 116     // where are this tracklets ???
 117     //if(ncls0 > ncls1) printf("%3d ESD[%3d] TRD[%3d|%3d]\n", itrk, ncls0, ncls1, n);
 118     if(track->GetNumberOfClustersRefit()){
 119       ref = track->GetTrackRef(0);
 120       esd = track->GetOuterParam();
 121       mom = ref ? ref->P(): esd->P();
 122       //printf("FOUND Id[%d] mom[%f]\n", itrk, mom);
 123       h->Fill(mom, 1.);
 124                         labelsacc[nTRD] = track->GetLabel();
 125       nTRD++;
 126       continue;
 127     }
 128
 129
 130     Int_t nrefs = track->GetNTrackRefs();
 131     if(nrefs<=1){
 132       nTPC--;
 133       // we don't have MC info.
 134       // we should discard this track
 135       continue;
 136     }
 137
 138     Float_t xmed, xleng;
 139     Int_t iref = 1; Bool_t found = kFALSE;
 140     while((ref = track->GetTrackRef(iref))){
 141       xmed = .5*(ref->LocalX() + track->GetTrackRef(iref-1)->LocalX());
 142       xleng= (ref->LocalX() - track->GetTrackRef(iref-1)->LocalX());
 143       if(TMath::Abs(xmed - 298.5) < .5 &&
 144         TMath::Abs(xleng - 3.7) < .1){
 145         found = kTRUE;
 146         break;
 147       }
 148       iref++;
 149     }
 150     if(!found){
 151       nTPC--;
 152       // track missing first layer. Maybe interesting for SA.
 153       continue;
 154     }
 155     nselect--;
 156     new ((*fMissed)[nMiss]) AliTRDtrackInfo(*track);
 157     nMiss++;
 158   }
 159   if(fDebugLevel>=1) printf("%3d Tracks: ESD[%3d] TPC[%3d] TRD[%3d | %5.2f%%] Off[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), fTracks->GetEntriesFast(), nTPC, nTRD, nTPC ? 1.E2*nTRD/float(nTPC) : 0., fMissed->GetEntriesFast());
 160
 161
 162   // Find double tracks
 163   Float_t threshold = 10.;
 164   AliTrackReference *refMiss = 0x0;
 165   AliExternalTrackParam *op = 0x0;
 166   AliTRDtrackInfo       *tt = 0x0;
 167   for(Int_t imiss=0; imiss<nMiss; imiss++){
 168     //printf("Searching missing %d ...\n", imiss);
 169
 170     // get outer param of missed
 171     tt = (AliTRDtrackInfo*)fMissed->UncheckedAt(imiss);
 172     op = tt->GetOuterParam();
 173     Double_t alpha = op->GetAlpha(), cosa = TMath::Cos(alpha), sina = TMath::Sin(alpha);
 174
 175     Double_t xyz[3], x0, y0, z0, x, y, z, dx, dy, dz, d;
 176
 177     Bool_t FOUND = kFALSE;
 178     for(Int_t iselect=0; iselect<nselect; iselect++){
 179       track = (AliTRDtrackInfo*)fTracks->UncheckedAt(selection[iselect]);
 180
 181       // check first MC ... if available
 182       d = 0;
 183       for(Int_t iref=0; iref<track->GetNTrackRefs(); iref++){
 184         if(!(ref = track->GetTrackRef(iref))) continue;
 185         if((refMiss = tt->GetTrackRef(iref))){
 186           dy = ref->LocalY() - refMiss->LocalY();
 187           dz = ref->Z() - refMiss->Z();
 188         } else {
 189           // compare missOP with refTrackRef in LTC
 190           x0 = ref->LocalX();
 191           op->GetYAt(x0, AliTracker::GetBz(), y0);
 192           op->GetZAt(x0, AliTracker::GetBz(), z0);
 193           dy = y0 - ref->LocalY();
 194           dz = z0 - ref->Z();
 195         }
 196         d += (dy*dy + dz*dz);
 197       }
 198       //printf("\td[%d] = %f N[%d]\n", selection[iselect], d, track->GetNTrackRefs());
 199       if((track->GetNTrackRefs())){
 200         d /= track->GetNTrackRefs();
 201         if(d < threshold){
 202           //printf("\t\tFound %2d in ref[%3d] : d[%f]\n", imiss, selection[iselect], d/track->GetNTrackRefs());
 203           FOUND = kTRUE; break;
 204         }
 205       }
 206
 207       // process outer param ... always available
 208       // compare missOP with OP in GTC
 209       esd = track->GetOuterParam();
 210       esd->GetXYZ(xyz);
 211       x0 = esd->GetX();
 212       op->GetYAt(x0, AliTracker::GetBz(), y0);
 213       op->GetZAt(x0, AliTracker::GetBz(), z0);
 214       x = x0*cosa - y0*sina;
 215       y = x0*sina + y0*cosa;
 216       z = z0;
 217       dx=xyz[0]-x;
 218       dy=xyz[1]-y;
 219       dz=xyz[2]-z;
 220       d = dx*dx+dy*dy+dz*dz;
 221       //printf("\td[%d] = %f op\n", selection[iselect], d);
 222       if(d < threshold){
 223         //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);
 224         FOUND = kTRUE; break;
 225       }
 226     }
 227     if(FOUND) nTPC--;
 228     else{
 229       ref = tt->GetTrackRef(0);
 230       mom = ref ? ref->P(): op->P();
 231       h->Fill(mom, 0.);
 232       printf("\tNOT FOUND Id[%d] Mom[%f]\n", tt->GetTrackId(), mom);
 233     }
 234   }
 235
 236   //if(fDebugLevel>=1)
 237   printf("%3d Tracks: ESD[%3d] TPC[%3d] TRD[%3d | %5.2f%%] Off[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), fTracks->GetEntriesFast(), nTPC, nTRD, nTPC ? 1.E2*nTRD/float(nTPC) : 0., fMissed->GetEntriesFast());
 238
 239   //fMissed->Delete();
 240         // check for double countings
 241         Int_t indices[10000]; memset(indices, 0, sizeof(Int_t) * 10000);
 242         TMath::Sort(nTRD, labelsacc, indices);
 243         if(fDebugLevel > 2){
 244         for(Int_t itk = 0; itk < nTRD - 1; itk++)
 245                 if(labelsacc[indices[itk]] ==labelsacc[indices[itk + 1]]) printf("Double counted MC track: %d\n", labelsacc[indices[itk]]);
 246         }
 247   PostData(0, fContainer);
 248 }
 249
 250 //____________________________________________________________________
 251 void AliTRDtrackingEfficiency::Terminate(Option_t *)
 252 {
 253   //
 254   // Terminate
 255   //
 256
 257   if(fDebugStream){
 258     delete fDebugStream;
 259     fDebugStream = 0x0;
 260     fDebugLevel = 0;
 261   }
 262
 263   fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
 264   if (!fContainer) {
 265     Printf("ERROR: list not available");
 266     return;
 267   }
 268
 269 }
 270