PWG3/hfe/AliHFEextraCuts.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 // Extra cuts implemented by the ALICE Heavy Flavour Electron Group
  17 // Cuts stored here:
  18 // - ITS pixels
  19 // - TPC cluster ratio
  20 // - TRD tracklets
  21 //
  22 // Authors:
  23 //   Markus Fasel <M.Fasel@gsi.de>
  24 //
  25 #include <TClass.h>
  26 #include <TH1F.h>
  27 #include <TH2F.h>
  28 #include <TList.h>
  29 #include <TString.h>
  30 #include <TMath.h>
  31
  32 #include "AliESDtrack.h"
  33 #include "AliMCParticle.h"
  34
  35 #include "AliHFEextraCuts.h"
  36
  37 ClassImp(AliHFEextraCuts)
  38
  39 //______________________________________________________
  40 AliHFEextraCuts::AliHFEextraCuts(const Char_t *name, const Char_t *title):
  41   AliCFCutBase(name, title),
  42   fCutCorrelation(0),
  43   fRequirements(0),
  44   fClusterRatioTPC(0.),
  45   fMinTrackletsTRD(0),
  46   fPixelITS(0),
  47   fCheck(kTRUE),
  48   fQAlist(0x0),
  49   fDebugLevel(0)
  50 {
  51   //
  52   // Default Constructor
  53   //
  54   memset(fImpactParamCut, 0, sizeof(Float_t) * 4);
  55 }
  56
  57 //______________________________________________________
  58 AliHFEextraCuts::AliHFEextraCuts(const AliHFEextraCuts &c):
  59   AliCFCutBase(c),
  60   fCutCorrelation(c.fCutCorrelation),
  61   fRequirements(c.fRequirements),
  62   fClusterRatioTPC(c.fClusterRatioTPC),
  63   fMinTrackletsTRD(c.fMinTrackletsTRD),
  64   fPixelITS(c.fPixelITS),
  65   fCheck(c.fCheck),
  66   fQAlist(0x0),
  67   fDebugLevel(c.fDebugLevel)
  68 {
  69   //
  70   // Copy constructor
  71   // Performs a deep copy
  72   //
  73   memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
  74   if(IsQAOn()){
  75     fIsQAOn = kTRUE;
  76     fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
  77     fQAlist->SetOwner();
  78   }
  79 }
  80
  81 //______________________________________________________
  82 AliHFEextraCuts &AliHFEextraCuts::operator=(const AliHFEextraCuts &c){
  83   //
  84   // Assignment operator
  85   //
  86   if(this != &c){
  87     AliCFCutBase::operator=(c);
  88     fCutCorrelation = c.fCutCorrelation;
  89     fRequirements = c.fRequirements;
  90     fClusterRatioTPC = c.fClusterRatioTPC;
  91     fMinTrackletsTRD = c.fMinTrackletsTRD;
  92     fPixelITS = c.fPixelITS;
  93     fCheck = c.fCheck;
  94     fDebugLevel = c.fDebugLevel;
  95
  96     memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
  97     if(IsQAOn()){
  98       fIsQAOn = kTRUE;
  99       fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
 100       fQAlist->SetOwner();
 101     }else fQAlist = 0x0;
 102   }
 103   return *this;
 104 }
 105
 106 //______________________________________________________
 107 AliHFEextraCuts::~AliHFEextraCuts(){
 108   //
 109   // Destructor
 110   //
 111   if(fQAlist){
 112     fQAlist->Delete();
 113     delete fQAlist;
 114   }
 115 }
 116
 117 //______________________________________________________
 118 Bool_t AliHFEextraCuts::IsSelected(TObject *o){
 119   //
 120   // Steering function for the track selection
 121   //
 122   if(TString(o->IsA()->GetName()).CompareTo("AliESDtrack") == 0){
 123     return CheckESDCuts(dynamic_cast<AliESDtrack *>(o));
 124   }
 125   return CheckMCCuts(dynamic_cast<AliMCParticle *>(o));
 126 }
 127
 128 //______________________________________________________
 129 Bool_t AliHFEextraCuts::CheckESDCuts(AliESDtrack *track){
 130   //
 131   // Checks cuts on reconstructed tracks
 132   // returns true if track is selected
 133   // QA histograms are filled before track selection and for
 134   // selected tracks after track selection
 135   //
 136   ULong64_t survivedCut = 0;    // Bitmap for cuts which are passed by the track, later to be compared with fRequirements
 137   if(IsQAOn()) FillQAhistosESD(track, kBeforeCuts);
 138   // Apply cuts
 139   Float_t impactR, impactZ, ratioTPC;
 140   track->GetImpactParameters(impactR, impactZ);
 141   // printf("Check TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
 142   ratioTPC = track->GetTPCNclsF() > 0. ? static_cast<Float_t>(track->GetTPCNcls())/static_cast<Float_t>(track->GetTPCNclsF()) : 1.;
 143   UChar_t trdTracklets;
 144   trdTracklets = track->GetTRDntrackletsPID();
 145   UChar_t itsPixel = track->GetITSClusterMap();
 146   Int_t det, status1, status2;
 147   Float_t xloc, zloc;
 148   track->GetITSModuleIndexInfo(0, det, status1, xloc, zloc);
 149   track->GetITSModuleIndexInfo(1, det, status2, xloc, zloc);
 150   if(TESTBIT(fRequirements, kMinImpactParamR)){
 151     // cut on min. Impact Parameter in Radial direction
 152     if(TMath::Abs(impactR) >= fImpactParamCut[0]) SETBIT(survivedCut, kMinImpactParamR);
 153   }
 154   if(TESTBIT(fRequirements, kMinImpactParamZ)){
 155     // cut on min. Impact Parameter in Z direction
 156     if(TMath::Abs(impactZ) >= fImpactParamCut[1]) SETBIT(survivedCut, kMinImpactParamZ);
 157   }
 158   if(TESTBIT(fRequirements, kMaxImpactParamR)){
 159     // cut on max. Impact Parameter in Radial direction
 160     if(TMath::Abs(impactR) <= fImpactParamCut[2]) SETBIT(survivedCut, kMaxImpactParamR);
 161   }
 162   if(TESTBIT(fRequirements, kMaxImpactParamZ)){
 163     // cut on max. Impact Parameter in Z direction
 164     if(TMath::Abs(impactZ) <= fImpactParamCut[3]) SETBIT(survivedCut, kMaxImpactParamZ);
 165   }
 166   if(TESTBIT(fRequirements, kClusterRatioTPC)){
 167     // cut on min ratio of found TPC clusters vs findable TPC clusters
 168     if(ratioTPC >= fClusterRatioTPC) SETBIT(survivedCut, kClusterRatioTPC);
 169   }
 170   if(TESTBIT(fRequirements, kMinTrackletsTRD)){
 171     // cut on minimum number of TRD tracklets
 172     if(fDebugLevel > 0){
 173       printf("Min TRD cut: [%d|%d]\n", fMinTrackletsTRD, trdTracklets);
 174     }
 175     if(trdTracklets >= fMinTrackletsTRD) SETBIT(survivedCut, kMinTrackletsTRD);
 176   }
 177   if(TESTBIT(fRequirements, kPixelITS)){
 178     // cut on ITS pixel layers
 179     if(fDebugLevel > 0){
 180       printf("ITS cluster Map: ");
 181       PrintBitMap(itsPixel);
 182     }
 183     switch(fPixelITS){
 184       case kFirst:
 185         if(!TESTBIT(itsPixel, 0)) {
 186           if(fCheck){
 187             if(!CheckITSstatus(status1)) {
 188               SETBIT(survivedCut, kPixelITS);
 189             }
 190           }
 191         }
 192         else {
 193           SETBIT(survivedCut, kPixelITS);
 194         }
 195                     break;
 196       case kSecond:
 197         if(!TESTBIT(itsPixel, 1)) {
 198           if(fCheck) {
 199             if(!CheckITSstatus(status2)) {
 200               SETBIT(survivedCut, kPixelITS);
 201             }
 202           }
 203         }
 204         else {
 205           SETBIT(survivedCut, kPixelITS);
 206         }
 207                     break;
 208       case kBoth:
 209         if(!(TESTBIT(track->GetITSClusterMap(),0))) {
 210           if(fCheck) {
 211             if(!CheckITSstatus(status1)) {
 212               if(!(TESTBIT(track->GetITSClusterMap(),1))) {
 213                 if(!CheckITSstatus(status2)) {
 214                   SETBIT(survivedCut, kPixelITS);
 215                 }
 216               }
 217               else SETBIT(survivedCut, kPixelITS);
 218             }
 219           }
 220         }
 221         else {
 222
 223           if(!(TESTBIT(track->GetITSClusterMap(),1))) {
 224             if(fCheck) {
 225               if(!CheckITSstatus(status2)) {
 226                 SETBIT(survivedCut, kPixelITS);
 227               }
 228             }
 229           }
 230           else SETBIT(survivedCut, kPixelITS);
 231
 232         }
 233                    break;
 234       case kAny:
 235         if((!TESTBIT(itsPixel, 0)) && (!TESTBIT(itsPixel, 1))){
 236           if(fCheck){
 237             if(!CheckITSstatus(status1) || (!CheckITSstatus(status2))) {
 238               SETBIT(survivedCut, kPixelITS);
 239             }
 240           }
 241         }
 242         else {
 243           SETBIT(survivedCut, kPixelITS);
 244         }
 245                     break;
 246       default: break;
 247     }
 248     if(fDebugLevel > 0)printf("Survived Cut: %s\n", TESTBIT(survivedCut, kPixelITS) ? "YES" : "NO");
 249   }
 250   if(fRequirements == survivedCut){
 251     //
 252     // Track selected
 253     //
 254     if(IsQAOn()) FillQAhistosESD(track, kAfterCuts);
 255     return kTRUE;
 256   }
 257   if(IsQAOn()) FillCutCorrelation(survivedCut);
 258   return kFALSE;
 259 }
 260
 261 //______________________________________________________
 262 Bool_t AliHFEextraCuts::CheckMCCuts(AliMCParticle */*track*/) const {
 263   //
 264   // Checks cuts on Monte Carlo tracks
 265   // returns true if track is selected
 266   // QA histograms are filled before track selection and for
 267   // selected tracks after track selection
 268   //
 269   return kTRUE; // not yet implemented
 270 }
 271
 272 //______________________________________________________
 273 void AliHFEextraCuts::FillQAhistosESD(AliESDtrack *track, UInt_t when){
 274   //
 275   // Fill the QA histograms for ESD tracks
 276   // Function can be called before cuts or after cut application (second argument)
 277   //
 278   TList *container = dynamic_cast<TList *>(fQAlist->At(when));
 279   Float_t impactR, impactZ;
 280   track->GetImpactParameters(impactR, impactZ);
 281   (dynamic_cast<TH1F *>(container->At(0)))->Fill(impactR);
 282   (dynamic_cast<TH1F *>(container->At(1)))->Fill(impactZ);
 283   // printf("TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
 284   (dynamic_cast<TH1F *>(container->At(2)))->Fill(track->GetTPCNclsF() > 0. ? static_cast<Float_t>(track->GetTPCNcls())/static_cast<Float_t>(track->GetTPCNclsF()) : 1.);
 285   (dynamic_cast<TH1F *>(container->At(3)))->Fill(track->GetTRDntrackletsPID());
 286   UChar_t itsPixel = track->GetITSClusterMap();
 287   TH1 *pixelHist = dynamic_cast<TH1F *>(container->At(4));
 288   Int_t firstEntry = pixelHist->GetXaxis()->GetFirst();
 289   if(!((itsPixel & BIT(0)) || (itsPixel & BIT(1))))
 290     pixelHist->Fill(firstEntry + 3);
 291   else{
 292     if(itsPixel & BIT(0)){
 293       pixelHist->Fill(firstEntry);
 294       if(itsPixel & BIT(1)) pixelHist->Fill(firstEntry + 2);
 295       else pixelHist->Fill(firstEntry + 4);
 296     }
 297     if(itsPixel & BIT(1)){
 298       pixelHist->Fill(firstEntry + 1);
 299       if(!(itsPixel & BIT(0))) pixelHist->Fill(firstEntry + 5);
 300     }
 301   }
 302 }
 303
 304 // //______________________________________________________
 305 // void AliHFEextraCuts::FillQAhistosMC(AliMCParticle *track, UInt_t when){
 306 //   //
 307 //   // Fill the QA histograms for MC tracks
 308 //   // Function can be called before cuts or after cut application (second argument)
 309 //   // Not yet implemented
 310 //   //
 311 // }
 312
 313 //______________________________________________________
 314 void AliHFEextraCuts::FillCutCorrelation(ULong64_t survivedCut){
 315   //
 316   // Fill cut correlation histograms for tracks that didn't pass cuts
 317   //
 318   TH2 *correlation = dynamic_cast<TH2F *>(fQAlist->At(2));
 319   for(Int_t icut = 0; icut < kNcuts; icut++){
 320     if(!TESTBIT(fRequirements, icut)) continue;
 321     for(Int_t jcut = icut; jcut < kNcuts; jcut++){
 322       if(!TESTBIT(fRequirements, jcut)) continue;
 323       if(TESTBIT(survivedCut, icut) && TESTBIT(survivedCut, jcut))
 324         correlation->Fill(icut, jcut);
 325     }
 326   }
 327 }
 328
 329 //______________________________________________________
 330 void AliHFEextraCuts::AddQAHistograms(TList *qaList){
 331   //
 332   // Add QA histograms
 333   // For each cut a histogram before and after track cut is created
 334   // Histos before respectively after cut are stored in different lists
 335   // Additionally a histogram with the cut correlation is created and stored
 336   // in the top directory
 337   //
 338   TList *histos[2];
 339   TH1 *histo1D = 0x0;
 340   TH2 *histo2D = 0x0;
 341   histos[0] = new TList();
 342   histos[0]->SetName("BeforeCut");
 343   histos[0]->SetOwner();
 344   histos[1] = new TList();
 345   histos[1]->SetName("AfterCut");
 346   histos[1]->SetOwner();
 347   TString cutstr[2] = {"before", "after"};
 348   for(Int_t icond = 0; icond < 2; icond++){
 349     histos[icond]->AddAt((histo1D = new TH1F(Form("impactParamR%s", cutstr[icond].Data()), "Radial Impact Parameter", 100, 0, 10)), 0);
 350     histo1D->GetXaxis()->SetTitle("Impact Parameter");
 351     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 352     histos[icond]->AddAt((histo1D = new TH1F(Form("impactParamZ%s", cutstr[icond].Data()), "Z Impact Parameter", 200, 0, 20)), 1);
 353     histo1D->GetXaxis()->SetTitle("Impact Parameter");
 354     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 355     histos[icond]->AddAt((histo1D = new TH1F(Form("tpcClr%s", cutstr[icond].Data()), "Cluster Ratio TPC", 10, 0, 1)), 2);
 356     histo1D->GetXaxis()->SetTitle("Cluster Ratio TPC");
 357     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 358     histos[icond]->AddAt((histo1D = new TH1F(Form("trdTracklets%s", cutstr[icond].Data()), "Number of TRD tracklets", 7, 0, 7)), 3);
 359     histo1D->GetXaxis()->SetTitle("Number of TRD Tracklets");
 360     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 361     histos[icond]->AddAt((histo1D = new TH1F(Form("itsPixel%s", cutstr[icond].Data()), "ITS Pixel Hits", 6, 0, 6)), 4);
 362     histo1D->GetXaxis()->SetTitle("ITS Pixel");
 363     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 364     Int_t first = histo1D->GetXaxis()->GetFirst();
 365     TString binNames[6] = { "First", "Second", "Both", "None", "Exclusive First", "Exclusive Second"};
 366     for(Int_t ilabel = 0; ilabel < 6; ilabel++)
 367       histo1D->GetXaxis()->SetBinLabel(first + ilabel, binNames[ilabel].Data());
 368   }
 369   fQAlist = new TList();
 370   fQAlist->SetOwner();
 371   fQAlist->SetName("HFelectronExtraCuts");
 372   fQAlist->AddAt(histos[0], 0);
 373   fQAlist->AddAt(histos[1], 1);
 374   // Add cut correlation
 375   fQAlist->AddAt((histo2D = new TH2F("cutcorrellation", "Cut Correlation", kNcuts, 0, kNcuts - 1, kNcuts, 0, kNcuts -1)), 2);
 376   TString labels[kNcuts] = {"MinImpactParamR", "MaxImpactParamR", "MinImpactParamZ", "MaxImpactParamZ", "ClusterRatioTPC", "MinTrackletsTRD", "ITSpixel"};
 377   Int_t firstx = histo2D->GetXaxis()->GetFirst(), firsty = histo2D->GetYaxis()->GetFirst();
 378   for(Int_t icut = 0; icut < kNcuts; icut++){
 379     histo2D->GetXaxis()->SetBinLabel(firstx + icut, labels[icut].Data());
 380     histo2D->GetYaxis()->SetBinLabel(firsty + icut, labels[icut].Data());
 381   }
 382   qaList->AddLast(fQAlist);
 383 }
 384
 385 //______________________________________________________
 386 void AliHFEextraCuts::PrintBitMap(Int_t bitmap){
 387   for(Int_t ibit = 32; ibit--; )
 388     printf("%d", bitmap & BIT(ibit) ? 1 : 0);
 389   printf("\n");
 390 }
 391
 392 //______________________________________________________
 393 Bool_t AliHFEextraCuts::CheckITSstatus(Int_t itsStatus) const {
 394   //
 395   // Check whether ITS area is dead
 396   //
 397   Bool_t status;
 398   switch(itsStatus){
 399     case 2: status = kFALSE; break;
 400     case 3: status = kFALSE; break;
 401     case 7: status = kFALSE; break;
 402     default: status = kTRUE;
 403   }
 404   return status;
 405 }