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