PWG2/runProtonCorrection.C

   1 Bool_t runProtonCorrection(Int_t stats = 0, const char* dataset = 0x0) {
   2   //macro used to extract the correction maps
   3   //using the official correction framework of ALICE
   4   //for protons and antiprotons
   5   //Author: Panos.Christakoglou@cern.ch
   6
   7   //________________________________________//
   8   //Connect to proof
   9   //TProof::Reset("proof://lxb6046.cern.ch");
  10   TProof::Open("proof://lxb6046.cern.ch");
  11
  12   // Enable the STEERBase Package
  13   gProof->UploadPackage("STEERBase.par");
  14   gProof->EnablePackage("STEERBase");
  15   // Enable the ESD Package
  16   gProof->UploadPackage("ESD.par");
  17   gProof->EnablePackage("ESD");
  18   // Enable the AOD Package
  19   gProof->UploadPackage("AOD.par");
  20   gProof->EnablePackage("AOD");
  21   // Enable the Analysis Package
  22   gProof->UploadPackage("ANALYSIS.par");
  23   gProof->EnablePackage("ANALYSIS");
  24   gProof->UploadPackage("ANALYSISalice.par");
  25   gProof->EnablePackage("ANALYSISalice");
  26   // Enable the CORRFW Package
  27   gProof->UploadPackage("CORRFW.par");
  28   gProof->EnablePackage("CORRFW");
  29
  30   gProof->Load("./AliProtonCorrectionTask.cxx+g");
  31
  32   //________________________________________//
  33   //Container definition
  34   //Variables of the GRID
  35   //For the time being: y-pT
  36   //Next step: add Vz
  37   //Setting up the container grid...
  38   //===============//
  39   //Global tracking//
  40   //===============//
  41   const Double_t ymin  = -1.0;
  42   const Double_t ymax  =  1.0;
  43   const Int_t nbin1 = 20; //bins in y
  44   const Double_t ptmin =  0.4;
  45   const Double_t ptmax =  3.1;
  46   const Int_t nbin2 = 27; //bins in pT
  47   //===============//
  48   //  TPC tracking //
  49   //===============//
  50   /*const Double_t ymin  = -0.5;
  51   const Double_t ymax  =  0.5;
  52   const Int_t nbin1 = 10; //bins in y
  53   const Double_t ptmin =  0.4;
  54   const Double_t ptmax =  0.9;
  55   const Int_t nbin2 = 15;*/ //bins in pT
  56
  57   //Setting up the container grid...
  58   UInt_t nstep = 4; //number of selection steps MC
  59   UInt_t iy  = 0;
  60   UInt_t ipT = 1;
  61   const Int_t nvar = 2; //number of variables on the grid:y-pT
  62     //arrays for the number of bins in each dimension
  63   Int_t iBin[nvar];
  64   iBin[0]=nbin1;
  65   iBin[1]=nbin2;
  66   //arrays for lower bounds :
  67   Double_t *binLim1=new Double_t[nbin1+1];
  68   Double_t *binLim2=new Double_t[nbin2+1];
  69   //values for bin lower bounds
  70   for(Int_t i=0; i<=nbin1; i++) binLim1[i]=(Double_t)ymin  + (ymax-ymin)  /nbin1*(Double_t)i;
  71   for(Int_t i=0; i<=nbin2; i++) binLim2[i]=(Double_t)ptmin + (ptmax-ptmin)/nbin2*(Double_t)i;
  72   //CF container for protons
  73   AliCFContainer* containerProtons = new AliCFContainer("containerProtons","container for protons",nstep,nvar,iBin);
  74   //setting the bin limits
  75   containerProtons->SetBinLimits(iy,binLim1);
  76   containerProtons->SetBinLimits(ipT,binLim2);
  77  //CF container for antiprotons
  78   AliCFContainer* containerAntiProtons = new AliCFContainer("containerAntiProtons","container for antiprotons",nstep,nvar,iBin);
  79   //setting the bin limits
  80   containerAntiProtons->SetBinLimits(iy,binLim1);
  81   containerAntiProtons->SetBinLimits(ipT,binLim2);
  82
  83   //________________________________________//
  84   // SET TLIST FOR QA HISTOS
  85   TList* qaList = new TList();
  86   //Cuts
  87   const Int_t    mintrackrefsTPC = 2;
  88   const Int_t    mintrackrefsITS = 3;
  89   const Int_t    chargeProtons  = 1;
  90   const Int_t    PDGProtons = 2212;
  91   const Int_t    chargeAntiProtons  = -1;
  92   const Int_t    PDGAntiProtons = -2212;
  93
  94   const Int_t    minclustersTPC = 50;
  95   const Float_t  maxChi2PerTPCCluster = 3.5;
  96   const Float_t  maxCov11 = 2.0;
  97   const Float_t  maxCov22 = 2.0;
  98   const Float_t  maxCov33 = 0.5;
  99   const Float_t  maxCov44 = 0.5;
 100   const Float_t  maxCov55 = 2.0;
 101   const Float_t  maxSigmaToVertexTPC = 2.5;
 102
 103   const Int_t    minclustersITS = 5;
 104   const Float_t  maxSigmaToVertex = 2.5;
 105
 106   // Gen-Level kinematic cuts
 107   AliCFTrackKineCuts *mcKineCutsProtons = new AliCFTrackKineCuts("mcKineCutsProtons",
 108                                                                  "MC-level kinematic cuts");
 109   mcKineCutsProtons->SetPtRange(ptmin,ptmax);
 110   mcKineCutsProtons->SetRapidityRange(ymin,ymax);
 111   mcKineCutsProtons->SetChargeMC(chargeProtons);
 112   mcKineCutsProtons->SetQAOn(qaList);
 113
 114   AliCFTrackKineCuts *mcKineCutsAntiProtons = new AliCFTrackKineCuts("mcKineCutsAntiProtons",
 115                                                                      "MC-level kinematic cuts");
 116   mcKineCutsAntiProtons->SetPtRange(ptmin,ptmax);
 117   mcKineCutsAntiProtons->SetRapidityRange(ymin,ymax);
 118   mcKineCutsAntiProtons->SetChargeMC(chargeAntiProtons);
 119   mcKineCutsAntiProtons->SetQAOn(qaList);
 120
 121   //Particle-Level cuts:
 122   AliCFParticleGenCuts* mcGenCuts = new AliCFParticleGenCuts("mcGenCuts",
 123                                                              "MC particle generation cuts");
 124   mcGenCuts->SetRequireIsPrimary();
 125   mcGenCuts->SetRequirePdgCode(PDGProtons);
 126   mcGenCuts->SetQAOn(qaList);
 127
 128   //Acceptance Cuts
 129   AliCFAcceptanceCuts *mcAccCuts = new AliCFAcceptanceCuts("mcAccCuts",
 130                                                            "MC acceptance cuts");
 131   mcAccCuts->SetMinNHitITS(mintrackrefsITS);
 132   mcAccCuts->SetMinNHitTPC(mintrackrefsTPC);
 133   mcAccCuts->SetQAOn(qaList);
 134
 135   // Rec-Level kinematic cuts
 136   AliCFTrackKineCuts *recKineCutsProtons = new AliCFTrackKineCuts("recKineCutsProtons",
 137                                                                   "rec-level kine cuts");
 138   recKineCutsProtons->SetPtRange(ptmin,ptmax);
 139   recKineCutsProtons->SetRapidityRange(ymin,ymax);
 140   recKineCutsProtons->SetChargeRec(chargeProtons);
 141   recKineCutsProtons->SetQAOn(qaList);
 142
 143   AliCFTrackKineCuts *recKineCutsAntiProtons = new AliCFTrackKineCuts("recKineCutsAntiProtons",
 144                                                                       "rec-level kine cuts");
 145   recKineCutsAntiProtons->SetPtRange(ptmin,ptmax);
 146   recKineCutsAntiProtons->SetRapidityRange(ymin,ymax);
 147   recKineCutsAntiProtons->SetChargeRec(chargeAntiProtons);
 148   recKineCutsAntiProtons->SetQAOn(qaList);
 149
 150   AliCFTrackQualityCuts *recQualityCuts = new AliCFTrackQualityCuts("recQualityCuts",
 151                                                                     "rec-level quality cuts");
 152   recQualityCuts->SetMinNClusterTPC(minclustersTPC);
 153   recQualityCuts->SetMaxChi2PerClusterTPC(maxChi2PerTPCCluster);
 154   recQualityCuts->SetMaxCovDiagonalElements(maxCov11,maxCov22,maxCov33,maxCov44,maxCov55);
 155   recQualityCuts->SetRequireTPCRefit(kTRUE);
 156
 157   //recQualityCuts->SetRequireITSRefit(kTRUE);
 158   recQualityCuts->SetQAOn(qaList);
 159
 160   AliCFTrackIsPrimaryCuts *recIsPrimaryCuts = new AliCFTrackIsPrimaryCuts("recIsPrimaryCuts",
 161                                                                           "rec-level isPrimary cuts");
 162   recIsPrimaryCuts->SetMaxNSigmaToVertex(maxSigmaToVertex);
 163   recIsPrimaryCuts->SetQAOn(qaList);
 164
 165   AliCFTrackCutPid* cutPID = new AliCFTrackCutPid("cutPID","ESD_PID");
 166   int n_species = AliPID::kSPECIES;
 167   Double_t* prior = new Double_t[n_species];
 168
 169   prior[0] = 0.0244519;
 170   prior[1] = 0.0143988;
 171   prior[2] = 0.805747 ;
 172   prior[3] = 0.0928785;
 173   prior[4] = 0.0625243;
 174
 175   cutPID->SetPriors(prior);
 176   cutPID->SetProbabilityCut(0.0);
 177   cutPID->SetDetectors("ITS TPC TOF");
 178   cutPID->SetParticleType(AliPID::kProton  , kTRUE);
 179   cutPID->SetQAOn(qaList);
 180
 181   //________________________________________//
 182   TObjArray* mcListProtons = new TObjArray(0);
 183   mcListProtons->AddLast(mcKineCutsProtons);
 184   mcListProtons->AddLast(mcGenCuts);
 185   TObjArray* mcListAntiProtons = new TObjArray(0);
 186   mcListAntiProtons->AddLast(mcKineCutsAntiProtons);
 187   mcListAntiProtons->AddLast(mcGenCuts);
 188
 189   printf("CREATE ACCEPTANCE CUTS\n");
 190   TObjArray* accList = new TObjArray(0);
 191   accList->AddLast(mcAccCuts);
 192
 193   printf("CREATE RECONSTRUCTION CUTS\n");
 194   TObjArray* recListProtons = new TObjArray(0);
 195   recListProtons->AddLast(recKineCutsProtons);
 196   recListProtons->AddLast(recQualityCuts);
 197   recListProtons->AddLast(recIsPrimaryCuts);
 198   TObjArray* recListAntiProtons = new TObjArray(0);
 199   recListAntiProtons->AddLast(recKineCutsAntiProtons);
 200   recListAntiProtons->AddLast(recQualityCuts);
 201   recListAntiProtons->AddLast(recIsPrimaryCuts);
 202
 203   printf("CREATE PID CUTS\n");
 204   TObjArray* fPIDCutList = new TObjArray(0);
 205   fPIDCutList->AddLast(cutPID);
 206
 207   //________________________________________//
 208   //CREATE THE INTERFACE TO CORRECTION FRAMEWORK USED IN THE TASK
 209   AliCFManager* manProtons = new AliCFManager();
 210   manProtons->SetParticleContainer(containerProtons);
 211   manProtons->SetParticleCutsList(AliCFManager::kPartGenCuts,mcListProtons);
 212   manProtons->SetParticleCutsList(AliCFManager::kPartAccCuts,accList);
 213   manProtons->SetParticleCutsList(AliCFManager::kPartRecCuts,recListProtons);
 214   manProtons->SetParticleCutsList(AliCFManager::kPartSelCuts,fPIDCutList);
 215
 216   AliCFManager* manAntiProtons = new AliCFManager();
 217   manAntiProtons->SetParticleContainer(containerAntiProtons);
 218   manAntiProtons->SetParticleCutsList(AliCFManager::kPartGenCuts,mcListAntiProtons);
 219   manAntiProtons->SetParticleCutsList(AliCFManager::kPartAccCuts,accList);
 220   manAntiProtons->SetParticleCutsList(AliCFManager::kPartRecCuts,recListAntiProtons);
 221   manAntiProtons->SetParticleCutsList(AliCFManager::kPartSelCuts,fPIDCutList);
 222
 223   //________________________________________//
 224   //CREATE THE TASK
 225   AliProtonCorrectionTask *task = new AliProtonCorrectionTask("AliProtonCorrectionTask");
 226   task->SetCFManagerProtons(manProtons); //here is set the CF manager
 227   task->SetCFManagerAntiProtons(manAntiProtons); //here is set the CF manager
 228   task->SetQAList(qaList);
 229
 230   //SETUP THE ANALYSIS MANAGER TO READ INPUT CHAIN AND WRITE DESIRED OUTPUTS
 231   printf("CREATE ANALYSIS MANAGER\n");
 232   // Make the analysis manager
 233   AliAnalysisManager *mgr = new AliAnalysisManager("TestManager");
 234
 235   AliMCEventHandler*  mcHandler = new AliMCEventHandler();
 236   AliESDInputHandler* esdHandler = new AliESDInputHandler();
 237   mgr->SetMCtruthEventHandler(mcHandler);
 238   mgr->SetInputEventHandler(esdHandler);
 239
 240   //------ input data ------
 241   AliAnalysisDataContainer *cinput0  = mgr->CreateContainer("cchain0",TChain::Class(),AliAnalysisManager::kInputContainer);
 242
 243   // ----- output data -----
 244   //slot 0 : default output tree (by default handled by AliAnalysisTaskSE)
 245   AliAnalysisDataContainer *coutput0 = mgr->CreateContainer("ctree0", TTree::Class(),
 246                                                             AliAnalysisManager::kOutputContainer,"corrections.root");
 247   // output TH1I for event counting
 248   AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("chist0", TH1I::Class(),
 249                                                             AliAnalysisManager::kOutputContainer,"corrections.root");
 250   // output Correction Framework Container (for acceptance & efficiency calculations)
 251   AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("ccontainer0",
 252                                                             AliCFContainer::Class(),
 253                                                             AliAnalysisManager::kOutputContainer,"corrections.root");
 254   // output Correction Framework Container (for acceptance & efficiency calculations)
 255   AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("ccontainer1",
 256                                                             AliCFContainer::Class(),
 257                                                             AliAnalysisManager::kOutputContainer,"corrections.root");
 258   // output QA histograms
 259   AliAnalysisDataContainer *coutput4 = mgr->CreateContainer("clist0",
 260                                                             TList::Class(),
 261                                                             AliAnalysisManager::kOutputContainer,"corrections.root");
 262
 263   mgr->AddTask(task);
 264   mgr->ConnectInput(task,0,cinput0);
 265   mgr->ConnectOutput(task,0,coutput0);
 266   mgr->ConnectOutput(task,1,coutput1);
 267   mgr->ConnectOutput(task,2,coutput2);
 268   mgr->ConnectOutput(task,3,coutput3);
 269   mgr->ConnectOutput(task,4,coutput4);
 270
 271   //________________________________________//
 272   if (mgr->InitAnalysis()) {
 273     if(dataset)
 274       mgr->StartAnalysis("proof",dataset,stats);
 275     else {
 276       // You should get this macro and the txt file from:
 277       // http://aliceinfo.cern.ch/Offline/Analysis/CAF/
 278       gROOT->LoadMacro("CreateESDChain.C");
 279       TChain* chain = 0x0;
 280       chain = CreateESDChain("ESD82XX_30K.txt",stats);
 281       chain->SetBranchStatus("*Calo*",0);
 282
 283       mgr->StartAnalysis("proof",chain);
 284     }
 285   }
 286
 287   return kTRUE;
 288 }
 289