PWG3/hfe/AliAnalysisElectronTask.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  * The analysis task:
  17  * Filling an AliCFContainer with the quantities pt, eta and phi
  18  * for tracks which survivied the particle cuts (MC resp. ESD tracks)
  19  * Track selection is done using the AliHFE package
  20  *
  21  * Author:
  22  *  Markus Fasel <M.Fasel@gsi.de>
  23  */
  24 #include <TAxis.h>
  25 #include <TCanvas.h>
  26 #include <TChain.h>
  27 #include <TH1F.h>
  28 #include <TH1I.h>
  29 #include <TH2F.h>
  30 #include <TIterator.h>
  31 #include <TList.h>
  32 #include <TLegend.h>
  33 #include <TMath.h>
  34 #include <TObjArray.h>
  35 #include <TParticle.h>
  36 #include <TProfile.h>
  37 #include <TTree.h>
  38
  39 #include "AliCFContainer.h"
  40 #include "AliCFManager.h"
  41 #include "AliESDEvent.h"
  42 #include "AliESDInputHandler.h"
  43 #include "AliESDtrack.h"
  44 #include "AliESDtrackCuts.h"
  45 #include "AliLog.h"
  46 #include "AliAnalysisManager.h"
  47 #include "AliMCEvent.h"
  48 #include "AliMCEventHandler.h"
  49 #include "AliMCParticle.h"
  50 #include "AliPID.h"
  51
  52 #include "AliHFEpid.h"
  53 #include "AliHFEcuts.h"
  54 #include "AliAnalysisElectronTask.h"
  55
  56 //____________________________________________________________
  57 AliAnalysisElectronTask::AliAnalysisElectronTask():
  58         AliAnalysisTask("PID efficiency Analysis", "")
  59         , fESD(0x0)
  60         , fMC(0x0)
  61         , fCFM(0x0)
  62         , fPID(0x0)
  63   , fCuts(0x0)
  64         , fNEvents(0x0)
  65         , fQA(0x0)
  66 {
  67         DefineInput(0, TChain::Class());
  68         DefineOutput(0, TH1I::Class());
  69         DefineOutput(1, AliCFContainer::Class());
  70         DefineOutput(2, TList::Class());
  71
  72   // Initialize cuts
  73   fCuts = new AliHFEcuts;
  74   fPID = new AliHFEpid;
  75 }
  76
  77 //____________________________________________________________
  78 AliAnalysisElectronTask::~AliAnalysisElectronTask(){
  79         if(fESD) delete fESD;
  80         if(fMC) delete fMC;
  81         if(fPID) delete fPID;
  82         if(fQA) delete fQA;
  83   if(fCuts) delete fCuts;
  84         if(fNEvents) delete fNEvents;
  85         fQA = 0x0; fMC = 0x0; fESD = 0x0; fCFM = 0x0; fPID = 0x0; fCuts = 0x0;
  86 }
  87
  88 //____________________________________________________________
  89 void AliAnalysisElectronTask::ConnectInputData(Option_t *){
  90         TTree *esdchain = dynamic_cast<TChain *>(GetInputData(0));
  91         if(!esdchain){
  92                 AliError("ESD chain empty");
  93                 return;
  94         } else {
  95                 esdchain->SetBranchStatus("Tracks", 1);
  96         }
  97         AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
  98         if(!esdH){
  99                 AliError("No ESD input handler");
 100                 return;
 101         } else {
 102                 fESD = esdH->GetEvent();
 103         }
 104         AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 105         if(!mcH){
 106                 AliError("No MC truth handler");
 107                 return;
 108         } else {
 109                 fMC = mcH->MCEvent();
 110         }
 111 }
 112
 113 //____________________________________________________________
 114 void AliAnalysisElectronTask::CreateOutputObjects(){
 115         fNEvents = new TH1I("nEvents", "Number of Events in the Analysis", 2, 0, 2); // Number of Events neccessary for the analysis and not a QA histogram
 116         // First Step: TRD alone
 117         if(!fQA) fQA = new TList;
 118         fQA->AddAt(new TProfile("conr", "Electron PID contamination", 20, 0, 20), 0);
 119         fQA->AddAt(new TH1F("alpha_rec", "Alpha from reconstructed tracks with TRD hits", 36, -TMath::Pi(), TMath::Pi()), 1);
 120         fQA->AddAt(new TH1F("alpha_sim", "Alpha from simulated electron tracks", 36, -TMath::Pi(), TMath::Pi()), 2);
 121         fQA->AddAt(new TH1F("nElectron", "Number of electrons", 100, 0, 100), 3);
 122         fQA->AddAt(new TProfile("pidquality", "TRD PID quality as function of momentum", 20, 0, 20), 4);
 123         fQA->AddAt(new TProfile("ntrdclusters", "Number of TRD clusters as function of momentum", 20, 0, 20), 5);
 124         fQA->AddAt(new TH1F("chi2TRD","#chi2 per TRD cluster", 20, 0, 20), 6);
 125
 126   // Initialize correction Framework and Cuts
 127   fCFM = new AliCFManager;
 128   MakeParticleContainer();
 129   // Temporary fix: Initialize particle cuts with 0x0
 130   for(Int_t istep = 0; istep < fCFM->GetParticleContainer()->GetNStep(); istep++)
 131     fCFM->SetParticleCutsList(istep, 0x0);
 132   if(IsQAOn()){
 133     fCuts->SetDebugMode();
 134     fQA->AddAt(fCuts->GetQAhistograms(), 7);
 135   }
 136   fCuts->CreateStandardCuts();
 137   fCuts->Initialize(fCFM);
 138
 139   // Initialize PID
 140   if(IsQAOn()){
 141     fPID->SetQAOn();
 142     fQA->AddAt(fPID->GetQAhistograms(), 8);
 143   }
 144   fPID->SetHasMCData(kTRUE);
 145   fPID->InitializePID("TRD");
 146 }
 147
 148 //____________________________________________________________
 149 void AliAnalysisElectronTask::Exec(Option_t *){
 150         //
 151         // Run the analysis
 152         //
 153         if(!fESD){
 154                 AliError("No ESD Event");
 155                 return;
 156         }
 157         if(!fMC){
 158                 AliError("No MC Event");
 159                 return;
 160         }
 161         fCFM->SetEventInfo(fMC);
 162   fPID->SetMCEvent(fMC);
 163
 164         //fCFM->CheckEventCuts(AliCFManager::kEvtGenCuts, fMC);
 165
 166         Double_t pt = 0;
 167         Double_t container[3];
 168
 169         // Loop over the Monte Carlo tracks to see whether we have overlooked any track
 170         AliMCParticle *mctrack = 0x0;
 171         Int_t nElectrons = 0;
 172         for(Int_t imc = fMC->GetNumberOfTracks(); imc--;){
 173                 mctrack = fMC->GetTrack(imc);
 174                 container[0] = mctrack->Pt();
 175                 container[1] = mctrack->Eta();
 176     container[2] = mctrack->Phi();
 177
 178                 if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) continue;
 179                 fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCGenerated);
 180                 (dynamic_cast<TH1F *>(fQA->At(2)))->Fill(mctrack->Phi() - TMath::Pi());
 181                 if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCInAcceptance, mctrack)) continue;
 182                 // find the label in the vector
 183                 fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCInAcceptance);
 184                 nElectrons++;
 185         }
 186         (dynamic_cast<TH1F *>(fQA->At(3)))->Fill(nElectrons);
 187
 188         // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
 189         AliESDtrack *track = 0x0;
 190         for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
 191                 track = fESD->GetTrack(itrack);
 192                 container[0] = track->Pt();
 193                 container[1] = track->Eta();
 194     container[2] = track->Phi();
 195                 if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKine, track)) continue;
 196                 fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecKine);
 197
 198                 // Check TRD criterions (outside the correction framework)
 199                 if(track->GetTRDncls()){
 200                         (dynamic_cast<TH1F *>(fQA->At(6)))->Fill(track->GetTRDchi2()/track->GetTRDncls());
 201                         (dynamic_cast<TH1F *>(fQA->At(1)))->Fill(track->GetAlpha());    // Check the acceptance without tight cuts
 202                         (dynamic_cast<TProfile *>(fQA->At(4)))->Fill(container[0], track->GetTRDpidQuality());
 203                         (dynamic_cast<TProfile *>(fQA->At(5)))->Fill(container[0], track->GetTRDncls());
 204                 }
 205                 if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim, track)) continue;
 206     fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecPrim);
 207     if(fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcuts, track)) continue;
 208     fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcuts);
 209     // track accepted, do PID
 210                 if(!fPID->IsSelected(track)) continue;
 211         fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcuts + 1);
 212         }
 213         fNEvents->Fill(1);
 214
 215         // Done!!!
 216         PostData(0, fNEvents);
 217         PostData(1, fCFM->GetParticleContainer());
 218         PostData(2, fQA);
 219 }
 220
 221 //____________________________________________________________
 222 void AliAnalysisElectronTask::Terminate(Option_t *){
 223         //
 224         // Terminate not implemented at the moment
 225         //
 226 }
 227
 228 //____________________________________________________________
 229 void AliAnalysisElectronTask::MakeParticleContainer(){
 230   //
 231   // Create the particle container for the correction framework manager and
 232   // link it
 233   //
 234   const Int_t nvar   = 3 ; //number of variables on the grid:pt,eta, phi
 235   const Double_t ptmin = 0., ptmax = 10.;
 236   const Double_t etamin = -0.9, etamax = 0.9;
 237   const Double_t phimin = 0., phimax = 2. * TMath::Pi();
 238
 239
 240   //arrays for the number of bins in each dimension
 241   Int_t iBin[nvar];
 242   iBin[0] = 20; //bins in pt
 243   iBin[1] =  8; //bins in eta
 244   iBin[2] = 18; // bins in phi
 245
 246   //arrays for lower bounds :
 247   Double_t *binLim1 = new Double_t[iBin[0] + 1];
 248   Double_t *binLim2 = new Double_t[iBin[1] + 1];
 249   Double_t *binLim3 = new Double_t[iBin[2] + 1];
 250
 251   //values for bin lower bounds
 252   for(Int_t i=0; i<=iBin[0]; i++) binLim1[i]=(Double_t)ptmin + (ptmax-ptmin)/iBin[0]*(Double_t)i;
 253   for(Int_t i=0; i<=iBin[1]; i++) binLim2[i]=(Double_t)etamin  + (etamax-etamin)/iBin[1]*(Double_t)i;
 254   for(Int_t i=0; i<=iBin[2]; i++) binLim3[i]=(Double_t)phimin  + (phimax-phimin)/iBin[2]*(Double_t)i;
 255
 256   //one "container" for MC
 257   AliCFContainer* container = new AliCFContainer("container","container for tracks", AliHFEcuts::kNcutSteps + 1, nvar, iBin);
 258   //setting the bin limits
 259   container -> SetBinLimits(0,binLim1);
 260   container -> SetBinLimits(1,binLim2);
 261   container -> SetBinLimits(2,binLim3);
 262   fCFM->SetParticleContainer(container);
 263 }