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