PWG3/vertexingHF/macros/AddTaskCFDStar.C

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