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