PWG3/hfe/AliHFEefficiency.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 // Task for Efficiency studies
  17 // Used for testing classes AliHFEcontainer and AliHFEfilter
  18 // Creates Efficiency Histograms
  19 //
  20 // Author:
  21 //   Markus Fasel <M.Fasel@gsi.de>
  22 //
  23 #include <TAxis.h>
  24 #include <TCanvas.h>
  25 #include <TFile.h>
  26 #include <TGraphErrors.h>
  27 #include <TH1F.h>
  28 #include <TH2.h>
  29 #include <TIterator.h>
  30 #include <TLegend.h>
  31 #include <TObjArray.h>
  32 #include <TPad.h>
  33
  34 #include "AliCFAcceptanceCuts.h"
  35 #include "AliCFContainer.h"
  36 #include "AliCFEffGrid.h"
  37 #include "AliESDEvent.h"
  38 #include "AliHFEcollection.h"
  39 #include "AliHFEcontainer.h"
  40 #include "AliHFEcutStep.h"
  41 #include "AliHFEefficiency.h"
  42 #include "AliHFEextraCuts.h"
  43 #include "AliHFEtrackFilter.h"
  44 #include "AliHFEtools.h"
  45 #include "AliMCEvent.h"
  46
  47 ClassImp(AliHFEefficiency)
  48
  49 AliHFEefficiency::AliHFEefficiency():
  50   AliAnalysisTaskSE(),
  51   fFilter(NULL),
  52   fMCcut(NULL),
  53   fAcceptanceCuts(NULL),
  54   fEfficiency(NULL),
  55   fOutput(NULL),
  56   fCutTRD(kFALSE)
  57 {
  58   //
  59   // Default constructor
  60   //
  61 }
  62
  63 AliHFEefficiency::AliHFEefficiency(const Char_t *name):
  64   AliAnalysisTaskSE(name),
  65   fFilter(NULL),
  66   fMCcut(NULL),
  67   fAcceptanceCuts(NULL),
  68   fEfficiency(NULL),
  69   fOutput(NULL),
  70   fCutTRD(kFALSE)
  71 {
  72   //
  73   // Main Constructor
  74   //
  75   DefineOutput(1, AliHFEcontainer::Class());
  76   DefineOutput(2, TList::Class());
  77 }
  78
  79 AliHFEefficiency::~AliHFEefficiency(){
  80   //
  81   // Destructor
  82   //
  83   if(fFilter) delete fFilter;
  84   if(fEfficiency) delete fEfficiency;
  85   if(fAcceptanceCuts) delete fAcceptanceCuts;
  86   if(fOutput) delete fOutput;
  87 }
  88
  89 void AliHFEefficiency::UserCreateOutputObjects(){
  90   //
  91   // Create the output objects
  92   //
  93   fEfficiency = new AliHFEcontainer("Efficiency");
  94   fEfficiency->SetNumberOfVariables(4);
  95   // InitializeVariables
  96   fEfficiency->SetVariableName(0, "pt");
  97   fEfficiency->MakeLogarithmicBinning(0, 40, 0.1, 10);
  98   fEfficiency->SetVariableName(1, "Eta");
  99   fEfficiency->MakeLinearBinning(1, 20, -0.9, 0.9);
 100   fEfficiency->SetVariableName(2, "phi");
 101   fEfficiency->MakeLinearBinning(2 , 18, 0, 2 * TMath::Pi());
 102   fEfficiency->SetVariableName(3, "Charge");
 103   fEfficiency->MakeLinearBinning(3, 2, -1.1, 1.1);
 104   // New container
 105   fEfficiency->CreateContainer("MCFilter", "Efficiency for MC Studies", 2);
 106   AliCFContainer *cont = fEfficiency->GetCFContainer("MCFilter");
 107   cont->SetStepTitle(0, "MC Signal");
 108   cont->SetStepTitle(1, "MC Signal in acceptance");
 109
 110   fFilter = new AliHFEtrackFilter("testfilter");
 111   fFilter->GenerateCutSteps();
 112   fFilter->InitCF(fEfficiency);
 113   fMCcut = fFilter->GetMCSignalCuts();
 114   if(fCutTRD){
 115     AliHFEcutStep *hfeTRD = fFilter->GetCutStep("HFETRD");
 116     AliHFEextraCuts *hfecuts = dynamic_cast<AliHFEextraCuts *>(hfeTRD->GetCut("HFETRDCuts"));
 117     hfecuts->SetMinTrackletsTRD(4);
 118   }
 119   AliHFEcutStep *hfeITS = fFilter->GetCutStep("HFEITS");
 120   AliHFEextraCuts *hfeitscuts = dynamic_cast<AliHFEextraCuts *>(hfeITS->GetCut("HFEPixelsCuts"));
 121   hfeitscuts->SetRequireITSpixel(AliHFEextraCuts::kFirst);
 122
 123   fAcceptanceCuts = new AliCFAcceptanceCuts("Acceptance", "MC Acceptance Cuts");
 124   fAcceptanceCuts->SetMinNHitITS(3);
 125   fAcceptanceCuts->SetMinNHitTPC(2);
 126   fAcceptanceCuts->SetMinNHitTRD(0);
 127
 128   // create additional histos for testing purpose
 129   fOutput = new AliHFEcollection("histos", "QA histograms");
 130   fOutput->CreateTH1F("hNtracks","Number of Tracks per Event", 100, 0, 100);
 131   fOutput->CreateTH1F("hPt","Pt of the tracks", 40, 0.1, 10);
 132   fOutput->BinLogAxis("hPt", 0);
 133   fOutput->CreateTH1F("kinkIndex", "Kink Index", 400, -200, 200);
 134   fOutput->CreateTH1F("itspixel", "ITS PIXEL", 2, 0, 2);
 135   fOutput->CreateTH1F("mcmother", "Mother PDG", 1000, -500, 500);
 136   fOutput->CreateTH2F("ptres", "Pt Resolution", 40, 0.1, 10, 200, -1.5, 1.5);
 137   fOutput->BinLogAxis("ptres", 0);
 138 }
 139
 140 void AliHFEefficiency::UserExec(Option_t *){
 141   //
 142   // Event Loop
 143   // Filter track, fill Efficiency container
 144   //
 145   fEfficiency->NewEvent();
 146   fFilter->SetRecEvent(fInputEvent);
 147   if(fMCEvent){
 148     fFilter->SetMC(fMCEvent);
 149     FilterMC();
 150   }
 151   fFilter->FilterTracks(dynamic_cast<AliESDEvent *>(fInputEvent));
 152   TObjArray *tracks = fFilter->GetFilteredTracks();
 153   TIterator *iter = tracks->MakeIterator();
 154   fOutput->Fill("hNtracks", tracks->GetEntriesFast());
 155   AliESDtrack *track = NULL;
 156   while((track = dynamic_cast<AliESDtrack *>(iter->Next()))){
 157     fOutput->Fill("hPt", track->Pt());
 158     fOutput->Fill("kinkIndex", track->GetKinkIndex(0));
 159     if(TESTBIT(track->GetITSClusterMap(), 0))
 160       fOutput->Fill("itspixel",1);
 161     else
 162       fOutput->Fill("itspixel",0);
 163     AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel())));
 164     AliMCParticle *mother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(mctrack->Particle()->GetFirstMother()));
 165     fOutput->Fill("mcmother", mother->Particle()->GetPdgCode());
 166     fOutput->Fill("ptres", mctrack->Pt(), (track->Pt() - mctrack->Pt())/mctrack->Pt());
 167   }
 168   delete iter;
 169   fFilter->Flush();
 170   PostData(1, fEfficiency);
 171   PostData(2, fOutput->GetList());
 172 }
 173
 174 void AliHFEefficiency::Terminate(Option_t *){
 175   //
 176   // Evaluate Results
 177   //
 178   fEfficiency = dynamic_cast<AliHFEcontainer *>(GetOutputData(1));
 179   if(!fEfficiency) return;
 180   PostProcess();
 181
 182   TList *l = dynamic_cast<TList *>(GetOutputData(2));
 183   DrawPtResolution(l);
 184 }
 185
 186 void AliHFEefficiency::FilterMC(){
 187   //
 188   // Cut MC tracks
 189   //
 190   Double_t cont[4] = {0., 0., 0., 0.};
 191   AliMCParticle *track = NULL;
 192   for(Int_t itrack = 0; itrack < fMCEvent->GetNumberOfTracks(); itrack++){
 193     track = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(itrack));
 194     if(!fMCcut->IsSelected(track)) continue;
 195     cont[0] = track->Pt();
 196     cont[1] = track->Eta();
 197     cont[2] = track->Phi();
 198     cont[3] = track->Charge()/3;
 199     //AliMCParticle *mother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(track->Particle()->GetFirstMother()));
 200     //if(TMath::Abs(mother->Particle()->GetPdgCode()) != 443) continue;
 201     fEfficiency->FillCFContainer("MCFilter", 0, cont);
 202     if(!fAcceptanceCuts->IsSelected(track)) continue;
 203     fEfficiency->FillCFContainer("MCFilter", 1, cont);
 204   }
 205 }
 206
 207 void AliHFEefficiency::Load(const char* filename){
 208   TFile *input = TFile::Open(filename);
 209   AliHFEcontainer *cin = dynamic_cast<AliHFEcontainer *>(input->Get("Efficiency"));
 210   fEfficiency = dynamic_cast<AliHFEcontainer *>(cin->Clone());
 211   input->Close();
 212   delete input;
 213 }
 214
 215 void AliHFEefficiency::PostProcess(){
 216   //
 217   // Calculate Efficiencies
 218   //
 219   fEfficiency->Print();
 220
 221   // Prepare containers
 222   AliCFContainer *crec  = fEfficiency->MakeMergedCFContainer("merged", "Merged Container", "MCFilter:testfilter_container_signal");
 223   AliCFContainer *cmc   = fEfficiency->MakeMergedCFContainer("merged", "Merged Container", "MCFilter:testfilter_container_signalMC");
 224   AliCFEffGrid *gridRec = new AliCFEffGrid("effrec", "Efficiency Calculation using rec vars", *crec);
 225   AliCFEffGrid *gridMC  = new AliCFEffGrid("effmc", "Efficiency Calculation using mc vars", *cmc);
 226
 227   TCanvas *ctmp = NULL;
 228   for(Int_t ivar = 0; ivar < crec->GetNVar(); ivar++){
 229     ctmp = new TCanvas(Form("cEffSignal%s", crec->GetVarTitle(ivar)), Form("Signal Efficiency (%s)", crec->GetVarTitle(ivar)));
 230     ctmp->cd();
 231     DrawSignalEfficiency(gridRec, crec, ivar);
 232
 233     ctmp = new TCanvas(Form("cCutEff%s", crec->GetVarTitle(ivar)), Form("Cut Step Efficiency (%s)", crec->GetVarTitle(ivar)));
 234     ctmp->cd();
 235     DrawCutEfficiency(gridRec, crec, ivar);
 236
 237     ctmp = new TCanvas(Form("cEffSignalMC%s", crec->GetVarTitle(ivar)), Form("Signal Efficiency (%s, MC)", crec->GetVarTitle(ivar)));
 238     ctmp->cd();
 239     DrawSignalEfficiency(gridMC, cmc, ivar);
 240
 241     ctmp = new TCanvas(Form("cCutEffMC%s", crec->GetVarTitle(ivar)), Form("Cut Step Efficiency (%s, MC)", crec->GetVarTitle(ivar)));
 242     ctmp->cd();
 243     DrawCutEfficiency(gridMC, cmc, ivar);
 244   }
 245
 246   CalculatePTsmearing();
 247
 248   delete gridRec;
 249   delete gridMC;
 250   delete crec;
 251   delete cmc;
 252 }
 253
 254 void AliHFEefficiency::DrawSignalEfficiency(AliCFEffGrid *eff, AliCFContainer *cont, Int_t var){
 255   //
 256   // Efficiency of a cut step with respect to the signal
 257   //
 258   TH1 *hEff = NULL;
 259   Bool_t kFirst = kTRUE;
 260   TLegend *leg = new TLegend(0.5, 0.7, 0.89, 0.89);
 261   for(Int_t istep = 1; istep < cont->GetNStep(); istep++){
 262     eff->CalculateEfficiency(istep, 0);
 263     hEff = eff->Project(var);
 264     hEff->SetTitle(Form("Signal Efficiency (%s)", cont->GetVarTitle(var)));
 265     hEff->GetYaxis()->SetRangeUser(0., 1.6);
 266     hEff->GetYaxis()->SetTitle("Efficiency");
 267     hEff->SetStats(kFALSE);
 268     hEff->SetMarkerStyle(20+istep);
 269     hEff->SetMarkerColor(kBlue - 5 + istep);
 270     hEff->Draw(kFirst ? "ep" : "epsame");
 271     kFirst = kFALSE;
 272     leg->AddEntry(hEff, cont->GetStepTitle(istep), "p");
 273   }
 274   leg->SetBorderSize(0);
 275   leg->SetFillStyle(0);
 276   leg->Draw();
 277   gPad->Update();
 278 }
 279
 280 void AliHFEefficiency::DrawCutEfficiency(AliCFEffGrid *eff, AliCFContainer *cont, Int_t var){
 281   //
 282   // Efficiency of a cut step with respect to the previous one
 283   //
 284   Bool_t kFirst = kTRUE;
 285   TH1 *hEff = NULL;
 286   TLegend *leg = new TLegend(0.5, 0.7, 0.89, 0.89);
 287   TString effTitle;
 288   Int_t start = 0, length = 0;
 289   for(Int_t istep = 1; istep < cont->GetNStep(); istep++){
 290     eff->CalculateEfficiency(istep, istep -1);
 291     hEff = eff->Project(var);
 292     effTitle = hEff->GetTitle();
 293     start = effTitle.Index("projected"); length = effTitle.Length() - start;
 294     effTitle.Remove(start - 1, length + 1);
 295     effTitle.ReplaceAll("Efficiency: ", "");
 296     hEff->SetTitle(Form("Cut Efficiency (%s)", cont->GetVarTitle(var)));
 297     hEff->GetYaxis()->SetRangeUser(0., 1.6);
 298     hEff->GetYaxis()->SetTitle("Efficiency");
 299     hEff->SetStats(kFALSE);
 300     hEff->SetMarkerStyle(20+istep);
 301     hEff->SetMarkerColor(kBlue - 5 + istep);
 302     hEff->Draw(kFirst ? "ep" : "epsame");
 303     kFirst = kFALSE;
 304     leg->AddEntry(hEff, effTitle.Data(), "p");
 305   }
 306   leg->SetBorderSize(0);
 307   leg->SetFillStyle(0);
 308   leg->Draw();
 309   gPad->Update();
 310 }
 311
 312 void AliHFEefficiency::CalculatePTsmearing(){
 313   //
 314   // Efficiency of a cut step with respect to the same cut step filled with MC
 315   // Information
 316   //
 317   AliCFContainer *cont = fEfficiency->MakeMergedCFContainer("merged", "Merged Container", "testfilter_container_signal:testfilter_container_signalMC");
 318   AliCFEffGrid *grid = new AliCFEffGrid("effc", "Efficiency Calculation", *cont);
 319   Bool_t kFirst = kTRUE;
 320   TH1 *hEff = NULL;
 321   TCanvas *c2 = new TCanvas("cSmear", "pT modification");
 322   c2->cd();
 323   TLegend *leg1 = new TLegend(0.5, 0.7, 0.89, 0.89);
 324   Int_t nStep = cont->GetNStep()/2;
 325   for(Int_t istep = 0; istep < nStep; istep++){
 326     grid->CalculateEfficiency(istep, istep + nStep);
 327     hEff = grid->Project(0);
 328     hEff->SetStats(kFALSE);
 329     hEff->SetMarkerStyle(20+istep);
 330     hEff->SetMarkerColor(kBlue - 5 + istep);
 331     hEff->Draw(kFirst ? "ep" : "epsame");
 332     kFirst = kFALSE;
 333     leg1->AddEntry(hEff, cont->GetStepTitle(istep), "p");
 334   }
 335   leg1->Draw();
 336   gPad->Update();
 337
 338   delete cont;
 339   delete grid;
 340 }
 341
 342 void AliHFEefficiency::DrawPtResolution(TList *l){
 343   //
 344   // Draw pt resolution
 345   //
 346   TH2 *h2 = dynamic_cast<TH2 *>(l->FindObject("ptres"));
 347   TGraphErrors *mean = new TGraphErrors(h2->GetNbinsX());
 348   TGraphErrors *rms = new TGraphErrors(h2->GetNbinsX());
 349
 350   Double_t pt = 0., pterr = 0., mymean = 0., myrms = 0., mymeanerr = 0., myrmserr = 0.;
 351   TH1 *htmp = NULL;
 352   for(Int_t imom = 1; imom <= h2->GetXaxis()->GetLast(); imom++){
 353     pt = h2->GetXaxis()->GetBinCenter(imom);
 354     pterr = h2->GetXaxis()->GetBinWidth(imom)/2.;
 355     htmp = h2->ProjectionY("py", imom, imom);
 356     mymean = htmp->GetMean();
 357     myrms = htmp->GetRMS();
 358     mymeanerr = htmp->GetMeanError();
 359     myrmserr = htmp->GetRMSError();
 360     delete htmp;
 361     mean->SetPoint(imom -1, pt, mymean);
 362     rms->SetPoint(imom -1, pt, myrms);
 363     mean->SetPointError(imom-1,pterr,mymeanerr);
 364     rms->SetPointError(imom-1,pterr,myrmserr);
 365   }
 366
 367   TCanvas *c1 = new TCanvas("cPtShift", "pT Shift");
 368   c1->cd();
 369   mean->SetMarkerStyle(22);
 370   mean->SetMarkerColor(kBlue);
 371   mean->Draw("ape");
 372   rms->SetMarkerStyle(22);
 373   rms->SetMarkerColor(kBlack);
 374   rms->Draw("pesame");
 375   TLegend *leg = new TLegend(0.5, 0.7, 0.89, 0.89);
 376   leg->AddEntry(mean, "Mean", "lp");
 377   leg->AddEntry(rms, "RMS", "lp");
 378   leg->Draw();
 379   gPad->Update();
 380 }
 381