PWG0/esdTrackCuts/ESDtrackQualityCuts.cxx

   1 #include "ESDtrackQualityCuts.h"
   2
   3 #include <Riostream.h>
   4
   5 //____________________________________________________________________
   6 ClassImp(ESDtrackQualityCuts);
   7
   8 //____________________________________________________________________
   9 ESDtrackQualityCuts::ESDtrackQualityCuts() {
  10
  11   //##############################################################################
  12   // setting default cuts
  13   SetMinNClustersTPC();
  14   SetMinNClustersITS();
  15   SetMaxChi2PerClusterTPC();
  16   SetMaxChi2PerClusterITS();
  17   SetMaxCovDiagonalElements();
  18   SetRequireTPCRefit();
  19   SetRequireITSRefit();
  20   SetAcceptKingDaughters();
  21   SetMinNsigmaToVertex();
  22   SetRequireSigmaToVertex();
  23
  24   SetHistogramsOn();
  25
  26   // set the cut names
  27   fCutNames[0]  = "require TPC refit";
  28   fCutNames[1]  = "require ITS refit";
  29   fCutNames[2]  = "n clusters TPC";
  30   fCutNames[3]  = "n clusters ITS";
  31   fCutNames[4]  = "#Chi^{2}/clusters TPC";
  32   fCutNames[5]  = "#Chi^{2}/clusters ITS";
  33   fCutNames[6]  = "cov 11";
  34   fCutNames[7]  = "cov 22";
  35   fCutNames[8]  = "cov 33";
  36   fCutNames[9]  = "cov 44";
  37   fCutNames[10] = "cov 55";
  38   fCutNames[11] = "trk-to-vtx";
  39   fCutNames[12] = "trk-to-vtx failed";
  40   fCutNames[13] = "kink daughters";
  41 }
  42
  43 //____________________________________________________________________
  44 Bool_t
  45 ESDtrackQualityCuts::AcceptTrack(AliESDtrack* esdTrack, AliESDVertex* esdVtx, Double_t field) {
  46
  47   // relate the track to the new vertex
  48   esdTrack->RelateToVertex(esdVtx, field, 999999);
  49
  50   return AcceptTrack(esdTrack);
  51 }
  52
  53 //____________________________________________________________________
  54 Bool_t
  55 ESDtrackQualityCuts::AcceptTrack(AliESDtrack* esdTrack, Double_t* vtx, Double_t* vtx_res, Double_t field) {
  56
  57   AliESDVertex* esdVtx = new AliESDVertex(vtx, vtx_res,"new vertex");
  58   esdTrack->RelateToVertex(esdVtx, field, 1e99);
  59   return AcceptTrack(esdTrack);
  60 }
  61
  62 //____________________________________________________________________
  63 Bool_t
  64 ESDtrackQualityCuts::AcceptTrack(AliESDtrack* esdTrack) {
  65
  66   UInt_t status = esdTrack->GetStatus();
  67
  68   // getting variables from ESD
  69   Int_t nClustersITS = esdTrack->GetITSclusters(fIdxInt);
  70   Int_t nClustersTPC = esdTrack->GetTPCclusters(fIdxInt);
  71
  72   Float_t chi2PerClusterITS = -1;
  73   Float_t chi2PerClusterTPC = -1;
  74   if (nClustersITS!=0)
  75     chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
  76   if (nClustersTPC!=0)
  77     chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
  78
  79   Double_t extCov[15];
  80   esdTrack->GetExternalCovariance(extCov);
  81
  82   Float_t b[2];
  83   Float_t bRes[2];
  84   Float_t bCov[3];
  85   esdTrack->GetImpactParameters(b,bCov);
  86   if (bCov[0]<=0 || bCov[2]<=0) {
  87     AliDebug(1, "Estimated b resolution zero!");
  88     bCov[0]=0; bCov[1]=0;
  89   }
  90   bRes[0] = TMath::Sqrt(bCov[0]);
  91   bRes[1] = TMath::Sqrt(bCov[2]);
  92
  93   //
  94   // FIX !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  95   //
  96   // this is not correct - it will not give n sigma!!!
  97   //
  98   //
  99   Float_t nSigmaToVertex = -1;
 100   if (bRes[0]!=0 && bRes[1]!=0)
 101     nSigmaToVertex = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 102
 103
 104   //########################################################################
 105   // cut the track?
 106
 107   Bool_t cuts[fNCuts];
 108   for (Int_t i=0; i<fNCuts; i++) cuts[i]=kFALSE;
 109
 110   if (fCut_RequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
 111     cuts[0]=kTRUE;
 112   if (fCut_RequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
 113     cuts[1]=kTRUE;
 114   if (nClustersTPC<fCut_MinNClusterTPC)
 115     cuts[2]=kTRUE;
 116   if (nClustersITS<fCut_MinNClusterITS)
 117     cuts[3]=kTRUE;
 118   if (chi2PerClusterTPC>fCut_MaxChi2PerClusterTPC)
 119     cuts[4]=kTRUE;
 120   if (chi2PerClusterITS>fCut_MaxChi2PerClusterITS)
 121     cuts[5]=kTRUE;
 122   if (extCov[0]  > fCut_MaxC11)
 123     cuts[6]=kTRUE;
 124   if (extCov[2]  > fCut_MaxC22)
 125     cuts[7]=kTRUE;
 126   if (extCov[5]  > fCut_MaxC33)
 127     cuts[8]=kTRUE;
 128   if (extCov[9]  > fCut_MaxC44)
 129     cuts[9]=kTRUE;
 130   if (extCov[14]  > fCut_MaxC55)
 131     cuts[10]=kTRUE;
 132   if (nSigmaToVertex > fCut_NsigmaToVertex)
 133     cuts[11] = kTRUE;
 134   // if n sigma could not be calculated
 135   if (nSigmaToVertex==-1 && fCut_SigmaToVertexRequired)
 136     cuts[12]=kTRUE;
 137   if (!fCut_AcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
 138     cuts[13]=kTRUE;
 139
 140   Bool_t cut=kFALSE;
 141   for (Int_t i=0; i<fNCuts; i++)
 142     if (cuts[i]) cut = kTRUE;
 143
 144   //########################################################################
 145   // filling histograms
 146   if (fHistogramsOn) {
 147     hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin("n tracks")));
 148
 149     if (cut)
 150       hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin("n cut tracks")));
 151
 152     for (Int_t i=0; i<fNCuts; i++) {
 153       if (cuts[i])
 154         hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin(fCutNames[i])));
 155
 156       for (Int_t j=i; j<fNCuts; j++) {
 157         if (cuts[i] && cuts[j]) {
 158           Float_t x = hCutCorrelation->GetXaxis()->GetBinCenter(hCutCorrelation->GetXaxis()->FindBin(fCutNames[i]));
 159           Float_t y = hCutCorrelation->GetYaxis()->GetBinCenter(hCutCorrelation->GetYaxis()->FindBin(fCutNames[j]));
 160           hCutCorrelation->Fill(x,y);
 161         }
 162       }
 163     }
 164
 165
 166     hNClustersITS[0]->Fill(nClustersITS);
 167     hNClustersTPC[0]->Fill(nClustersTPC);
 168     hChi2PerClusterITS[0]->Fill(chi2PerClusterITS);
 169     hChi2PerClusterTPC[0]->Fill(chi2PerClusterTPC);
 170
 171     hC11[0]->Fill(extCov[0]);
 172     hC22[0]->Fill(extCov[2]);
 173     hC33[0]->Fill(extCov[5]);
 174     hC44[0]->Fill(extCov[9]);
 175     hC55[0]->Fill(extCov[14]);
 176
 177     hDZ[0]->Fill(b[1]);
 178     hDXY[0]->Fill(b[0]);
 179     hDXYvsDZ[0]->Fill(b[1],b[0]);
 180
 181     if (bRes[0]!=0 && bRes[1]!=0) {
 182       hDZNormalized[0]->Fill(b[1]/bRes[1]);
 183       hDXYNormalized[0]->Fill(b[0]/bRes[0]);
 184       hDXYvsDZNormalized[0]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
 185     }
 186   }
 187
 188   //########################################################################
 189   // cut the track!
 190   if (cut) return kFALSE;
 191
 192   //########################################################################
 193   // filling histograms after cut
 194   if (fHistogramsOn) {
 195     hNClustersITS[1]->Fill(nClustersITS);
 196     hNClustersTPC[1]->Fill(nClustersTPC);
 197     hChi2PerClusterITS[1]->Fill(chi2PerClusterITS);
 198     hChi2PerClusterTPC[1]->Fill(chi2PerClusterTPC);
 199
 200     hC11[1]->Fill(extCov[0]);
 201     hC22[1]->Fill(extCov[2]);
 202     hC33[1]->Fill(extCov[5]);
 203     hC44[1]->Fill(extCov[9]);
 204     hC55[1]->Fill(extCov[14]);
 205
 206     hDZ[1]->Fill(b[1]);
 207     hDXY[1]->Fill(b[0]);
 208     hDXYvsDZ[1]->Fill(b[1],b[0]);
 209
 210     hDZNormalized[1]->Fill(b[1]/bRes[1]);
 211     hDXYNormalized[1]->Fill(b[0]/bRes[0]);
 212     hDXYvsDZNormalized[1]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
 213   }
 214
 215   return kTRUE;
 216 }
 217
 218 //____________________________________________________________________
 219 void
 220 ESDtrackQualityCuts::DefineHistograms(Int_t color) {
 221
 222   fHistogramsOn=kTRUE;
 223
 224   //###################################################################################
 225   // defining histograms
 226
 227   hCutStatistics = new TH1F("cut_statistics","cut statistics",fNCuts+4,-0.5,fNCuts+3.5);
 228
 229   hCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
 230   hCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
 231
 232   hCutCorrelation = new TH2F("cut_correlation","cut correlation",fNCuts,-0.5,fNCuts-0.5,fNCuts,-0.5,fNCuts-0.5);;
 233
 234   for (Int_t i=0; i<fNCuts; i++) {
 235     hCutStatistics->GetXaxis()->SetBinLabel(i+4,fCutNames[i]);
 236     hCutCorrelation->GetXaxis()->SetBinLabel(i+1,fCutNames[i]);
 237     hCutCorrelation->GetYaxis()->SetBinLabel(i+1,fCutNames[i]);
 238   }
 239
 240   hCutStatistics  ->SetLineColor(color);
 241   hCutCorrelation ->SetLineColor(color);
 242   hCutStatistics  ->SetLineWidth(2);
 243   hCutCorrelation ->SetLineWidth(2);
 244
 245
 246   hNClustersITS        = new TH1F*[2];
 247   hNClustersTPC        = new TH1F*[2];
 248   hChi2PerClusterITS   = new TH1F*[2];
 249   hChi2PerClusterTPC   = new TH1F*[2];
 250
 251   hC11                 = new TH1F*[2];
 252   hC22                 = new TH1F*[2];
 253   hC33                 = new TH1F*[2];
 254   hC44                 = new TH1F*[2];
 255   hC55                 = new TH1F*[2];
 256
 257   hDXY                 = new TH1F*[2];
 258   hDZ                  = new TH1F*[2];
 259   hDXYvsDZ             = new TH2F*[2];
 260
 261   hDXYNormalized       = new TH1F*[2];
 262   hDZNormalized        = new TH1F*[2];
 263   hDXYvsDZNormalized   = new TH2F*[2];
 264
 265
 266
 267   Char_t str[256];
 268   for (Int_t i=0; i<2; i++) {
 269     if (i==0) sprintf(str,"");
 270     else sprintf(str,"_cut");
 271
 272     hNClustersITS[i]        = new TH1F(Form("nClustersITS%s",str),"",8,-0.5,7.5);
 273     hNClustersTPC[i]        = new TH1F(Form("nClustersTPC%s",str),"",165,-0.5,164.5);
 274     hChi2PerClusterITS[i]   = new TH1F(Form("chi2PerClusterITS%s",str),"",500,0,10);
 275     hChi2PerClusterTPC[i]   = new TH1F(Form("chi2PerClusterTPC%s",str),"",500,0,10);
 276
 277     hC11[i]                 = new TH1F(Form("covMatrixDiagonal11%s",str),"",1000,0,5);
 278     hC22[i]                 = new TH1F(Form("covMatrixDiagonal22%s",str),"",1000,0,5);
 279     hC33[i]                 = new TH1F(Form("covMatrixDiagonal33%s",str),"",1000,0,0.5);
 280     hC44[i]                 = new TH1F(Form("covMatrixDiagonal44%s",str),"",1000,0,5);
 281     hC55[i]                 = new TH1F(Form("covMatrixDiagonal55%s",str),"",1000,0,5);
 282
 283     hDXY[i]                 = new TH1F(Form("dXY%s",str),"",500,-10,10);
 284     hDZ[i]                  = new TH1F(Form("dZ%s",str),"",500,-10,10);
 285     hDXYvsDZ[i]             = new TH2F(Form("dXYvsDZ%s",str),"",200,-10,10,200,-10,10);
 286
 287     hDXYNormalized[i]       = new TH1F(Form("dXYNormalized%s",str),"",500,-10,10);
 288     hDZNormalized[i]        = new TH1F(Form("dZNormalized%s",str),"",500,-10,10);
 289     hDXYvsDZNormalized[i]   = new TH2F(Form("dXYvsDZNormalized%s",str),"",200,-10,10,200,-10,10);
 290
 291
 292     hNClustersITS[i]        ->SetXTitle("n ITS clusters");
 293     hNClustersTPC[i]        ->SetXTitle("n TPC clusters");
 294     hChi2PerClusterITS[i]   ->SetXTitle("#Chi^{2} per ITS cluster");
 295     hChi2PerClusterTPC[i]   ->SetXTitle("#Chi^{2} per TPC cluster");
 296
 297     hC11[i]                 ->SetXTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]");
 298     hC22[i]                 ->SetXTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]");
 299     hC33[i]                 ->SetXTitle("cov 33 : #sigma_{sin(#phi)}^{2}");
 300     hC44[i]                 ->SetXTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}");
 301     hC55[i]                 ->SetXTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]");
 302
 303     hDXY[i]                 ->SetXTitle("transverse impact parameter");
 304     hDZ[i]                  ->SetXTitle("longitudinal impact parameter");
 305     hDXYvsDZ[i]             ->SetXTitle("longitudinal impact parameter");
 306     hDXYvsDZ[i]             ->SetYTitle("transverse impact parameter");
 307
 308     hDXYNormalized[i]       ->SetXTitle("normalized trans impact par");
 309     hDZNormalized[i]        ->SetXTitle("normalized long impact par");
 310     hDXYvsDZNormalized[i]   ->SetXTitle("normalized long impact par");
 311     hDXYvsDZNormalized[i]   ->SetYTitle("normalized trans impact par");
 312
 313     hNClustersITS[i]        ->SetLineColor(color);   hNClustersITS[i]        ->SetLineWidth(2);
 314     hNClustersTPC[i]        ->SetLineColor(color);   hNClustersTPC[i]        ->SetLineWidth(2);
 315     hChi2PerClusterITS[i]   ->SetLineColor(color);   hChi2PerClusterITS[i]   ->SetLineWidth(2);
 316     hChi2PerClusterTPC[i]   ->SetLineColor(color);   hChi2PerClusterTPC[i]   ->SetLineWidth(2);
 317
 318     hC11[i]                 ->SetLineColor(color);   hC11[i]                 ->SetLineWidth(2);
 319     hC22[i]                 ->SetLineColor(color);   hC22[i]                 ->SetLineWidth(2);
 320     hC33[i]                 ->SetLineColor(color);   hC33[i]                 ->SetLineWidth(2);
 321     hC44[i]                 ->SetLineColor(color);   hC44[i]                 ->SetLineWidth(2);
 322     hC55[i]                 ->SetLineColor(color);   hC55[i]                 ->SetLineWidth(2);
 323
 324     hDXY[i]                 ->SetLineColor(color);   hDXY[i]                 ->SetLineWidth(2);
 325     hDZ[i]                  ->SetLineColor(color);   hDZ[i]                  ->SetLineWidth(2);
 326
 327     hDXYNormalized[i]       ->SetLineColor(color);   hDXYNormalized[i]       ->SetLineWidth(2);
 328     hDZNormalized[i]        ->SetLineColor(color);   hDZNormalized[i]        ->SetLineWidth(2);
 329
 330   }
 331 }
 332
 333 //____________________________________________________________________
 334 void
 335 ESDtrackQualityCuts::SaveHistograms(Char_t* dir) {
 336
 337   if (!fHistogramsOn) {
 338     AliDebug(0, "Histograms not on - cannot save histograms!!!");
 339     return;
 340   }
 341
 342   gDirectory->mkdir(dir);
 343   gDirectory->cd(dir);
 344
 345   gDirectory->mkdir("before_cuts");
 346   gDirectory->mkdir("after_cuts");
 347
 348   hCutStatistics->Write();
 349   hCutCorrelation->Write();
 350
 351   for (Int_t i=0; i<2; i++) {
 352     if (i==0)
 353       gDirectory->cd("before_cuts");
 354     else
 355       gDirectory->cd("after_cuts");
 356
 357     hNClustersITS[i]        ->Write();
 358     hNClustersTPC[i]        ->Write();
 359     hChi2PerClusterITS[i]   ->Write();
 360     hChi2PerClusterTPC[i]   ->Write();
 361
 362     hC11[i]                 ->Write();
 363     hC22[i]                 ->Write();
 364     hC33[i]                 ->Write();
 365     hC44[i]                 ->Write();
 366     hC55[i]                 ->Write();
 367
 368     hDXY[i]                 ->Write();
 369     hDZ[i]                  ->Write();
 370     hDXYvsDZ[i]             ->Write();
 371
 372     hDXYNormalized[i]       ->Write();
 373     hDZNormalized[i]        ->Write();
 374     hDXYvsDZNormalized[i]   ->Write();
 375
 376     gDirectory->cd("../");
 377   }
 378
 379   gDirectory->cd("../");
 380 }
 381
 382
 383