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