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