PWG3/hfe/AliHFEtrackFilter.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 // Track filter class
  17 // Apply cut steps to all tracks in one event and returns a list of
  18 // filtered tracks
  19 //
  20 // Author:
  21 //   Markus Fasel <M.Fasel@gsi.de>
  22 //
  23 #include <TMath.h>
  24 #include <TString.h>
  25
  26 #include "AliAnalysisCuts.h"
  27 #include "AliCFContainer.h"
  28 #include "AliCFParticleGenCuts.h"
  29 #include "AliCFTrackKineCuts.h"
  30 #include "AliCFTrackIsPrimaryCuts.h"
  31 #include "AliCFTrackQualityCuts.h"
  32 #include "AliESDEvent.h"
  33 #include "AliESDtrack.h"
  34 #include "AliMCEvent.h"
  35 #include "AliMCParticle.h"
  36
  37 #include "AliHFEcontainer.h"
  38 #include "AliHFEcutStep.h"
  39 #include "AliHFEextraCuts.h"
  40 #include "AliHFEtrackFilter.h"
  41 #include "AliHFEtools.h"
  42
  43 ClassImp(AliHFEtrackFilter)
  44
  45 //__________________________________________________________________
  46 AliHFEtrackFilter::AliHFEtrackFilter(const Char_t *name) :
  47   TNamed(name, ""),
  48   fFilteredTracks(NULL),
  49   fCutSteps(NULL),
  50   fEfficiencyContainers(NULL),
  51   fMC(NULL),
  52   fMCsignal(NULL),
  53   fPtBins(0),
  54   fEtaBins(0),
  55   fPhiBins(0),
  56   fPtBinning(NULL),
  57   fEtaBinning(NULL),
  58   fPhiBinning(NULL)
  59 {
  60   //
  61   // Constructor
  62   //
  63   fCutSteps = new TObjArray;
  64   fFilteredTracks = new TObjArray;
  65   fEfficiencyContainers = new TObjArray(4);
  66 }
  67
  68 //__________________________________________________________________
  69 AliHFEtrackFilter::AliHFEtrackFilter(const AliHFEtrackFilter &ref):
  70   TNamed(ref),
  71   fFilteredTracks(NULL),
  72   fCutSteps(NULL),
  73   fEfficiencyContainers(NULL),
  74   fMC(NULL),
  75   fMCsignal(NULL),
  76   fPtBins(0),
  77   fEtaBins(0),
  78   fPhiBins(0),
  79   fPtBinning(NULL),
  80   fEtaBinning(NULL),
  81   fPhiBinning(NULL)
  82 {
  83   //
  84   // Copy constructor
  85   //
  86   ref.Copy(*this);
  87 }
  88
  89 //__________________________________________________________________
  90 AliHFEtrackFilter &AliHFEtrackFilter::operator=(const AliHFEtrackFilter &o){
  91   //
  92   // Assignment operator
  93   //
  94   if(this != &o){
  95     o.Copy(*this);
  96   }
  97   return *this;
  98 }
  99
 100 //__________________________________________________________________
 101 AliHFEtrackFilter::~AliHFEtrackFilter(){
 102   //
 103   // Destructor
 104   //
 105
 106   // remove correction framework content
 107   if(fPtBinning) delete fPtBinning;
 108   if(fEtaBinning) delete fEtaBinning;
 109   if(fPhiBinning) delete fPhiBinning;
 110   if(TestBit(kOwnCFContainers))
 111     fEfficiencyContainers->Delete();
 112   delete fEfficiencyContainers;
 113   // remove default content
 114   delete fFilteredTracks;
 115   delete fCutSteps;
 116 }
 117
 118 //__________________________________________________________________
 119 void AliHFEtrackFilter::Copy(TObject &o) const{
 120   //
 121   // Copy content into the object o
 122   //
 123   TNamed::Copy(o);
 124   AliHFEtrackFilter &target = dynamic_cast<AliHFEtrackFilter &>(o);
 125
 126   // Make copy
 127   target.fFilteredTracks = dynamic_cast<TObjArray *>(fFilteredTracks->Clone());
 128   target.fCutSteps = dynamic_cast<TObjArray *>(fCutSteps->Clone());
 129   target.fEfficiencyContainers = dynamic_cast<TObjArray *>(fEfficiencyContainers->Clone());
 130   target.fMC = fMC;
 131   target.fMCsignal = dynamic_cast<AliHFEcutStep *>(fMCsignal->Clone());
 132
 133   target.fPtBins = fPtBins;
 134   target.fEtaBins = fEtaBins;
 135   target.fPhiBins = fPhiBins;
 136   target.fPtBinning = new Double_t[fPtBins];
 137   memcpy(target.fPtBinning, fPtBinning, sizeof(Double_t) * fPtBins);
 138   target.fEtaBinning = new Double_t[fEtaBins];
 139   memcpy(target.fEtaBinning, fEtaBinning, sizeof(Double_t) * fEtaBins);
 140   target.fPhiBinning = new Double_t[fPhiBins];
 141   memcpy(target.fPhiBinning, fPhiBinning, sizeof(Double_t) * fPhiBins);
 142 }
 143
 144 //__________________________________________________________________
 145 void AliHFEtrackFilter::InitCF(){
 146   //
 147   // Initialize correction framework container
 148   // No HFE container given
 149   // Only for testing purpose
 150   //
 151    const Char_t *cnames[4] = {"_container", "_container_signal", "_container_MC", "_container_signalMC"};
 152   const Char_t *ctitlesAppends[4] = {"(All Rec)", "(All Signals)", "(All Rec (MC))", "(All Signals(MC))"};
 153
 154   // Create the binning if not done from outside
 155   if(!fPtBins){
 156     fPtBins = 40;
 157     fPtBinning = AliHFEtools::MakeLogarithmicBinning(fPtBins, 0.1, 10);
 158   }
 159   if(!fEtaBins){
 160     fEtaBins = 8;
 161     fEtaBinning = AliHFEtools::MakeLinearBinning(fEtaBins, -0.9, 0.9);
 162   }
 163   if(!fPhiBins){
 164     fPhiBins = 18;
 165     fPhiBinning = AliHFEtools::MakeLinearBinning(fPhiBins, 0, 2*TMath::Pi());
 166   }
 167   Double_t chargeBins[3] = {-1.1, 0., 1.1};
 168   Int_t nStep = fCutSteps->GetEntriesFast()+1;
 169   Int_t nBins[4] = {fPtBins, fEtaBins, fPhiBins + 2};
 170   AliCFContainer *ctmp = NULL;
 171   for(Int_t icont = 0; icont < 4; icont++){
 172     TString containername =  GetName() + TString(cnames[icont]);
 173     TString containertitle = TString("Container for filter ") + GetName() + TString(ctitlesAppends[icont]);
 174     ctmp = new AliCFContainer(containername.Data(), containertitle.Data(), nStep, 4, nBins);
 175     SetBit(kOwnCFContainers, kTRUE);
 176
 177     // Set the binning
 178     ctmp->SetBinLimits(0, fPtBinning);
 179     ctmp->SetBinLimits(1, fEtaBinning);
 180     ctmp->SetBinLimits(2, fPhiBinning);
 181     ctmp->SetBinLimits(3, chargeBins);
 182
 183     // Label variable names
 184     ctmp->SetVarTitle(0, "pt");
 185     ctmp->SetVarTitle(1, "eta");
 186     ctmp->SetVarTitle(2, "phi");
 187     ctmp->SetVarTitle(3, "charge");
 188
 189     // Label step name
 190     ctmp->SetStepTitle(0, "No Cuts");
 191     AliHFEcutStep *cutStep = NULL;
 192     for(Int_t istep = 0; istep < fCutSteps->GetEntriesFast(); istep++){
 193       cutStep = dynamic_cast<AliHFEcutStep *>(fCutSteps->UncheckedAt(istep));
 194       if(cutStep) ctmp->SetStepTitle(istep + 1, cutStep->GetName());
 195     }
 196     fEfficiencyContainers->AddAt(ctmp, 0);
 197   }
 198   OwnContainers();
 199 }
 200
 201 //__________________________________________________________________
 202 void AliHFEtrackFilter::InitCF(AliHFEcontainer *cont){
 203   //
 204   // Initialize Correction Framework container
 205   // Use the given HFE container
 206   // HFE standard way to use the correction Framework
 207   //
 208
 209   const Char_t *cnames[4] = {"_container", "_container_signal", "_container_MC", "_container_signalMC"};
 210   const Char_t *ctitlesAppends[4] = {"(All Rec)", "(All Signals)", "(All Rec (MC))", "(All Signals(MC))"};
 211   Int_t nStep = fCutSteps->GetEntriesFast()+1;
 212   AliCFContainer *ctmp = NULL;
 213   for(Int_t icont = 0; icont < 4; icont++){
 214     TString contname = GetName() + TString(cnames[icont]);
 215     TString contTitle = TString("Container for filter ") + GetName() + TString(ctitlesAppends[icont]);
 216     //printf("Adding container %s: %s\n", contname.Data(), contTitle.Data());
 217     cont->CreateContainer(contname.Data(), contTitle.Data(), nStep);
 218     fEfficiencyContainers->AddAt((ctmp = cont->GetCFContainer(contname.Data())), icont);
 219
 220     // Label step name
 221     ctmp->SetStepTitle(0, "No Cuts");
 222     AliHFEcutStep *cutStep = NULL;
 223     for(Int_t istep = 0; istep < fCutSteps->GetEntriesFast(); istep++){
 224       cutStep = dynamic_cast<AliHFEcutStep *>(fCutSteps->UncheckedAt(istep));
 225       if(cutStep) ctmp->SetStepTitle(istep + 1, cutStep->GetName());
 226     }
 227   }
 228   ReleaseContainers();
 229 }
 230
 231 //__________________________________________________________________
 232 void AliHFEtrackFilter::FilterTracks(AliESDEvent * const event){
 233   //
 234   // Perform track filtering
 235   // Check each cut step one by one and select tracks which pass
 236   // all cuts. If the correction framework is initialized, a correction#
 237   // framework container will be filled for each cut step.
 238   //
 239   AliHFEcutStep *cutStep = NULL;
 240   AliESDtrack *track = NULL;
 241   AliMCParticle *mctrack = NULL;
 242   Double_t cont[4] = {0., 0., 0., 0.}, contMC[4] = {0., 0., 0., 0.};
 243   Bool_t goodTrack = kTRUE, signal = kFALSE;
 244   Int_t nStep = fCutSteps->GetEntriesFast();
 245   AliCFContainer *call = dynamic_cast<AliCFContainer *>(fEfficiencyContainers->At(0));
 246   AliCFContainer *cSignal = dynamic_cast<AliCFContainer *>(fEfficiencyContainers->At(1));
 247   AliCFContainer *callMC = dynamic_cast<AliCFContainer *>(fEfficiencyContainers->At(2));
 248   AliCFContainer *cSignalMC = dynamic_cast<AliCFContainer *>(fEfficiencyContainers->At(3));
 249   for(Int_t itrack = 0; itrack < event->GetNumberOfTracks(); itrack++){
 250     signal = kFALSE;
 251     track = event->GetTrack(itrack);
 252     // check Monte-Carlo Information if available
 253     if(fMC){
 254       mctrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(TMath::Abs(track->GetLabel())));
 255       if(mctrack){
 256         //AliMCParticle *mother = dynamic_cast<AliMCParticle *>(fMC->GetTrack(mctrack->Particle()->GetFirstMother()));
 257         //AliInfo(Form("Label %d, Mother %d", track->GetLabel(), mother->Particle()->GetPdgCode()));
 258         contMC[0] = mctrack->Pt();
 259         contMC[1] = mctrack->Eta();
 260         contMC[2] = mctrack->Phi();
 261         contMC[3] = mctrack->Charge();
 262         if(fMCsignal->IsSelected(mctrack)) signal = kTRUE;
 263         //if(TMath::Abs(mother->Particle()->GetPdgCode()) != 443) signal = kFALSE;
 264         //AliInfo(Form("Signal? %s", signal ? "Yes": "No"));
 265       }
 266     }
 267
 268     // Fill Array without cut
 269     cont[0] = track->Pt();
 270     cont[1] = track->Eta();
 271     cont[2] = track->Phi();
 272     cont[3] = track->Charge();
 273     if(call) call->Fill(cont, 0);
 274     if(callMC) callMC->Fill(contMC, 0);
 275     if(signal){
 276       if(cSignal) cSignal->Fill(cont, 0);
 277       if(cSignalMC) cSignalMC->Fill(contMC, 0);
 278     }
 279     // cut the track
 280     goodTrack = kTRUE;
 281     for(Int_t icut = 0; icut < nStep; icut++){
 282       cutStep = dynamic_cast<AliHFEcutStep *>(fCutSteps->UncheckedAt(icut));
 283       if(!cutStep->IsSelected(track)){
 284         // track cut away
 285         goodTrack = kFALSE;
 286         break;
 287       }
 288
 289       // Track survived cut step: Fill container
 290       if(call) call->Fill(cont, icut + 1);
 291       if(callMC) callMC->Fill(contMC, icut + 1);
 292       if(signal){
 293         if(cSignal) cSignal->Fill(cont, icut + 1);
 294         if(cSignalMC) cSignalMC->Fill(contMC, icut + 1);
 295       }
 296     }
 297
 298     // Append track to the list of filtered tracks
 299     if(goodTrack) fFilteredTracks->Add(track);
 300   }
 301 }
 302
 303 //__________________________________________________________________
 304 void AliHFEtrackFilter::GenerateCutSteps(){
 305   //
 306   // Make the default cuts
 307   //
 308   MakeCutStepRecKineITSTPC();
 309   MakeCutStepPrimary();
 310   MakeCutStepHFEITS();
 311   MakeCutStepHFETRD();
 312
 313   MakeMCSignalCuts();
 314 }
 315
 316 //__________________________________________________________________
 317 void AliHFEtrackFilter::AddCutStep(AliHFEcutStep *step){
 318   //
 319   // Add new cut step to the filter
 320   //
 321   if(!fCutSteps->FindObject(step->GetName())) fCutSteps->Add(step);
 322 }
 323
 324 //__________________________________________________________________
 325 AliHFEcutStep *AliHFEtrackFilter::GetCutStep(Int_t istep){
 326   //
 327   // Getter for single cut step
 328   //
 329   return dynamic_cast<AliHFEcutStep *>(fCutSteps->At(istep));
 330 }
 331
 332 //__________________________________________________________________
 333 AliHFEcutStep *AliHFEtrackFilter::GetCutStep(const Char_t *name){
 334   //
 335   // Getter for single cut step (by step name)
 336   //
 337   return dynamic_cast<AliHFEcutStep *>(fCutSteps->FindObject(name));
 338 }
 339
 340 //__________________________________________________________________
 341 void AliHFEtrackFilter::Flush(){
 342   //
 343   // Empty track container
 344   //
 345   fFilteredTracks->Clear();
 346 }
 347
 348 //__________________________________________________________________
 349 void AliHFEtrackFilter::SetPtBins(Int_t nBins, Double_t *binning){
 350   //
 351   // User defined pt binning
 352   //
 353   fPtBins = nBins;
 354   fPtBinning = new Double_t[fPtBins + 1];
 355   memcpy(fPtBinning, binning, sizeof(Double_t) * nBins);
 356 }
 357
 358 //__________________________________________________________________
 359 void AliHFEtrackFilter::SetEtaBins(Int_t nBins, Double_t *binning){
 360   //
 361   // User defined eta binning
 362   //
 363   fEtaBins = nBins;
 364   fEtaBinning = new Double_t[fEtaBins + 1];
 365   memcpy(fEtaBinning, binning, sizeof(Double_t) * nBins);
 366 }
 367
 368 //__________________________________________________________________
 369 void AliHFEtrackFilter::SetPhiBins(Int_t nBins, Double_t *binning){
 370   //
 371   // User defined phi binning
 372   //
 373   fPhiBins = nBins;
 374   fPhiBinning = new Double_t[fPhiBins + 1];
 375   memcpy(fPhiBinning, binning, sizeof(Double_t) * nBins);
 376 }
 377
 378 //__________________________________________________________________
 379 AliHFEcutStep *AliHFEtrackFilter::MakeCutStepRecKineITSTPC(){
 380   //
 381   // Make the cut step for Rec Kine
 382   // Cut step is already included in the filter
 383   //
 384   AliHFEcutStep *fCutStep = new AliHFEcutStep("RecKineITSTPC");
 385
 386   AliCFTrackQualityCuts *trackQuality = new AliCFTrackQualityCuts((Char_t *)"QualityRec", (Char_t *)"REC Track Quality Cuts");
 387   trackQuality->SetMinNClusterTPC(80);
 388   trackQuality->SetMaxChi2PerClusterTPC(3.5);
 389   trackQuality->SetStatus(AliESDtrack::kTPCrefit | AliESDtrack::kITSrefit);
 390   trackQuality->SetMaxCovDiagonalElements(2., 2., 0.5, 0.5, 2);
 391   fCutStep->AddCut(trackQuality);
 392
 393   AliHFEextraCuts *hfecuts = new AliHFEextraCuts("HFETPC","Extra cuts from the HFE group");
 394   hfecuts->SetClusterRatioTPC(0.6);
 395   fCutStep->AddCut(hfecuts);
 396
 397   AliCFTrackKineCuts *kineCuts = new AliCFTrackKineCuts((Char_t *)"RecKine", (Char_t *)"REC Kine Cuts");
 398   kineCuts->SetPtRange(0.1, 20);
 399   kineCuts->SetEtaRange(-0.8, 0.8);
 400   fCutStep->AddCut(kineCuts);
 401
 402   AddCutStep(fCutStep);
 403   return fCutStep;
 404 }
 405
 406 //__________________________________________________________________
 407 AliHFEcutStep *AliHFEtrackFilter::MakeCutStepPrimary(){
 408   //
 409   // Make cut on primaries
 410   // Cut step is already included in the filter
 411   //
 412   AliHFEcutStep *fCutStep = new AliHFEcutStep("Primary");
 413
 414   AliCFTrackIsPrimaryCuts *primaryCut = new AliCFTrackIsPrimaryCuts((Char_t *)"PrimaryCuts", (Char_t *)"REC Primary Cuts");
 415   primaryCut->SetMaxDCAToVertexXY(0.5);
 416   primaryCut->SetMaxDCAToVertexZ(5.);
 417   primaryCut->SetAcceptKinkDaughters(kFALSE);
 418   fCutStep->AddCut(primaryCut);
 419
 420   AddCutStep(fCutStep);
 421   return fCutStep;
 422 }
 423
 424 //__________________________________________________________________
 425 AliHFEcutStep *AliHFEtrackFilter::MakeCutStepHFEITS(){
 426   //
 427   // Add special ITS cuts
 428   // Cut step is already included in the filter
 429   //
 430   AliHFEcutStep *fCutStep = new AliHFEcutStep("HFEITS");
 431
 432   AliHFEextraCuts *hfecuts = new AliHFEextraCuts((Char_t *)"HFEPixelsCuts",(Char_t *)"Extra cuts from the HFE group");
 433   hfecuts->SetRequireITSpixel(AliHFEextraCuts::kFirst);
 434   //hfecuts->SetCheckITSstatus(kTRUE);
 435   fCutStep->AddCut(hfecuts);
 436
 437   AddCutStep(fCutStep);
 438   return fCutStep;
 439 }
 440
 441 //__________________________________________________________________
 442 AliHFEcutStep *AliHFEtrackFilter::MakeCutStepHFETRD(){
 443   //
 444   // Add special TRD cut
 445   // Cut step is already included in the filter
 446   //
 447   AliHFEcutStep *fCutStep = new AliHFEcutStep("HFETRD");
 448
 449   AliHFEextraCuts *hfecuts = new AliHFEextraCuts("HFETRDCuts","Extra cuts from the HFE group");
 450   hfecuts->SetMinTrackletsTRD(0);
 451   fCutStep->AddCut(hfecuts);
 452
 453   AddCutStep(fCutStep);
 454   return fCutStep;
 455 }
 456
 457 //__________________________________________________________________
 458 AliHFEcutStep *AliHFEtrackFilter::MakeMCSignalCuts(){
 459   //
 460   // Define MC Signal
 461   // Cut step is already included in the filter
 462   //
 463   fMCsignal = new AliHFEcutStep("MCSignal");
 464   AliCFParticleGenCuts *genCuts = new AliCFParticleGenCuts((Char_t *)"fCutsGenMC", (Char_t *)"Particle Generation Cuts");
 465   genCuts->SetRequireIsCharged();
 466   genCuts->SetRequireIsPrimary();
 467   genCuts->SetProdVtxRange2D();
 468   genCuts->SetProdVtxRangeX(0, 1);
 469   genCuts->SetProdVtxRangeY(0, 1);
 470   genCuts->SetRequirePdgCode(11, kTRUE);
 471   fMCsignal->AddCut(genCuts);
 472
 473   AliCFTrackKineCuts *kineMCcuts = new AliCFTrackKineCuts((Char_t *)"fCutsKineMC",(Char_t *)"MC Kine Cuts");
 474   kineMCcuts->SetPtRange(0.1, 20.);
 475   kineMCcuts->SetEtaRange(-0.8, 0.8);
 476   fMCsignal->AddCut(kineMCcuts);
 477
 478   return fMCsignal;
 479 }
 480
 481 //__________________________________________________________________
 482 void AliHFEtrackFilter::SetMC(AliMCEvent * const mc){
 483   //
 484   // Publish MC event to the single cut steps
 485   //
 486   fMC = mc;
 487   fMCsignal->SetMC(fMC);
 488   AliHFEcutStep *cs = NULL;
 489   for(Int_t icut = 0; icut < fCutSteps->GetEntriesFast(); icut++)
 490     if((cs = dynamic_cast<AliHFEcutStep *>(fCutSteps->UncheckedAt(icut)))) cs->SetMC(fMC);
 491 }
 492
 493 //__________________________________________________________________
 494 void AliHFEtrackFilter::SetRecEvent(AliVEvent *rec){
 495   //
 496   // Publish MC event to the single cut steps
 497   //
 498   AliHFEcutStep *cs = NULL;
 499   for(Int_t icut = 0; icut < fCutSteps->GetEntriesFast(); icut++)
 500     if((cs = dynamic_cast<AliHFEcutStep *>(fCutSteps->UncheckedAt(icut)))) cs->SetRecEvent(rec);
 501 }
 502
 503
 504 //__________________________________________________________________
 505 AliCFContainer *AliHFEtrackFilter::GetEfficiencyContainer(Int_t icont){
 506   //
 507   // return EfficiencyContainer
 508   //
 509   if(icont >= 4 || icont < 0) return NULL;
 510   return dynamic_cast<AliCFContainer *>(fEfficiencyContainers->At(icont));
 511 }