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