CORRFW/test/AliCFTaskForUnfolding.C

   1 //DEFINITION OF A FEW CONSTANTS
   2 const Double_t ptmin  =   0.0 ;
   3 const Double_t ptmax  =   1.0 ;
   4 const Double_t etamin =  -1.5 ;
   5 const Double_t etamax =   1.5 ;
   6 const Int_t    PDG    =    211;
   7 const Int_t    minclustersTPC = 40 ;
   8 //----------------------------------------------------
   9
  10 /*
  11
  12 Macro to prepare the necessary objects to perform spectrum unfolding
  13 (to be used by the class AliCFUnfolding)
  14
  15 Note that the efficiency calculation (and therfore the container filling)
  16 has to be done using MC values (and not reconstructed values)
  17
  18 This macro creates the following objects :
  19 - (AliCFContainer) container    : used to calculate the efficiency (see AliCFEffGrid)
  20 - (THnSparseD)     correlation  : this is the response matrix
  21
  22 Once you have run this macro, you may launch unfolding using
  23 the example macro CORRFW/test/testUnfolding.C
  24
  25 */
  26
  27 Bool_t AliCFTaskForUnfolding()
  28 {
  29
  30   TBenchmark benchmark;
  31   benchmark.Start("AliSingleTrackTask");
  32
  33   AliLog::SetGlobalDebugLevel(0);
  34
  35   Load() ; //load the required libraries
  36
  37   TChain * analysisChain ;
  38   printf("\n\nRunning on local file, please check the path\n\n");
  39
  40   analysisChain = new TChain("esdTree");
  41   analysisChain->Add("your_data_path/001/AliESDs.root");
  42   analysisChain->Add("your_data_path/002/AliESDs.root");
  43
  44
  45   Info("AliCFTaskForUnfolding",Form("CHAIN HAS %d ENTRIES",(Int_t)analysisChain->GetEntries()));
  46
  47   //CONTAINER DEFINITION
  48   Info("AliCFTaskForUnfolding","SETUP CONTAINER");
  49   //the sensitive variables (2 in this example), their indices
  50   UInt_t ipt = 0;
  51   UInt_t ieta  = 1;
  52   //Setting up the container grid...
  53   UInt_t nstep = 3 ; //number of selection steps MC
  54   const Int_t nvar   = 2 ; //number of variables on the grid:pt,eta
  55   const Int_t nbin[nvar] = {20,20} ;
  56
  57   //arrays for the number of bins in each dimension
  58   Int_t iBin[nvar];
  59   iBin[0]=nbin[0];
  60   iBin[1]=nbin[1];
  61
  62   //arrays for lower bounds :
  63   Double_t *binLim0=new Double_t[nbin[0]+1];
  64   Double_t *binLim1=new Double_t[nbin[1]+1];
  65
  66   //values for bin lower bounds
  67   for(Int_t i=0; i<=nbin[0]; i++) binLim0[i]=(Double_t)ptmin  + ( ptmax- ptmin)/nbin[0]*(Double_t)i ;
  68   for(Int_t i=0; i<=nbin[1]; i++) binLim1[i]=(Double_t)etamin + (etamax-etamin)/nbin[1]*(Double_t)i ;
  69
  70   //one "container" for MC
  71   AliCFContainer* container = new AliCFContainer("container","container for tracks",nstep,nvar,iBin);
  72   //setting the bin limits
  73   container -> SetBinLimits(ipt,binLim0);
  74   container -> SetBinLimits(ieta ,binLim1);
  75   container -> SetVarTitle(ipt,"pT");
  76   container -> SetVarTitle(ieta,"#eta");
  77
  78   // Gen-Level kinematic cuts
  79   AliCFTrackKineCuts *mcKineCuts = new AliCFTrackKineCuts("mcKineCuts","MC-level kinematic cuts");
  80   mcKineCuts->SetPtRange (ptmin ,ptmax );
  81   mcKineCuts->SetEtaRange(etamin,etamax);
  82
  83   //Particle-Level cuts:
  84   AliCFParticleGenCuts* mcGenCuts = new AliCFParticleGenCuts("mcGenCuts","MC particle generation cuts");
  85   mcGenCuts->SetRequireIsPrimary();
  86   mcGenCuts->SetRequirePdgCode(PDG);
  87
  88   // Rec-Level kinematic cuts
  89   AliCFTrackKineCuts *recKineCuts = new AliCFTrackKineCuts("recKineCuts","rec-level kine cuts");
  90   recKineCuts->SetPtRange( ptmin, ptmax);
  91   recKineCuts->SetPtRange(etamin,etamax);
  92
  93   AliCFTrackQualityCuts *recQualityCuts = new AliCFTrackQualityCuts("recQualityCuts","rec-level quality cuts");
  94   recQualityCuts->SetMinNClusterTPC(minclustersTPC);
  95
  96   AliCFTrackIsPrimaryCuts *recIsPrimaryCuts = new AliCFTrackIsPrimaryCuts("recIsPrimaryCuts","rec-level isPrimary cuts");
  97   recIsPrimaryCuts->SetMaxNSigmaToVertex(3);
  98
  99   printf("CREATE MC KINE CUTS\n");
 100   TObjArray* mcList = new TObjArray(0) ;
 101   mcList->AddLast(mcKineCuts);
 102   mcList->AddLast(mcGenCuts);
 103
 104   printf("CREATE RECONSTRUCTION CUTS\n");
 105   TObjArray* recList = new TObjArray(0) ;
 106   recList->AddLast(recKineCuts);
 107   recList->AddLast(recQualityCuts);
 108   recList->AddLast(recIsPrimaryCuts);
 109
 110   TObjArray* emptyList = new TObjArray(0);
 111
 112   //CREATE THE INTERFACE TO CORRECTION FRAMEWORK USED IN THE TASK
 113   printf("CREATE INTERFACE AND CUTS\n");
 114   AliCFManager* man = new AliCFManager() ;
 115   man->SetNStepEvent(0);
 116   man->SetParticleContainer(container);
 117   man->SetParticleCutsList(0,mcList);
 118   man->SetParticleCutsList(1,recList);
 119   man->SetParticleCutsList(2,emptyList); // this is special step for monte carlo spectrum
 120
 121   //CREATE THE TASK
 122   printf("CREATE TASK\n");
 123   // create the task
 124   AliCFTaskForUnfolding *task = new AliCFTaskForUnfolding("AliCFTaskForUnfolding");
 125   task->SetCFManager(man); //here is set the CF manager
 126
 127   //create correlation matrix for unfolding
 128   Int_t thnDim[2*nvar];
 129   for (int k=0; k<nvar; k++) {
 130     //first half  : reconstructed
 131     //second half : MC
 132     thnDim[k]      = nbin[k];
 133     thnDim[k+nvar] = nbin[k];
 134   }
 135   THnSparseD* correlation = new THnSparseD("correlation","THnSparse with correlations",2*nvar,thnDim);
 136   Double_t** binEdges = new Double_t[nvar];
 137   for (int k=0; k<nvar; k++) {
 138     binEdges[k]=new Double_t[nbin[k]+1];
 139     container->GetBinLimits(k,binEdges[k]);
 140     correlation->SetBinEdges(k,binEdges[k]);
 141     correlation->SetBinEdges(k+nvar,binEdges[k]);
 142   }
 143   correlation->Sumw2();
 144   task->SetCorrelationMatrix(correlation);
 145   //---
 146
 147   //SETUP THE ANALYSIS MANAGER TO READ INPUT CHAIN AND WRITE DESIRED OUTPUTS
 148   printf("CREATE ANALYSIS MANAGER\n");
 149   // Make the analysis manager
 150   AliAnalysisManager *mgr = new AliAnalysisManager("TestManager");
 151   mgr->SetAnalysisType(AliAnalysisManager::kLocalAnalysis);
 152
 153   AliMCEventHandler*  mcHandler = new AliMCEventHandler();
 154   mgr->SetMCtruthEventHandler(mcHandler);
 155
 156   AliInputEventHandler* dataHandler = new AliESDInputHandler();
 157   mgr->SetInputEventHandler(dataHandler);
 158
 159   // Create and connect containers for input/output
 160
 161   //------ input data ------
 162   AliAnalysisDataContainer *cinput0  = mgr->CreateContainer("cchain0",TChain::Class(),AliAnalysisManager::kInputContainer);
 163
 164   // ----- output data -----
 165
 166   //slot 0 : default output tree (by default handled by AliAnalysisTaskSE)
 167   AliAnalysisDataContainer *coutput0 = mgr->CreateContainer("ctree0", TTree::Class(),AliAnalysisManager::kOutputContainer,"output.root");
 168
 169   //now comes user's output objects :
 170
 171   // output TH1I for event counting
 172   AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("chist0", TH1I::Class(),AliAnalysisManager::kOutputContainer,"output.root");
 173   // output Correction Framework Container (for acceptance & efficiency calculations)
 174   AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("ccontainer0", AliCFContainer::Class(),AliAnalysisManager::kOutputContainer,"output.root");
 175   AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("corr0", THnSparseD::Class(),AliAnalysisManager::kOutputContainer,"output.root");
 176
 177   cinput0->SetData(analysisChain);
 178
 179   mgr->AddTask(task);
 180   mgr->ConnectInput(task,0,mgr->GetCommonInputContainer());
 181   mgr->ConnectOutput(task,0,coutput0);
 182   mgr->ConnectOutput(task,1,coutput1);
 183   mgr->ConnectOutput(task,2,coutput2);
 184   mgr->ConnectOutput(task,3,coutput3);
 185
 186   printf("READY TO RUN\n");
 187   //RUN !!!
 188   if (mgr->InitAnalysis()) {
 189     mgr->PrintStatus();
 190     mgr->StartAnalysis("local",analysisChain);
 191   }
 192
 193   benchmark.Stop("AliSingleTrackTask");
 194   benchmark.Show("AliSingleTrackTask");
 195
 196   return kTRUE ;
 197 }
 198
 199 void Load() {
 200
 201   //load the required aliroot libraries
 202   gSystem->Load("libANALYSIS") ;
 203   gSystem->Load("libANALYSISalice") ;
 204   gSystem->Load("libCORRFW.so") ;
 205
 206   //compile online the task class
 207   gSystem->SetIncludePath("-I. -I$ALICE_ROOT/include -I$ROOTSYS/include");
 208   gROOT->LoadMacro("./AliCFTaskForUnfolding.cxx+");
 209 }