PWG3/vertexingHF/AddTaskCFMultiVarMultiStep.C

   1 //DEFINITION OF A FEW CONSTANTS
   2 const Double_t ymin  = -2.1 ;
   3 const Double_t ymax  =  2.1 ;
   4 const Double_t ptmin_0_4 =  0.0 ;
   5 const Double_t ptmax_0_4 =  4.0 ;
   6 const Double_t ptmin_4_8 =  4.0 ;
   7 const Double_t ptmax_4_8 =  8.0 ;
   8 const Double_t ptmin_8_10 =  8.0 ;
   9 const Double_t ptmax_8_10 =  10.0 ;
  10 const Double_t cosmin = -1.05;
  11 const Double_t cosmax =  1.05;
  12 const Double_t cTmin = 0;  // micron
  13 const Double_t cTmax = 500;  // micron
  14 const Double_t dcamin = 0;  // micron
  15 const Double_t dcamax = 500;  // micron
  16 const Double_t d0min = -1000;  // micron
  17 const Double_t d0max = 1000;  // micron
  18 const Double_t d0xd0min = -100000;  // micron
  19 const Double_t d0xd0max = 100000;  // micron
  20 const Double_t phimin = 0.0;
  21 //const Double_t phimax = 2Pi;  // defined in the macro!!!!!!!!!!!!!!
  22 const Int_t    mintrackrefsTPC = 2 ;
  23 const Int_t    mintrackrefsITS = 3 ;
  24 const Int_t    charge  = 1 ;
  25 const Int_t    PDG = 421;
  26 const Int_t    minclustersTPC = 50 ;
  27 // cuts
  28 const Double_t ptmin = 0.1;
  29 const Double_t ptmax = 9999.;
  30 const Double_t etamin = -0.9;
  31 const Double_t etamax = 0.9;
  32 const Int_t    minITSClusters = 5;
  33
  34 //----------------------------------------------------
  35
  36 AliCFHeavyFlavourTaskMultiVarMultiStep *AddTaskCFMultiVarMultiStep()
  37 {
  38
  39         //CONTAINER DEFINITION
  40         Info("AliCFHeavyFlavourTaskMultiVarMultiStep","SETUP CONTAINER");
  41         //the sensitive variables (6 in this example), their indices
  42         UInt_t ipt = 0;
  43         UInt_t iy  = 1;
  44         UInt_t icosThetaStar  = 2;
  45         UInt_t ipTpi  = 3;
  46         UInt_t ipTk  = 4;
  47         UInt_t icT  = 5;
  48         UInt_t idca  = 6;
  49         UInt_t id0pi  = 7;
  50         UInt_t id0K  = 8;
  51         UInt_t id0xd0  = 9;
  52         UInt_t ipointing  = 10;
  53         UInt_t iphi  = 11;
  54
  55         const Double_t phimax = 2*TMath::Pi();
  56
  57         //Setting up the container grid...
  58         UInt_t nstep = 8; //number of selection steps: MC, Acceptance, Vertex, Refit, Reco (no cuts), RecoAcceptance, RecoITSClusters (RecoAcceptance included), RecoPPR (RecoAcceptance+RecoITSCluster included)
  59         const Int_t nvar   = 12 ; //number of variables on the grid:pt, y, cosThetaStar, pTpi, pTk, cT, dca, d0pi, d0K, d0xd0, cosPointingAngle, phi
  60         const Int_t nbin0_0_4  = 8 ; //bins in pt from 0 to 4 GeV
  61         const Int_t nbin0_4_8  = 4 ; //bins in pt from 4 to 8 GeV
  62         const Int_t nbin0_8_10  = 1 ; //bins in pt from 8 to 10 GeV
  63         const Int_t nbin1  = 42 ; //bins in y
  64         const Int_t nbin2  = 42 ; //bins in cosThetaStar
  65         const Int_t nbin3_0_4  = 8 ; //bins in ptPi from 0 to 4 GeV
  66         const Int_t nbin3_4_8  = 4 ; //bins in ptPi from 4 to 8 GeV
  67         const Int_t nbin3_8_10  = 1 ; //bins in ptPi from 8 to 10 GeV
  68         const Int_t nbin4_0_4  = 8 ; //bins in ptKa from 0 to 4 GeV
  69         const Int_t nbin4_4_8  = 4 ; //bins in ptKa from 4 to 8 GeV
  70         const Int_t nbin4_8_10  = 1 ; //bins in ptKa from 8 to 10 GeV
  71         const Int_t nbin5  = 24 ; //bins in cT
  72         const Int_t nbin6  = 24 ; //bins in dca
  73         const Int_t nbin7  = 100 ; //bins in d0pi
  74         const Int_t nbin8  = 100 ; //bins in d0K
  75         const Int_t nbin9  = 80 ; //bins in d0xd0
  76         const Int_t nbin10  = 1050 ; //bins in cosPointingAngle
  77         const Int_t nbin11  = 20 ; //bins in Phi
  78
  79         //arrays for the number of bins in each dimension
  80         Int_t iBin[nvar];
  81         iBin[0]=nbin0_0_4+nbin0_4_8+nbin0_8_10;
  82         iBin[1]=nbin1;
  83         iBin[2]=nbin2;
  84         iBin[3]=nbin3_0_4+nbin3_4_8+nbin3_8_10;
  85         iBin[4]=nbin4_0_4+nbin4_4_8+nbin4_8_10;
  86         iBin[5]=nbin5;
  87         iBin[6]=nbin6;
  88         iBin[7]=nbin7;
  89         iBin[8]=nbin8;
  90         iBin[9]=nbin9;
  91         iBin[10]=nbin10;
  92         iBin[11]=nbin11;
  93
  94         //arrays for lower bounds :
  95         Double_t *binLim0=new Double_t[iBin[0]+1];
  96         Double_t *binLim1=new Double_t[iBin[1]+1];
  97         Double_t *binLim2=new Double_t[iBin[2]+1];
  98         Double_t *binLim3=new Double_t[iBin[3]+1];
  99         Double_t *binLim4=new Double_t[iBin[4]+1];
 100         Double_t *binLim5=new Double_t[iBin[5]+1];
 101         Double_t *binLim6=new Double_t[iBin[6]+1];
 102         Double_t *binLim7=new Double_t[iBin[7]+1];
 103         Double_t *binLim8=new Double_t[iBin[8]+1];
 104         Double_t *binLim9=new Double_t[iBin[9]+1];
 105         Double_t *binLim10=new Double_t[iBin[10]+1];
 106         Double_t *binLim11=new Double_t[iBin[11]+1];
 107
 108         // checking limits
 109         if (ptmax_0_4 != ptmin_4_8) {
 110                 Error("AliCFHeavyFlavourTaskMultiVarMultiStep","max lim 1st range != min lim 2nd range, please check!");
 111         }
 112         if (ptmax_4_8 != ptmin_8_10) {
 113                 Error("AliCFHeavyFlavourTaskMultiVarMultiStep","max lim 2nd range != min lim 3rd range, please check!");
 114         }
 115
 116         // values for bin lower bounds
 117         // pt
 118         for(Int_t i=0; i<=nbin0_0_4; i++) binLim0[i]=(Double_t)ptmin_0_4 + (ptmax_0_4-ptmin_0_4)/nbin0_0_4*(Double_t)i ;
 119         if (binLim0[nbin0_0_4] != ptmin_4_8)  {
 120                 Error("AliCFHeavyFlavourTaskMultiVarMultiStep","Calculated bin lim for pt - 1st range - differs from expected!\n");
 121         }
 122         for(Int_t i=0; i<=nbin0_4_8; i++) binLim0[i+nbin0_0_4]=(Double_t)ptmin_4_8 + (ptmax_4_8-ptmin_4_8)/nbin0_4_8*(Double_t)i ;
 123         if (binLim0[nbin0_0_4+nbin0_4_8] != ptmin_8_10)  {
 124                 Error("AliCFHeavyFlavourTaskMultiVarMultiStep","Calculated bin lim for pt - 2nd range - differs from expected!\n");
 125         }
 126         for(Int_t i=0; i<=nbin0_8_10; i++) binLim0[i+nbin0_0_4+nbin0_4_8]=(Double_t)ptmin_8_10 + (ptmax_8_10-ptmin_8_10)/nbin0_8_10*(Double_t)i ;
 127
 128         // y
 129         for(Int_t i=0; i<=nbin1; i++) binLim1[i]=(Double_t)ymin  + (ymax-ymin)  /nbin1*(Double_t)i ;
 130
 131         // cosThetaStar
 132         for(Int_t i=0; i<=nbin2; i++) binLim2[i]=(Double_t)cosmin  + (cosmax-cosmin)  /nbin2*(Double_t)i ;
 133
 134         // ptPi
 135         for(Int_t i=0; i<=nbin3_0_4; i++) binLim3[i]=(Double_t)ptmin_0_4 + (ptmax_0_4-ptmin_0_4)/nbin3_0_4*(Double_t)i ;
 136         if (binLim3[nbin3_0_4] != ptmin_4_8)  {
 137                 Error("AliCFHeavyFlavourTaskMultiVarMultiStep","Calculated bin lim for ptPi - 1st range - differs from expected!");
 138         }
 139         for(Int_t i=0; i<=nbin3_4_8; i++) binLim3[i+nbin3_0_4]=(Double_t)ptmin_4_8 + (ptmax_4_8-ptmin_4_8)/nbin3_4_8*(Double_t)i ;
 140         if (binLim3[nbin3_0_4+nbin3_4_8] != ptmin_8_10)  {
 141                 Error("AliCFHeavyFlavourTaskMultiVarMultiStep","Calculated bin lim for ptPi - 2nd range - differs from expected!\n");
 142         }
 143         for(Int_t i=0; i<=nbin3_8_10; i++) binLim3[i+nbin3_0_4+nbin3_4_8]=(Double_t)ptmin_8_10 + (ptmax_8_10-ptmin_8_10)/nbin3_8_10*(Double_t)i ;
 144
 145         // ptKa
 146         for(Int_t i=0; i<=nbin4_0_4; i++) binLim4[i]=(Double_t)ptmin_0_4 + (ptmax_0_4-ptmin_0_4)/nbin4_0_4*(Double_t)i ;
 147         if (binLim4[nbin4_0_4] != ptmin_4_8)  {
 148                 Error("AliCFHeavyFlavourTaskMultiVarMultiStep","Calculated bin lim for ptKa - 1st range - differs from expected!");
 149         }
 150         for(Int_t i=0; i<=nbin4_4_8; i++) binLim4[i+nbin4_0_4]=(Double_t)ptmin_4_8 + (ptmax_4_8-ptmin_4_8)/nbin4_4_8*(Double_t)i ;
 151         if (binLim4[nbin4_0_4+nbin4_4_8] != ptmin_8_10)  {
 152                 Error("AliCFHeavyFlavourTaskMultiVarMultiStep","Calculated bin lim for ptKa - 2nd range - differs from expected!\n");
 153         }
 154         for(Int_t i=0; i<=nbin4_8_10; i++) binLim4[i+nbin4_0_4+nbin4_4_8]=(Double_t)ptmin_8_10 + (ptmax_8_10-ptmin_8_10)/nbin4_8_10*(Double_t)i ;
 155
 156         // cT
 157         for(Int_t i=0; i<=nbin5; i++) binLim5[i]=(Double_t)cTmin  + (cTmax-cTmin)  /nbin5*(Double_t)i ;
 158
 159         // dca
 160         for(Int_t i=0; i<=nbin6; i++) binLim6[i]=(Double_t)dcamin  + (dcamax-dcamin)  /nbin6*(Double_t)i ;
 161
 162         // d0pi
 163         for(Int_t i=0; i<=nbin7; i++) binLim7[i]=(Double_t)d0min  + (d0max-d0min)  /nbin7*(Double_t)i ;
 164
 165         // d0K
 166         for(Int_t i=0; i<=nbin8; i++) binLim8[i]=(Double_t)d0min  + (d0max-d0min)  /nbin8*(Double_t)i ;
 167
 168         // d0xd0
 169         for(Int_t i=0; i<=nbin9; i++) binLim9[i]=(Double_t)d0xd0min  + (d0xd0max-d0xd0min)  /nbin9*(Double_t)i ;
 170
 171         // cosPointingAngle
 172         for(Int_t i=0; i<=nbin10; i++) binLim10[i]=(Double_t)cosmin  + (cosmax-cosmin)  /nbin10*(Double_t)i ;
 173
 174         // Phi
 175         for(Int_t i=0; i<=nbin11; i++) binLim11[i]=(Double_t)phimin  + (phimax-phimin)  /nbin11*(Double_t)i ;
 176
 177         // debugging printings
 178         //Info("AliCFHeavyFlavourTaskMultiVarMultiStep","Printing lower limits for bins in pt");
 179         //for (Int_t i =0; i<= iBin[0]; i++){
 180         //      Info("AliCFHeavyFlavourTaskMultiVarMultiStep",Form("i-th bin, lower limit = %f", binLim0[i]));
 181         //}
 182         //Info("Printing lower limits for bins in ptPi");
 183         //for (Int_t i =0; i<= iBin[3]; i++){
 184         //      Info("AliCFHeavyFlavourTaskMultiVarMultiStep",Form("i-th bin, lower limit = %f", binLim3[i]));
 185         //}
 186         //Info("Printing lower limits for bins in ptKa");
 187         //for (Int_t i =0; i<= iBin[4]; i++){
 188         //      Info("AliCFHeavyFlavourTaskMultiVarMultiStep",Form("i-th bin, lower limit = %f", binLim4[i]));
 189         //      }
 190
 191         //one "container" for MC
 192         AliCFContainer* container = new AliCFContainer("container","container for tracks",nstep,nvar,iBin);
 193         //setting the bin limits
 194         container -> SetBinLimits(ipt,binLim0);
 195         container -> SetBinLimits(iy,binLim1);
 196         container -> SetBinLimits(icosThetaStar,binLim2);
 197         container -> SetBinLimits(ipTpi,binLim3);
 198         container -> SetBinLimits(ipTk,binLim4);
 199         container -> SetBinLimits(icT,binLim5);
 200         container -> SetBinLimits(idca,binLim6);
 201         container -> SetBinLimits(id0pi,binLim7);
 202         container -> SetBinLimits(id0K,binLim8);
 203         container -> SetBinLimits(id0xd0,binLim9);
 204         container -> SetBinLimits(ipointing,binLim10);
 205         container -> SetBinLimits(iphi,binLim11);
 206
 207         //CREATE THE  CUTS -----------------------------------------------
 208
 209         // Gen-Level kinematic cuts
 210         AliCFTrackKineCuts *mcKineCuts = new AliCFTrackKineCuts("mcKineCuts","MC-level kinematic cuts");
 211
 212         //Particle-Level cuts:
 213         AliCFParticleGenCuts* mcGenCuts = new AliCFParticleGenCuts("mcGenCuts","MC particle generation cuts");
 214         mcGenCuts->SetRequirePdgCode(PDG, kTRUE);  // kTRUE set in order to include D0_bar
 215         mcGenCuts->SetAODMC(1); //special flag for reading MC in AOD tree (important)
 216
 217         // Acceptance cuts:
 218         AliCFAcceptanceCuts* accCuts = new AliCFAcceptanceCuts("accCuts", "Acceptance cuts");
 219         AliCFTrackKineCuts *kineAccCuts = new AliCFTrackKineCuts("kineAccCuts","Kine-Acceptance cuts");
 220         kineAccCuts->SetPtRange(ptmin,ptmax);
 221         kineAccCuts->SetEtaRange(etamin,etamax);
 222
 223         // Rec-Level kinematic cuts
 224         AliCFTrackKineCuts *recKineCuts = new AliCFTrackKineCuts("recKineCuts","rec-level kine cuts");
 225
 226         AliCFTrackQualityCuts *recQualityCuts = new AliCFTrackQualityCuts("recQualityCuts","rec-level quality cuts");
 227
 228         AliCFTrackIsPrimaryCuts *recIsPrimaryCuts = new AliCFTrackIsPrimaryCuts("recIsPrimaryCuts","rec-level isPrimary cuts");
 229
 230         printf("CREATE MC KINE CUTS\n");
 231         TObjArray* mcList = new TObjArray(0) ;
 232         mcList->AddLast(mcKineCuts);
 233         mcList->AddLast(mcGenCuts);
 234
 235         printf("CREATE ACCEPTANCE CUTS\n");
 236         TObjArray* accList = new TObjArray(0) ;
 237         accList->AddLast(kineAccCuts);
 238
 239         printf("CREATE RECONSTRUCTION CUTS\n");
 240         TObjArray* recList = new TObjArray(0) ;   // not used!!
 241         recList->AddLast(recKineCuts);
 242         recList->AddLast(recQualityCuts);
 243         recList->AddLast(recIsPrimaryCuts);
 244
 245         TObjArray* emptyList = new TObjArray(0);
 246
 247         //CREATE THE INTERFACE TO CORRECTION FRAMEWORK USED IN THE TASK
 248         printf("CREATE INTERFACE AND CUTS\n");
 249         AliCFManager* man = new AliCFManager() ;
 250         man->SetParticleContainer     (container);
 251         man->SetParticleCutsList(0 , mcList); // MC
 252         man->SetParticleCutsList(1 , accList); // Acceptance
 253         man->SetParticleCutsList(2 , emptyList); // Vertex
 254         man->SetParticleCutsList(3 , emptyList); // Refit
 255         man->SetParticleCutsList(4 , emptyList); // AOD
 256         man->SetParticleCutsList(5 , emptyList); // AOD in Acceptance
 257         man->SetParticleCutsList(6 , emptyList); // AOD with required n. of ITS clusters
 258         man->SetParticleCutsList(7 , emptyList); // AOD Reco (PPR cuts implemented in Task)
 259
 260         // Get the pointer to the existing analysis manager via the static access method.
 261         //==============================================================================
 262         AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 263         if (!mgr) {
 264           ::Error("AddTaskCompareHF", "No analysis manager to connect to.");
 265           return NULL;
 266         }
 267         //CREATE THE TASK
 268         printf("CREATE TASK\n");
 269         // create the task
 270         AliCFHeavyFlavourTaskMultiVarMultiStep *task = new AliCFHeavyFlavourTaskMultiVarMultiStep("AliCFHeavyFlavourTaskMultiVarMultiStep");
 271         task->SetFillFromGenerated(kFALSE);
 272         task->SetMinITSClusters(minITSClusters);
 273         task->SetCFManager(man); //here is set the CF manager
 274
 275         //-----------------------------------------------------------//
 276         //   create correlation matrix for unfolding - only eta-pt   //
 277         //-----------------------------------------------------------//
 278
 279         Bool_t AcceptanceUnf = kTRUE; // unfold at acceptance level, otherwise PPR
 280
 281         Int_t thnDim[4];
 282
 283         //first half  : reconstructed
 284         //second half : MC
 285
 286         thnDim[0] = iBin[0];
 287         thnDim[2] = iBin[0];
 288         thnDim[1] = iBin[1];
 289         thnDim[3] = iBin[1];
 290
 291         THnSparseD* correlation = new THnSparseD("correlation","THnSparse with correlations",4,thnDim);
 292         Double_t** binEdges = new Double_t[2];
 293
 294         // set bin limits
 295
 296         binEdges[0]= binLim0;
 297         binEdges[1]= binLim1;
 298
 299         correlation->SetBinEdges(0,binEdges[0]);
 300         correlation->SetBinEdges(2,binEdges[0]);
 301
 302         correlation->SetBinEdges(1,binEdges[1]);
 303         correlation->SetBinEdges(3,binEdges[1]);
 304
 305         correlation->Sumw2();
 306
 307         // correlation matrix ready
 308         //------------------------------------------------//
 309
 310         task->SetCorrelationMatrix(correlation); // correlation matrix for unfolding
 311
 312         // Create and connect containers for input/output
 313
 314         // ------ input data ------
 315         AliAnalysisDataContainer *cinput0  = mgr->GetCommonInputContainer();
 316
 317         // ----- output data -----
 318
 319         //slot 0 : default output tree (by default handled by AliAnalysisTaskSE)
 320         AliAnalysisDataContainer *coutput0 = mgr->CreateContainer("ctree0", TTree::Class(),AliAnalysisManager::kOutputContainer,"output.root");
 321
 322         //now comes user's output objects :
 323
 324         // output TH1I for event counting
 325         AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("chist0", TH1I::Class(),AliAnalysisManager::kOutputContainer,"output.root");
 326         // output Correction Framework Container (for acceptance & efficiency calculations)
 327         AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("ccontainer0", AliCFContainer::Class(),AliAnalysisManager::kOutputContainer,"output.root");
 328 // Unfolding - correlation matrix
 329         AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("corr0", THnSparseD::Class(),AliAnalysisManager::kOutputContainer,"output.root");
 330
 331         mgr->AddTask(task);
 332
 333         mgr->ConnectInput(task,0,mgr->GetCommonInputContainer());
 334         mgr->ConnectOutput(task,0,coutput0);
 335         mgr->ConnectOutput(task,1,coutput1);
 336         mgr->ConnectOutput(task,2,coutput2);
 337         mgr->ConnectOutput(task,3,coutput3);
 338         return task;
 339 }
 340