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