TRD/qaRec/AliTRDtrackingEfficiencyCombined.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: AliTRDtrackingEfficiencyCombined.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 <TObjArray.h>
  28 #include <TProfile.h>
  29 #include <TMath.h>
  30 #include <TCanvas.h>
  31 #include "TTreeStream.h"
  32
  33 #include "AliMagFMaps.h"
  34 #include "AliTracker.h"
  35 #include "AliTrackReference.h"
  36 #include "AliAnalysisManager.h"
  37
  38 #include "AliTRDseedV1.h"
  39 #include "AliTRDtrackV1.h"
  40 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
  41 #include "AliTRDtrackingEfficiencyCombined.h"
  42
  43 ClassImp(AliTRDtrackingEfficiencyCombined)
  44
  45 //_____________________________________________________________________________
  46 AliTRDtrackingEfficiencyCombined::AliTRDtrackingEfficiencyCombined()
  47   :AliTRDrecoTask("TrackingEffMC", "Combined Tracking Efficiency")
  48 {
  49   //
  50   // Default constructor
  51   //
  52 }
  53
  54
  55 //_____________________________________________________________________________
  56 void AliTRDtrackingEfficiencyCombined::CreateOutputObjects(){
  57   //
  58   // Create output objects
  59   //
  60
  61   OpenFile(0, "RECREATE");
  62
  63   const Int_t nbins = 11;
  64   Float_t xbins[nbins+1] = {.5, .7, .9, 1.3, 1.7, 2.4, 3.5, 4.5, 5.5, 7., 9., 11.};
  65
  66   fContainer = new TObjArray();
  67   fContainer->Add(new TProfile("trEffComb", "Combined Tracking Efficiency", nbins, xbins));
  68   fContainer->Add(new TProfile("trContComb", "Combined Tracking Contamination", nbins, xbins));
  69 }
  70
  71 //_____________________________________________________________________________
  72 void AliTRDtrackingEfficiencyCombined::Exec(Option_t *){
  73   //
  74   // Do it
  75   //
  76
  77         const Float_t kAlpha = 0.349065850;
  78         Int_t naccepted = 0, nrejected = 0, ndoublecounted = 0;
  79         Int_t labelsacc[10000];
  80         Int_t labelsrej[10000];
  81         Float_t momacc[10000];
  82         Float_t momrej[10000];
  83         if(!AliTracker::GetFieldMap()){
  84                 AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., AliMagFMaps::k5kG);
  85                 AliTracker::SetFieldMap(field, kTRUE);
  86         }
  87         TProfile *efficiency = (TProfile *)fContainer->At(0);
  88         TProfile *contamination = (TProfile *)fContainer->At(1);
  89
  90         Int_t nTrackInfos = fTracks->GetEntriesFast();
  91         Double_t mom = 0;
  92         AliTRDtrackV1 *TRDtrack = 0x0;
  93         AliTRDtrackInfo *trkInf = 0x0;
  94         AliTrackReference *trackRef = 0x0;
  95         for(Int_t itinf = 0; itinf < nTrackInfos; itinf++){
  96                 mom = 0.;
  97                 trkInf = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(itinf));
  98                 if(!trkInf) continue;
  99                 if((TRDtrack = trkInf->GetTRDtrack()) || trkInf->GetNumberOfClustersRefit()){
 100                         // check if allready found by the tracker
 101                         Bool_t found = kFALSE;
 102                         for(Int_t il = 0; il < naccepted; il++){
 103                                 if(labelsacc[il] == trkInf->GetLabel()) found = kTRUE;
 104                         }
 105                         if(found){
 106                                 mom =  trackRef ? trackRef->P() : TRDtrack->P();
 107                                 contamination->Fill(mom, 1);
 108                                 ndoublecounted++;
 109                                 continue;
 110                         }
 111                         if(trkInf->GetNTrackRefs()){
 112                                 Int_t iref = 0;
 113                                 while(!(trackRef = trkInf->GetTrackRef(iref++)));
 114                         }
 115                         if(!trackRef) printf("Error: Track Reference missing for Track %d\n", TRDtrack->GetLabel());
 116                         mom =  trackRef ? trackRef->P() : TRDtrack->P();
 117 //           Accept track
 118                         if(fDebugLevel > 3)printf("Accept track\n");
 119                         momacc[naccepted] = mom;
 120                         labelsacc[naccepted++] = trkInf->GetLabel();
 121 /*                      printf("Reconstructed: event %3d Tracks: MC[%d] ESD[%d] NRefs[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), trkInf->GetLabel(),  trkInf->GetTrackId(), trkInf->GetNTrackRefs());*/
 122     } else{
 123                         if(fDebugLevel>10) printf("Analysing Track References\n");
 124       // Check if track is findable
 125                         Float_t xmin = 10000.0, xmax = 0.0;
 126                         Float_t ymin = 0.0, ymax = 0.0;
 127                         Float_t zmin = 0.0, zmax = 0.0;
 128                         Float_t lastx = 0.0, x = 0.0;
 129                         Int_t nLayers = 0;
 130 /*                      trackRef = trkInf->GetTrackRef(0);*/
 131 /*                      xmin = trackRef->LocalX(); xmax = trackRef->LocalX();
 132       ymin = trackRef->LocalY(); ymax = trackRef->LocalY();
 133       mom = trackRef->P();*/
 134       Int_t *sector = new Int_t[trkInf->GetNTrackRefs()];
 135       for(Int_t itr = 0; itr < trkInf->GetNTrackRefs(); itr++){
 136         trackRef = trkInf->GetTrackRef(itr);
 137                                 if(fDebugLevel>10) printf("%d. x[%f], y[%f], z[%f]\n", itr, trackRef->LocalX(), trackRef->LocalY(), trackRef->Z());
 138                                 x = trackRef->LocalX();
 139                                 if(x < 250. || x > 370.) continue;      // Be Sure that we are inside TRD
 140         sector[itr] = Int_t(trackRef->Alpha()/kAlpha);
 141         if(x < xmin){
 142           xmin = trackRef->LocalX();
 143           ymin = trackRef->LocalY();
 144           zmin = trackRef->Z();
 145           mom = trackRef->P();
 146         }
 147         if(x > xmax){
 148           xmax = trackRef->LocalX();
 149           ymax = trackRef->LocalY();
 150           zmax = trackRef->Z();
 151         }
 152                                 if(itr > 0){
 153                                         Float_t dist = TMath::Abs(x - lastx);
 154                                         if(fDebugLevel>10) printf("x = %f, lastx = %f, dist = %f\n", x, lastx, dist);
 155                                         if(TMath::Abs(dist - 3.7) < 0.1) nLayers++;     // ref(i+1) has to be larger than ref(i)
 156                                 }
 157                                 lastx = x;
 158       }
 159       // Apply cuts
 160       Bool_t findable = kTRUE;
 161       if(trkInf->GetNTrackRefs() > 2 && xmax > xmin){
 162         if(mom < 0.55) findable = kFALSE;                                                                       // momentum cut at 0.6
 163                                 Double_t yangle = (ymax -ymin)/(xmax - xmin);
 164                                 Double_t zangle = (zmax -zmin)/(xmax - xmin);
 165                                 if(fDebugLevel>10) printf("track: y-Angle = %f, z-Angle = %f\n", yangle, zangle);
 166                                 if(fDebugLevel>10) printf("nLayers = %d\n", nLayers);
 167                                 if(TMath::ATan(TMath::Abs(yangle)) > 45.) findable = kFALSE;
 168                                 if(TMath::ATan(TMath::Abs(zangle)) > 45.) findable = kFALSE;
 169                                 if(nLayers < 4) findable = kFALSE;
 170         if(!trkInf->IsPrimary()) findable = kFALSE;
 171         Bool_t samesec = kTRUE;
 172         for(Int_t iref = 1; iref < trkInf->GetNTrackRefs(); iref++)
 173           if(sector[iref] != sector[0]) samesec = kFALSE;
 174         if(!samesec) findable = kFALSE;         // Discard all tracks which are passing more than one sector
 175         if(fDebugLevel){
 176           Double_t trackAngle = TMath::ATan(yangle);
 177           Bool_t primary = trkInf->IsPrimary();
 178           (*fDebugStream) << "NotFoundTrack"
 179             << "Momentum="      << mom
 180             << "trackAngle="<< trackAngle
 181             << "NLayers="       << nLayers
 182             << "Primary="       << primary
 183             << "\n";
 184         }
 185       }
 186       else
 187         findable = kFALSE;
 188       delete[] sector;
 189       if(findable){
 190                                 momrej[nrejected] = mom;
 191                                 labelsrej[nrejected++] = trkInf->GetLabel();
 192 /*                        printf("Not Reconstructed: event %3d Tracks: MC[%d] ESD[%d] NRefs[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), trkInf->GetLabel(),  trkInf->GetTrackId(), trkInf->GetNTrackRefs());*/
 193       }
 194     }
 195   }
 196         for(Int_t itk = 0; itk < naccepted; itk++){
 197                 if(fDebugLevel > 2)printf("Accepted MC track: %d\n", labelsacc[itk]);
 198                 efficiency->Fill(momacc[itk], 1);
 199                 contamination->Fill(momacc[itk], 0);
 200         }
 201         Int_t nall = naccepted;
 202         for(Int_t imis = 0; imis < nrejected; imis++){
 203                 Bool_t found = kFALSE;
 204                 for(Int_t ifound = 0; ifound < naccepted; ifound++){
 205                         if(labelsacc[ifound] == labelsrej[imis]){
 206                                 found = kTRUE;
 207                                 break;
 208                         }
 209                 }
 210                 if(!found){
 211                         efficiency->Fill(momrej[imis], 0);
 212                         contamination->Fill(momrej[imis], 0);
 213                         if(fDebugLevel > 2)printf("Rejected MC track: %d\n", labelsrej[imis]);
 214                         nall++;
 215                 }
 216         }
 217   //if(fDebugLevel>=1)
 218   printf("%3d Tracks: MC[%3d] TRD[%3d | %5.2f%%] \n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nall, naccepted, 1.E2*Float_t(naccepted)/Float_t(nall));
 219   printf("%3d Tracks: ALL[%3d] DoubleCounted[%3d | %5.2f%%] \n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nall + ndoublecounted, ndoublecounted, 1.E2*Float_t(ndoublecounted)/Float_t(nall + ndoublecounted));
 220
 221   PostData(0, fContainer);
 222 }
 223
 224 //_____________________________________________________________________________
 225 void AliTRDtrackingEfficiencyCombined::Terminate(Option_t *)
 226 {
 227   //
 228   // Termination
 229   //
 230
 231   if(fDebugStream){
 232     delete fDebugStream;
 233     fDebugStream = 0x0;
 234     fDebugLevel = 0;
 235   }
 236
 237   fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
 238   if (!fContainer) {
 239     Printf("ERROR: list not available");
 240     return;
 241   }
 242
 243   /*TProfile *hEff = (TProfile*)fContainer->At(0);
 244   TProfile *hEffCont = (TProfile*)fContainer->At(1);
 245   Printf("Eff[%p] EffCont[%p]\n", (void*)hEff, (void*)hEffCont);
 246
 247
 248   TCanvas *c2 = new TCanvas("c2","",800,400);
 249   c2->Divide(2,1);
 250
 251   c2->cd(1);
 252   hEff->DrawCopy("e1");
 253   c2->cd(2);
 254   hEffCont->DrawCopy("e1");*/
 255 }