PWG3/hfe/AliAnalysisTaskDCA.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 // impact parameter resolution and pull study
  18 // for tracks which survivied the particle cuts
  19 //
  20 //
  21 // Authors:
  22 //  Hongyan Yang <hongyan@physi.uni-heidelberg.de>
  23 //  Carlo Bombonati <carlo.bombonati@cern.ch>
  24 //
  25 #include <TChain.h>
  26 #include <TFile.h>
  27 #include <TH1F.h>
  28 #include <TH1I.h>
  29 #include <TList.h>
  30 #include <TMath.h>
  31 #include <TObjArray.h>
  32 #include <TParticle.h>
  33 #include <TString.h>
  34
  35 #include <TCanvas.h>
  36
  37 #include "AliCFManager.h"
  38 #include "AliMCEvent.h"
  39 #include "AliESDInputHandler.h"
  40 #include "AliESDtrack.h"
  41 #include "AliAnalysisManager.h"
  42 #include "AliMCEventHandler.h"
  43 #include "AliMCParticle.h"
  44
  45 #include "AliHFEcuts.h"
  46 #include "AliHFEdca.h"
  47
  48 #include "AliAnalysisTaskDCA.h"
  49
  50
  51 //____________________________________________________________
  52 AliAnalysisTaskDCA::AliAnalysisTaskDCA():
  53   AliAnalysisTask("Impact Parameter Resolution and Pull Analysis", "")
  54   , fPlugins(0)
  55   , fESD(0x0)
  56   , fMC(0x0)
  57   , fCuts(0x0)
  58   , fCFM(0x0)
  59   , fDCA(0x0)
  60   , fPixelStatus(0x0)
  61   , fNclustersITS(0x0)
  62   , fNEvents(0x0)
  63   , fResidualList(0x0)
  64   , fPullList(0x0)
  65   , fOutput(0x0)
  66 {
  67   //
  68   // Dummy constructor
  69   //
  70   DefineInput(0, TChain::Class());
  71   DefineOutput(0, TH1I::Class());   // event
  72   DefineOutput(1, TList::Class());  // output
  73
  74   SetHasMCData();
  75
  76 }
  77
  78 //____________________________________________________________
  79 AliAnalysisTaskDCA::AliAnalysisTaskDCA(const char * name):
  80   AliAnalysisTask(name, "")
  81   , fPlugins(0)
  82   , fESD(0x0)
  83   , fMC(0x0)
  84   , fCuts(0x0)
  85   , fCFM(0x0)
  86   , fDCA(0x0)
  87   , fPixelStatus(0x0)
  88   , fNclustersITS(0x0)
  89   , fNEvents(0x0)
  90   , fResidualList(0x0)
  91   , fPullList(0x0)
  92   , fOutput(0x0)
  93 {
  94   //
  95   // Default constructor
  96   //
  97   DefineInput(0, TChain::Class());
  98   DefineOutput(0, TH1I::Class());
  99   DefineOutput(1, TList::Class());
 100
 101   SetHasMCData();
 102 }
 103
 104 //____________________________________________________________
 105 AliAnalysisTaskDCA::AliAnalysisTaskDCA(const AliAnalysisTaskDCA &ref):
 106   AliAnalysisTask(ref)
 107   , fPlugins(ref.fPlugins)
 108   , fESD(ref.fESD)
 109   , fMC(ref.fMC)
 110   , fCuts(ref.fCuts)
 111   , fCFM(ref.fCFM)
 112   , fDCA(ref.fDCA)
 113   , fPixelStatus(ref.fPixelStatus)
 114   , fNclustersITS(ref.fNclustersITS)
 115   , fNEvents(ref.fNEvents)
 116   , fResidualList(ref.fResidualList)
 117   , fPullList(ref.fPullList)
 118   , fOutput(ref.fOutput)
 119 {
 120   //
 121   // Copy Constructor
 122   //
 123 }
 124
 125 //____________________________________________________________
 126 AliAnalysisTaskDCA &AliAnalysisTaskDCA::operator=(const AliAnalysisTaskDCA &ref){
 127   //
 128   // Assignment operator
 129   //
 130   if(this == &ref) return *this;
 131   AliAnalysisTask::operator=(ref);
 132   fPlugins = ref.fPlugins;
 133   fESD = ref.fESD;
 134   fMC = ref.fMC;
 135   fCuts = ref.fCuts;
 136   fCFM = ref.fCFM;
 137   fDCA = ref.fDCA;
 138   fPixelStatus = ref.fPixelStatus;
 139   fNclustersITS = ref.fNclustersITS;
 140   fNEvents = ref.fNEvents;
 141   fOutput = ref.fOutput;
 142   fResidualList = ref.fResidualList;
 143   fPullList = ref.fPullList;
 144
 145   return *this;
 146 }
 147
 148 //____________________________________________________________
 149 AliAnalysisTaskDCA::~AliAnalysisTaskDCA(){
 150   //
 151   // Destructor
 152   //
 153
 154   if(fESD) delete fESD;
 155   if(fMC) delete fMC;
 156
 157   if(fCFM) delete fCFM;
 158
 159   if(fDCA) delete fDCA;
 160
 161   if(fOutput){
 162     fOutput->Clear();
 163     delete fOutput;
 164   }
 165
 166   if(fResidualList){
 167     fResidualList->Clear();
 168     delete fResidualList;
 169   }
 170
 171   if(fPullList){
 172     fPullList->Clear();
 173     delete fPullList;
 174   }
 175
 176   if(fNEvents) delete fNEvents;
 177 }
 178
 179 //____________________________________________________________
 180 void AliAnalysisTaskDCA::ConnectInputData(Option_t *){
 181   //
 182   // Connecting the input
 183
 184   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 185   if(!esdH){
 186     AliError("No ESD input handler");
 187     return;
 188   } else {
 189     fESD = esdH->GetEvent();
 190   }
 191   AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 192   if(!mcH){
 193     AliError("No MC truth handler");
 194     return;
 195   } else {
 196     fMC = mcH->MCEvent();
 197   }
 198 }
 199
 200 //____________________________________________________________
 201 void AliAnalysisTaskDCA::CreateOutputObjects(){
 202   // create output objects
 203   // fNEvents
 204   // residual and pull
 205
 206   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
 207   if(!fOutput) fOutput = new TList;
 208   // Initialize correction Framework and Cuts
 209   fCFM = new AliCFManager;
 210   MakeParticleContainer();
 211   // Temporary fix: Initialize particle cuts with 0x0
 212   for(Int_t istep = 0; istep < fCFM->GetParticleContainer()->GetNStep(); istep++)
 213     fCFM->SetParticleCutsList(istep, 0x0);
 214   if(!fCuts){
 215     AliWarning("Cuts not available. Default cuts will be used");
 216     fCuts = new AliHFEcuts;
 217     fCuts->CreateStandardCuts();
 218   }
 219
 220   fCuts->SetCutITSpixel(fPixelStatus);
 221   fCuts->Initialize(fCFM);
 222
 223
 224   // dca study----------------------------------
 225   if(!fDCA) fDCA = new AliHFEdca;
 226   if(!fResidualList) fResidualList = new TList();
 227   if(!fPullList) fPullList = new TList();
 228
 229   fDCA->CreateHistogramsResidual(fResidualList);
 230   fDCA->CreateHistogramsPull(fPullList);
 231
 232   // add output objects to the List
 233   fOutput->AddAt(fResidualList,0);
 234   fOutput->AddAt(fPullList,1);
 235
 236 }
 237
 238 //____________________________________________________________
 239 void AliAnalysisTaskDCA::Exec(Option_t *){
 240   //
 241   // Run the analysis
 242   //
 243
 244   AliDebug(3, "Processing ESD events");
 245
 246   if(!fESD){
 247     AliError("No ESD Event");
 248     return;
 249   }
 250   if(!fMC){
 251     AliError("No MC Event");
 252     return;
 253   }
 254   if(!fCuts){
 255     AliError("HFE cuts not available");
 256     return;
 257   }
 258
 259   //
 260   // Loop ESD
 261   //
 262   AliESDtrack *track = 0x0;
 263   fCFM->SetRecEventInfo(fESD);
 264   // event cut level
 265   if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fESD)) return;
 266
 267   for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
 268
 269     track = fESD->GetTrack(itrack);
 270
 271     // RecPrim: primary cuts
 272     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim, track)) continue;
 273     // RecKine: ITSTPC cuts
 274     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
 275     // HFEcuts: ITS layers cuts
 276     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsITS, track)) continue;
 277
 278     if(track->GetITSclusters(0)<=fNclustersITS) continue;  // require 6 hits on all pixel layers
 279
 280     if(GetPlugin(kImpactPar)) {
 281       //      Printf("analysis on impact parameter is ON");
 282       fDCA->FillHistograms(fESD, track, fMC);
 283     }
 284
 285   }
 286   fNEvents->Fill(1);
 287
 288   PostData(0, fNEvents);
 289   PostData(1, fOutput);
 290 }
 291
 292 //____________________________________________________________
 293 void AliAnalysisTaskDCA::Terminate(Option_t *){
 294   //
 295   // Terminate not implemented at the moment
 296   //
 297
 298   if(GetPlugin(kPostProcess)){
 299     fOutput = dynamic_cast<TList *>(GetOutputData(1));
 300     if(!fOutput){
 301       AliError("Results not available");
 302       return;
 303     }
 304     PostProcess();
 305   }
 306
 307 }
 308
 309
 310 //____________________________________________________________
 311 void AliAnalysisTaskDCA::Load(TString filename){
 312   // no need for postprocessing for the moment
 313   TFile *input = TFile::Open(filename.Data());
 314   if(!input || input->IsZombie()){
 315     AliError("Cannot read file");
 316     return;
 317   }
 318
 319   /*
 320   TH1 *htmp = dynamic_cast<TH1I *>(input->Get("nEvents"));
 321   if(htmp)
 322     fNEvents = dynamic_cast<TH1I *>(htmp->Clone());
 323   else
 324     AliError("Event Counter histogram not found");
 325   */
 326   input->Close();
 327   delete input;
 328
 329
 330 }
 331
 332 //____________________________________________________________
 333 void AliAnalysisTaskDCA::PostProcess(){
 334   // do post processing
 335   // should do fitting here for dca resolution
 336   // moved to an external macro to do the job
 337
 338   Load("impactPar.root");
 339   TCanvas *c1 = new TCanvas("c1", "number of analyzed events", 300, 400);
 340   fNEvents->Draw();
 341   c1->SaveAs("temp.png");
 342
 343 }
 344
 345
 346
 347
 348 //____________________________________________________________
 349 void AliAnalysisTaskDCA::PrintStatus() const {
 350
 351   //
 352   // Print Analysis status
 353   //
 354   printf("\n\tAnalysis Settings\n\t========================================\n");
 355   printf("\timpact parameter analysis is %s\n", GetPlugin(kImpactPar)?"ON":"OFF");
 356   printf("\tcuts: %s\n", (fCuts != NULL) ? "YES" : "NO");
 357   printf("\t ");
 358   printf("\n");
 359 }
 360
 361 //__________________________________________
 362 void AliAnalysisTaskDCA::SwitchOnPlugin(Int_t plug){
 363   //
 364   // Switch on Plugin
 365   // Available:
 366   //  - analyze impact parameter
 367   //  - Post Processing
 368
 369   switch(plug){
 370   case kImpactPar:
 371     SETBIT(fPlugins, plug);
 372     break;
 373   case kPostProcess:
 374     SETBIT(fPlugins, plug);
 375     break;
 376   default:
 377     AliError("Unknown Plugin");
 378   };
 379 }
 380
 381
 382 //____________________________________________________________
 383 void AliAnalysisTaskDCA::MakeParticleContainer(){
 384   //
 385   // Create the particle container (borrowed from AliAnalysisTaskHFE)
 386   //
 387   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
 388   const Double_t kPtmin = 0.1, kPtmax = 10.;
 389   const Double_t kEtamin = -0.9, kEtamax = 0.9;
 390   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
 391
 392   //arrays for the number of bins in each dimension
 393   Int_t iBin[kNvar];
 394   iBin[0] = 40; //bins in pt
 395   iBin[1] =  8; //bins in eta
 396   iBin[2] = 18; // bins in phi
 397
 398   //arrays for lower bounds :
 399   Double_t* binEdges[kNvar];
 400   for(Int_t ivar = 0; ivar < kNvar; ivar++)
 401     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
 402
 403   //values for bin lower bounds
 404   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)TMath::Power(10,TMath::Log10(kPtmin) + (TMath::Log10(kPtmax)-TMath::Log10(kPtmin))/iBin[0]*(Double_t)i);
 405   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
 406   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
 407
 408   //one "container" for MC
 409   AliCFContainer* container = new AliCFContainer("container","container for tracks", (AliHFEcuts::kNcutStepsTrack + 1 + 2*(AliHFEcuts::kNcutStepsESDtrack + 1)), kNvar, iBin);
 410
 411   //setting the bin limits
 412   for(Int_t ivar = 0; ivar < kNvar; ivar++)
 413     container -> SetBinLimits(ivar, binEdges[ivar]);
 414   fCFM->SetParticleContainer(container);
 415
 416   //create correlation matrix for unfolding
 417   Int_t thnDim[2*kNvar];
 418   for (int k=0; k<kNvar; k++) {
 419     //first half  : reconstructed
 420     //second half : MC
 421     thnDim[k]      = iBin[k];
 422     thnDim[k+kNvar] = iBin[k];
 423   }
 424
 425
 426   // add more containers for correction purpose
 427
 428
 429 }