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