PWGHF/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 <Riostream.h>
  26 #include <TChain.h>
  27 #include <TFile.h>
  28 #include <TH1F.h>
  29 #include <TH1I.h>
  30 #include <TList.h>
  31 #include <TMath.h>
  32 #include <TObjArray.h>
  33 #include <TParticle.h>
  34 #include <TString.h>
  35
  36 #include <TCanvas.h>
  37
  38 #include "AliAnalysisManager.h"
  39
  40 #include "AliCFManager.h"
  41
  42 #include "AliAODEvent.h"
  43 #include "AliESDInputHandler.h"
  44 #include "AliESDtrack.h"
  45 #include "AliVertexerTracks.h"
  46 #include "AliESDVertex.h"
  47
  48 #include "AliMCEventHandler.h"
  49 #include "AliMCEvent.h"
  50 #include "AliMCParticle.h"
  51 #include "AliPIDResponse.h"
  52
  53 #include "AliHFEpid.h"
  54 #include "AliHFEcuts.h"
  55 #include "AliHFEdca.h"
  56 #include "AliHFEtools.h"
  57
  58 #include "AliAnalysisTaskDCA.h"
  59
  60
  61 //____________________________________________________________
  62 AliAnalysisTaskDCA::AliAnalysisTaskDCA():
  63   AliAnalysisTaskSE("Impact Parameter Resolution and Pull Analysis")
  64   , fPlugins(0)
  65   , fCuts(0x0)
  66   , fHFEpid(0x0)
  67   , fPIDdetectors("")
  68   , fPIDstrategy(0)
  69   , fCFM(0x0)
  70   , fDCA(0x0)
  71   , fNclustersITS(0x0)
  72   , fMinNprimVtxContrbutor(0x0)
  73   , fNEvents(0x0)
  74   , fResidualList(0x0)
  75   , fPullList(0x0)
  76   , fDcaList(0x0)
  77   , fKfDcaList(0x0)
  78   , fMcVertexList(0x0)
  79   , fDataDcaList(0x0)
  80   , fDataVertexList(0x0)
  81   , fDataPullList(0x0)
  82   , fMcPidList(0x0)
  83   , fDataPidList(0x0)
  84   , fHfeDcaList(0x0)
  85   , fHfeDataDcaList(0x0)
  86   , fOutput(0x0)
  87 {
  88   //
  89   // Dummy constructor
  90   //
  91 }
  92
  93 //____________________________________________________________
  94 AliAnalysisTaskDCA::AliAnalysisTaskDCA(const char * name):
  95   AliAnalysisTaskSE(name)
  96   , fPlugins(0)
  97   , fCuts(0x0)
  98   , fHFEpid(0x0)
  99   , fPIDdetectors("")
 100   , fPIDstrategy(0)
 101   , fCFM(0x0)
 102   , fDCA(0x0)
 103   , fNclustersITS(0x0)
 104   , fMinNprimVtxContrbutor(0x0)
 105   , fNEvents(0x0)
 106   , fResidualList(0x0)
 107   , fPullList(0x0)
 108   , fDcaList(0x0)
 109   , fKfDcaList(0x0)
 110   , fMcVertexList(0x0)
 111   , fDataDcaList(0x0)
 112   , fDataVertexList(0x0)
 113   , fDataPullList(0x0)
 114   , fMcPidList(0x0)
 115   , fDataPidList(0x0)
 116   , fHfeDcaList(0x0)
 117   , fHfeDataDcaList(0x0)
 118   , fOutput(0x0)
 119 {
 120   //
 121   // Default constructor
 122   //
 123   DefineInput(0, TChain::Class());
 124   DefineOutput(1, TH1I::Class());
 125   DefineOutput(2, TList::Class());
 126
 127
 128   //-CUTS SETTING-//
 129   Int_t nMinTPCcluster = 100;
 130   Float_t maxDcaXY = 0.5;
 131   Float_t maxDcaZ = 1.0;
 132   //--------------//
 133   AliHFEcuts *hfecuts = new AliHFEcuts;
 134   hfecuts->CreateStandardCuts();
 135   hfecuts->SetMinNClustersTPC(nMinTPCcluster);
 136   hfecuts->SetCutITSpixel(AliHFEextraCuts::kFirst);
 137   hfecuts->SetCheckITSLayerStatus(kFALSE);
 138   hfecuts->SetMaxImpactParam(maxDcaXY, maxDcaZ);
 139   SetHFECuts(hfecuts);
 140
 141   fHFEpid = new AliHFEpid("PIDforDCAanalysis");
 142
 143 }
 144
 145 //____________________________________________________________
 146 AliAnalysisTaskDCA::AliAnalysisTaskDCA(const AliAnalysisTaskDCA &ref):
 147   AliAnalysisTaskSE(ref)
 148   , fPlugins(ref.fPlugins)
 149   , fCuts(ref.fCuts)
 150   , fHFEpid(ref.fHFEpid)
 151   , fPIDdetectors(ref.fPIDdetectors)
 152   , fPIDstrategy(ref.fPIDstrategy)
 153   , fCFM(ref.fCFM)
 154   , fDCA(ref.fDCA)
 155   , fNclustersITS(ref.fNclustersITS)
 156   , fMinNprimVtxContrbutor(ref.fMinNprimVtxContrbutor)
 157   , fNEvents(ref.fNEvents)
 158   , fResidualList(ref.fResidualList)
 159   , fPullList(ref.fPullList)
 160   , fDcaList(ref.fDcaList)
 161   , fKfDcaList(ref.fKfDcaList)
 162   , fMcVertexList(ref.fMcVertexList)
 163   , fDataDcaList(ref.fDataDcaList)
 164   , fDataVertexList(ref.fDataVertexList)
 165   , fDataPullList(ref.fDataPullList)
 166   , fMcPidList(ref.fMcPidList)
 167   , fDataPidList(ref.fDataPidList)
 168   , fHfeDcaList(ref.fHfeDcaList)
 169   , fHfeDataDcaList(ref.fHfeDataDcaList)
 170   , fOutput(ref.fOutput)
 171 {
 172   //
 173   // Copy Constructor
 174   //
 175   AliInfo("Copy Constructor");
 176   ref.Copy(*this);
 177 }
 178
 179 //____________________________________________________________
 180 AliAnalysisTaskDCA &AliAnalysisTaskDCA::operator=(const AliAnalysisTaskDCA &ref){
 181   //
 182   // Assignment operator
 183   //
 184   if(this == &ref) return *this;
 185   AliAnalysisTaskSE::operator=(ref);
 186   fPlugins = ref.fPlugins;
 187   fCuts = ref.fCuts;
 188   fHFEpid = ref.fHFEpid;
 189   fPIDdetectors = ref.fPIDdetectors;
 190   fPIDstrategy = ref.fPIDstrategy;
 191   fCFM = ref.fCFM;
 192   fDCA = ref.fDCA;
 193   fNclustersITS = ref.fNclustersITS;
 194   fMinNprimVtxContrbutor = ref.fMinNprimVtxContrbutor;
 195   fNEvents = ref.fNEvents;
 196   fResidualList = ref.fResidualList;
 197   fPullList = ref.fPullList;
 198   fDcaList = ref.fDcaList;
 199   fKfDcaList = ref.fKfDcaList;
 200   fMcVertexList = ref.fMcVertexList;
 201   fDataDcaList = ref.fDataDcaList;
 202   fDataVertexList = ref.fDataVertexList;
 203   fDataPullList = ref.fDataPullList;
 204   fMcPidList = ref.fMcPidList;
 205   fDataPidList = ref.fDataPidList;
 206   fHfeDcaList = ref.fHfeDcaList;
 207   fHfeDataDcaList = ref.fHfeDataDcaList;
 208   fOutput = ref.fOutput;
 209
 210   return *this;
 211 }
 212
 213 //____________________________________________________________
 214 AliAnalysisTaskDCA::~AliAnalysisTaskDCA(){
 215   //
 216   // Destructor
 217   //
 218
 219   if(fHFEpid) delete fHFEpid;
 220   if(fCFM) delete fCFM;
 221   if(fDCA) delete fDCA;
 222   if(fNEvents) delete fNEvents;
 223   if(fResidualList){
 224     fResidualList->Clear();
 225     delete fResidualList;
 226   }
 227
 228   if(fPullList){
 229     fPullList->Clear();
 230     delete fPullList;
 231   }
 232
 233   if(fDcaList){
 234     fDcaList->Clear();
 235     delete fDcaList;
 236   }
 237   if(fKfDcaList){
 238     fKfDcaList->Clear();
 239     delete fKfDcaList;
 240   }
 241
 242   if(fMcVertexList){
 243     fMcVertexList->Clear();
 244     delete   fMcVertexList;
 245   }
 246
 247   if(fDataDcaList){
 248     fDataDcaList->Clear();
 249     delete fDataDcaList;
 250   }
 251
 252   if(fDataVertexList){
 253     fDataVertexList->Clear();
 254     delete   fDataVertexList;
 255   }
 256   if(fDataPullList){
 257     fDataPullList->Clear();
 258     delete fDataPullList;
 259   }
 260
 261   if(fMcPidList){
 262     fMcPidList -> Clear();
 263     delete fMcPidList;
 264   }
 265   if(fDataPidList){
 266     fDataPidList -> Clear();
 267     delete fDataPidList;
 268   }
 269
 270   if(fHfeDcaList) {
 271     fHfeDcaList->Clear();
 272     delete fHfeDcaList;
 273   }
 274
 275   if(fHfeDataDcaList) {
 276     fHfeDataDcaList->Clear();
 277     delete fHfeDataDcaList;
 278   }
 279
 280   if(fOutput){
 281     fOutput->Clear();
 282     delete fOutput;
 283   }
 284
 285 }
 286
 287 //____________________________________________________________
 288 void AliAnalysisTaskDCA::UserCreateOutputObjects(){
 289   // create output objects
 290   // fNEvents
 291   // residual and pull
 292   //printf("\n=====UserCreateOutputObjects=====\n");
 293
 294   // Automatic determination of the analysis mode
 295   AliVEventHandler *inputHandler = dynamic_cast<AliVEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 296   if(!inputHandler){
 297     AliError("NoEvent Handler available");
 298     return;
 299   }
 300
 301   if(!TString(inputHandler->IsA()->GetName()).CompareTo("AliAODInputHandler")){
 302     SetAODAnalysis();
 303   } else {
 304     SetESDAnalysis();
 305     if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
 306       SetHasMCData();
 307   }
 308
 309
 310   fNEvents = new TH1I("nEvents", "Number of Events in the Analysis", 5, -0.5, 4.5); // Number of Events neccessary for the analysis and not a QA histogram
 311   if(!fOutput) fOutput = new TList;
 312   // Initialize correction Framework and Cuts
 313   fCFM = new AliCFManager;
 314   MakeParticleContainer();
 315   // Temporary fix: Initialize particle cuts with 0x0
 316   for(Int_t istep = 0; istep < fCFM->GetParticleContainer()->GetNStep(); istep++)
 317     fCFM->SetParticleCutsList(istep, 0x0);
 318   if(!fCuts){
 319     AliWarning("Cuts not available. Default cuts will be used");
 320     fCuts = new AliHFEcuts;
 321     fCuts->CreateStandardCuts();
 322   }
 323
 324   fCuts->Initialize(fCFM);
 325
 326   if(!fHFEpid) AliWarning("Hello, fHFEpid is not available");
 327   cout<<"  Hello this is a cout "<<endl<<endl;
 328
 329   if(GetPlugin(kHFEpid)) {
 330     fHFEpid->SetHasMCData(HasMCData());
 331     fHFEpid->AddDetector("TOF", 0);
 332     fHFEpid->AddDetector("TPC", 1);
 333     cout<<endl<<" ---> TPC and TOF added to the PID"<<endl;
 334                 fHFEpid->ConfigureTOF();
 335     fHFEpid->ConfigureTPCdefaultCut();
 336     fHFEpid->InitializePID();
 337   }
 338
 339   // dca study----------------------------------
 340
 341
 342   if(!fDCA) fDCA = new AliHFEdca;
 343   if(!fResidualList) fResidualList = new TList();
 344   if(!fPullList) fPullList = new TList();
 345   if(!fDcaList) fDcaList = new TList();
 346   if(!fKfDcaList) fKfDcaList = new TList();
 347   if(!fMcVertexList) fMcVertexList = new TList();
 348   if(!fDataDcaList) fDataDcaList = new TList();
 349   if(!fDataVertexList) fDataVertexList = new TList();
 350   if(!fDataPullList) fDataPullList = new TList();
 351   if(!fMcPidList) fMcPidList = new TList();
 352   if(!fDataPidList) fDataPidList = new TList();
 353
 354   if(!fHfeDcaList) fHfeDcaList = new TList();
 355   if(!fHfeDataDcaList) fHfeDataDcaList = new TList();
 356
 357   if(HasMCData()) {
 358     if(GetPlugin(kImpactPar) ) {
 359       fDCA->CreateHistogramsResidual(fResidualList);
 360       fDCA->CreateHistogramsPull(fPullList);
 361       fDCA->CreateHistogramsDca(fDcaList);
 362       fOutput->AddAt(fResidualList,0);
 363       fOutput->AddAt(fPullList,1);
 364       fOutput->AddAt(fDcaList,2);
 365     }
 366     if(GetPlugin(kKFdca)){
 367       fDCA->CreateHistogramsKfDca(fKfDcaList);
 368       fOutput->AddAt(fDcaList,3);
 369     }
 370     if(GetPlugin(kPrimVtx)){//<---
 371       fDCA->CreateHistogramsVertex(fMcVertexList);
 372       fOutput->AddAt(fMcVertexList,4);
 373     }
 374     if(GetPlugin(kCombinedPid)){//<---
 375       fDCA->CreateHistogramsPid(fMcPidList);
 376       fOutput->AddAt(fMcPidList, 5);
 377     }
 378     if(GetPlugin(kHFEpid)){//<---
 379       fDCA->CreateHistogramsHfeDca(fHfeDcaList);
 380       fOutput->AddAt(fHfeDcaList, 6);
 381     }
 382   } // mc case
 383
 384   if(!HasMCData())  {
 385
 386     if(GetPlugin(kPrimVtx)){
 387       fDCA->CreateHistogramsDataVertex(fDataVertexList);
 388       fOutput->AddAt(fDataVertexList,0);
 389     }
 390
 391     if(GetPlugin(kCombinedPid)){
 392       fDCA->CreateHistogramsDataDca(fDataDcaList);
 393       fDCA->CreateHistogramsDataPull(fDataPullList);
 394       fDCA->CreateHistogramsDataPid(fDataPidList);
 395       fOutput->AddAt(fDataDcaList,1);
 396       fOutput->AddAt(fDataPullList,2);
 397       fOutput->AddAt(fDataPidList, 3);
 398     }
 399     if(GetPlugin(kHFEpid)){
 400       fDCA->CreateHistogramsHfeDataDca(fHfeDataDcaList);
 401       fOutput->AddAt(fHfeDataDcaList, 4);
 402     }
 403
 404
 405
 406   }  // data case
 407
 408 }
 409
 410 //____________________________________________________________
 411 void AliAnalysisTaskDCA::UserExec(Option_t *){
 412   //
 413   // Run the analysis
 414   //
 415   //printf("\n=====UserExec=====\n");
 416   if(HasMCData()) printf("WITH MC!\n");
 417
 418   AliDebug(3, "Processing ESD events");
 419
 420   if(!fInputEvent){
 421     AliError("Reconstructed Event not available");
 422     return;
 423   }
 424   if(HasMCData()){
 425     AliDebug(4, Form("MC Event: %p", fMCEvent));
 426     if(!fMCEvent){
 427       AliError("No MC Event, but MC Data required");
 428       return;
 429     }
 430   }
 431
 432   if(!fCuts){
 433     AliError("HFE cuts not available");
 434     return;
 435   }
 436
 437
 438   // protection
 439   if(IsESDanalysis() && HasMCData()){
 440     // Protect against missing MC trees
 441     AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 442     if(!mcH){
 443       AliError("No MC Event Handler available");
 444       return;
 445     }
 446     if(!mcH->InitOk()) return;
 447     if(!mcH->TreeK()) return;
 448     if(!mcH->TreeTR()) return;
 449   }
 450
 451
 452   AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
 453   if(!pidResponse){
 454     AliDebug(1, "Using default PID Response");
 455     pidResponse = AliHFEtools::GetDefaultPID(HasMCData(), fInputEvent->IsA() == AliAODEvent::Class());
 456   }
 457   fHFEpid->SetPIDResponse(pidResponse);
 458   ProcessDcaAnalysis();
 459
 460   PostData(1, fNEvents);
 461   PostData(2, fOutput);
 462 }
 463 //____________________________________________________________
 464 void AliAnalysisTaskDCA::ProcessDcaAnalysis(){
 465
 466   //printf("\n=====ProcessDcaAnalysis=====\n");
 467
 468   //
 469   // Loop ESD
 470   //
 471   AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
 472   if(!fESD){
 473     AliError("ESD Event required for ESD Analysis");
 474       return;
 475   }
 476
 477   AliMCEvent *fMC = 0x0;
 478   if(HasMCData()){
 479     fMC = dynamic_cast<AliMCEvent*>(fMCEvent);
 480     if(!fMC){
 481       AliError("MC Event required for Analysis");
 482       return;
 483     }
 484   }
 485
 486   fNEvents->Fill(1);  // original event number before cut
 487   fDCA->ApplyExtraCuts(fESD,fMinNprimVtxContrbutor);  // cut on primVtx contributors
 488   fNEvents->Fill(3);  // events number after cut
 489   fCFM->SetRecEventInfo(fESD);
 490
 491   // event cut level
 492   if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fESD)) return;
 493
 494   AliESDtrack *track = 0x0;
 495   AliMCParticle *mctrack = 0x0;
 496   AliESDVertex *vtxESDSkip = 0x0;
 497   AliHFEpidObject hfetrack;
 498
 499   for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
 500
 501     track = fESD->GetTrack(itrack);
 502     if(HasMCData()) mctrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(TMath::Abs(track->GetLabel())));
 503
 504     // RecPrim: primary cuts
 505     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim + AliHFEcuts::kNcutStepsMCTrack, track)) continue;
 506     // RecKine: ITSTPC cuts
 507     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack, track)) continue;
 508     // HFEcuts: ITS layers cuts
 509     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsITS + AliHFEcuts::kNcutStepsMCTrack, track)) continue;
 510
 511     if(track->GetITSclusters(0)<=fNclustersITS) continue;  // require number of ITS clusters
 512
 513     // track accepted, do PID
 514     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
 515     hfetrack.SetRecTrack(track);
 516     if(HasMCData()) hfetrack.SetMCTrack(mctrack);
 517
 518     //printf("Track %d passed all the cuts!\n",itrack);
 519
 520     if(HasMCData()){
 521       if(GetPlugin(kPrimVtx))
 522         fDCA->FillHistogramsVtx(fESD, fMC);
 523       if(GetPlugin(kImpactPar))
 524         fDCA->FillHistogramsDca(fESD, track, fMC);
 525       if(GetPlugin(kKFdca))
 526         fDCA->FillHistogramsKfDca(fESD, track, fMC);
 527       if(GetPlugin(kCombinedPid))
 528         fDCA->FillHistogramsPid(track, fMC);
 529       if(GetPlugin(kHFEpid)) { // data-like
 530         if(fHFEpid->IsSelected(&hfetrack)){
 531
 532           //      printf("Found an electron in p+p collision! from HFE pid \n");
 533           if(!vtxESDSkip){
 534             // method from Andrea D 28.05.2010
 535             AliVertexerTracks *vertexer = new AliVertexerTracks(fESD->GetMagneticField());
 536             vertexer->SetITSMode();
 537             vertexer->SetMinClusters(fNclustersITS);
 538             Int_t skipped[2];
 539             skipped[0] = (Int_t)track->GetID();
 540             vertexer->SetSkipTracks(1,skipped);
 541             vtxESDSkip = (AliESDVertex*)vertexer->FindPrimaryVertex(fESD);
 542             delete vertexer; vertexer = NULL;
 543             if(vtxESDSkip->GetNContributors()<fMinNprimVtxContrbutor) continue;
 544           }
 545           //printf("\n[ABOUT TO FILL HFE DCA: MC!]\n");
 546           fDCA->FillHistogramsHfeDataDca(fESD, track, vtxESDSkip);
 547         }
 548       } // plugin for hfepid
 549     }  // MC
 550
 551     if(!HasMCData()){
 552       if(GetPlugin(kPrimVtx))
 553         fDCA->FillHistogramsDataVtx(fESD);
 554       if(GetPlugin(kCombinedPid)) {
 555
 556         // method from Andrea D 28.05.2010
 557         AliVertexerTracks *vertexer = new AliVertexerTracks(fESD->GetMagneticField());
 558         vertexer->SetITSMode();
 559         vertexer->SetMinClusters(fNclustersITS);
 560         Int_t skipped[2];
 561         skipped[0] = (Int_t)track->GetID();
 562         vertexer->SetSkipTracks(1,skipped);
 563         vtxESDSkip = (AliESDVertex*)vertexer->FindPrimaryVertex(fESD);
 564         delete vertexer; vertexer = NULL;
 565         if(vtxESDSkip->GetNContributors()<fMinNprimVtxContrbutor) continue;
 566
 567         fDCA->FillHistogramsDataDca(fESD, track, vtxESDSkip);
 568         fDCA->FillHistogramsDataPid(track);
 569       }
 570       if(GetPlugin(kHFEpid)) {
 571         if(fHFEpid->IsSelected(&hfetrack)) {
 572           //      printf("Found an electron in p+p collision! from HFE pid \n");
 573           if(!vtxESDSkip){
 574             // method from Andrea D 28.05.2010
 575             AliVertexerTracks *vertexer = new AliVertexerTracks(fESD->GetMagneticField());
 576             vertexer->SetITSMode();
 577             vertexer->SetMinClusters(fNclustersITS);
 578             Int_t skipped[2];
 579             skipped[0] = (Int_t)track->GetID();
 580             vertexer->SetSkipTracks(1,skipped);
 581             vtxESDSkip = (AliESDVertex*)vertexer->FindPrimaryVertex(fESD);
 582             delete vertexer; vertexer = NULL;
 583             if(vtxESDSkip->GetNContributors()<fMinNprimVtxContrbutor) continue;
 584           }
 585         printf("\n[ABOUT TO FILL HFE DCA: DATA!]\n");
 586           fDCA->FillHistogramsHfeDataDca(fESD, track,vtxESDSkip);
 587         }
 588       } // plugin for hfepid
 589     }  // data case
 590
 591   } // track loop
 592
 593 }
 594
 595
 596 //____________________________________________________________
 597 void AliAnalysisTaskDCA::Terminate(Option_t *){
 598   //
 599   // Terminate not implemented at the moment
 600   //
 601   //printf("\n=====Terminate=====\n");
 602
 603   if(GetPlugin(kPostProcess)){
 604     fOutput = dynamic_cast<TList *>(GetOutputData(1));
 605     if(!fOutput){
 606       AliError("Results not available");
 607       return;
 608     }
 609     PostProcess();
 610   }
 611
 612 }
 613
 614
 615 //____________________________________________________________
 616 void AliAnalysisTaskDCA::Load(TString filename){
 617
 618   //printf("\n=====Load=====\n");
 619
 620   // no need for postprocessing for the moment
 621   TFile *input = TFile::Open(filename.Data());
 622   if(!input || input->IsZombie()){
 623     AliError("Cannot read file");
 624     return;
 625   }
 626
 627   input->Close();
 628   delete input;
 629
 630
 631 }
 632
 633 //____________________________________________________________
 634 void AliAnalysisTaskDCA::PostProcess(){
 635   // do post processing
 636   // should do fitting here for dca resolution
 637   // moved to an external macro to do the job
 638
 639   //printf("\n=====PostProcess=====\n");
 640   Load("HFEdca.root");
 641   TCanvas *c1 = new TCanvas("c1", "number of analyzed events", 300, 400);
 642   fNEvents->Draw();
 643   c1->SaveAs("temp.png");
 644
 645 }
 646
 647
 648
 649
 650 //____________________________________________________________
 651 void AliAnalysisTaskDCA::PrintStatus() const {
 652
 653   //
 654   // Print Analysis status
 655   //
 656   printf("\n\tAnalysis Settings\n\t========================================\n");
 657   printf("\t Running on %s\n", !HasMCData()?"p+p collision data":"MC sample");
 658   printf("\t Cuts: %s\n", (fCuts != NULL) ? "YES" : "NO");
 659   printf("\t Impact parameter analysis is %s\n", GetPlugin(kImpactPar)?"ON":"OFF");
 660   printf("\t Using AliKFParticle for analysis? %s\n", GetPlugin(kKFdca)?"ON":"OFF");
 661   printf("\t Primary vertex analysis is %s\n", GetPlugin(kPrimVtx)?"ON":"OFF");
 662   printf("\t Combined pid analysis is %s\n", GetPlugin(kCombinedPid)?"ON":"OFF");
 663   printf("\t HFE pid analysis is %s\n", GetPlugin(kHFEpid)?"ON":"OFF");
 664   printf("\t Post process analysis is %s\n", GetPlugin(kPostProcess)?"ON":"OFF");
 665   printf("\t ");
 666   printf("\n");
 667 }
 668
 669 //__________________________________________
 670 void AliAnalysisTaskDCA::SwitchOnPlugin(Int_t plug){
 671   //
 672   // Switch on Plugin
 673   // Available:
 674   //  - analyze impact parameter
 675   //  - Post Processing
 676
 677   AliDebug(2,Form("SwitchOnPlugin %d",plug));
 678
 679   switch(plug){
 680   case kPostProcess:
 681     SETBIT(fPlugins, plug);
 682     break;
 683   case kImpactPar:
 684     SETBIT(fPlugins, plug);
 685     break;
 686   case kPrimVtx:
 687     SETBIT(fPlugins, plug);
 688     break;
 689   case kCombinedPid:
 690     SETBIT(fPlugins, plug);
 691     break;
 692   case kHFEpid:
 693     SETBIT(fPlugins, plug);
 694     break;
 695   case kKFdca:
 696     SETBIT(fPlugins, plug);
 697     break;
 698   default:
 699     AliError("Unknown Plugin");
 700   };
 701 }
 702
 703
 704 //____________________________________________________________
 705 void AliAnalysisTaskDCA::MakeParticleContainer(){
 706
 707   //printf("\n=====MakeParticleContainer=====\n");
 708   //
 709   // Create the particle container (borrowed from AliAnalysisTaskHFE)
 710   //
 711   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
 712   const Double_t kPtmin = 0.1, kPtmax = 10.;
 713   const Double_t kEtamin = -0.9, kEtamax = 0.9;
 714   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
 715
 716   //arrays for the number of bins in each dimension
 717   Int_t iBin[kNvar];
 718   iBin[0] = 40; //bins in pt
 719   iBin[1] =  8; //bins in eta
 720   iBin[2] = 18; // bins in phi
 721
 722   //arrays for lower bounds :
 723   Double_t* binEdges[kNvar];
 724   for(Int_t ivar = 0; ivar < kNvar; ivar++)
 725     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
 726
 727   //values for bin lower bounds
 728   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);
 729   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
 730   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
 731
 732   //one "container" for MC
 733   const Int_t kNcutStepsESDtrack = AliHFEcuts::kNcutStepsRecTrack + 1;
 734   const Int_t kNcutStepsTrack = AliHFEcuts::kNcutStepsMCTrack + kNcutStepsESDtrack;
 735   AliCFContainer* container = new AliCFContainer("container","container for tracks", (kNcutStepsTrack + 1 + 2*(kNcutStepsESDtrack + 1)), kNvar, iBin);
 736
 737   //setting the bin limits
 738   for(Int_t ivar = 0; ivar < kNvar; ivar++)
 739     container -> SetBinLimits(ivar, binEdges[ivar]);
 740   fCFM->SetParticleContainer(container);
 741
 742   //create correlation matrix for unfolding
 743   Int_t thnDim[2*kNvar];
 744   for (int k=0; k<kNvar; k++) {
 745     //first half  : reconstructed
 746     //second half : MC
 747     thnDim[k]      = iBin[k];
 748     thnDim[k+kNvar] = iBin[k];
 749   }
 750
 751
 752 }
 753
 754 //____________________________________________________________
 755 void AliAnalysisTaskDCA::AddPIDdetector(TString detector){
 756
 757   //
 758   // Adding PID detector to the task
 759   //
 760   //printf("\n=====AddPIDdetector=====\n");
 761
 762   if(!fPIDdetectors.Length())
 763     fPIDdetectors = detector;
 764   else
 765     fPIDdetectors += ":" + detector;
 766 }
 767