PWGHF/correlationHF/AliAnalysisTaskSEmcCorr.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2009, 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 /////////////////////////////////////////////////////////////
  17 //
  18 // Class AliAnalysisTaskSEmcCorr
  19 // AliAnalysisTaskSE for studying HF-(hadron,electrons) and hadron-hadron correlations
  20 //     at MC level
  21 // Author: Andrea Rossi, andrea.rossi@cern.ch
  22 /////////////////////////////////////////////////////////////
  23
  24
  25 #include <TH1F.h>
  26 #include <TH2F.h>
  27 #include <TAxis.h>
  28 #include <THnSparse.h>
  29 #include <TDatabasePDG.h>
  30 #include <TMath.h>
  31 #include <TROOT.h>
  32 #include <TArrayI.h>
  33 #include <TString.h>
  34 #include "AliAODEvent.h"
  35 #include "AliAODRecoDecayHF2Prong.h"
  36 #include "AliAODRecoDecayHF.h"
  37 #include "AliAODRecoDecay.h"
  38 #include "AliAnalysisDataSlot.h"
  39 #include "AliAnalysisDataContainer.h"
  40 #include "AliAODTrack.h"
  41 #include "AliAODHandler.h"
  42 #include "AliESDtrack.h"
  43 #include "AliAODVertex.h"
  44 #include "AliESDVertex.h"
  45 #include "AliVertexerTracks.h"
  46 #include "AliAODMCParticle.h"
  47 #include "AliAODPid.h"
  48 #include "AliTPCPIDResponse.h"
  49 #include "AliAODMCHeader.h"
  50 #include "AliAnalysisVertexingHF.h"
  51 #include "AliRDHFCutsD0toKpi.h"
  52 #include "AliAODInputHandler.h"
  53 #include "AliAnalysisManager.h"
  54 #include "AliNormalizationCounter.h"
  55 //#include "/Users/administrator/ALICE/CHARM/HFCJ/TestsForProduction/AliAnalysisTaskSEmcCorr.h"
  56 #include "AliAnalysisTaskSEmcCorr.h"
  57 #include "AliGenEventHeader.h"
  58
  59 class TCanvas;
  60 class TTree;
  61 class TChain;
  62 class AliAnalysisTaskSE;
  63
  64
  65 ClassImp(AliAnalysisTaskSEmcCorr)
  66
  67 AliAnalysisTaskSEmcCorr::AliAnalysisTaskSEmcCorr()
  68 : AliAnalysisTaskSE(),
  69   fReadMC(kTRUE),
  70   fFastSimul(kFALSE),
  71   fCheckDecay(kTRUE),
  72   fDoHFCorrelations(kTRUE),
  73   fOnlyKine(kTRUE),
  74   fDoHadronHadron(kFALSE),
  75   fNentries(0),
  76   fhNprongsD0(0),
  77   fhNprongsD0chargedOnly(0),
  78   fhNprongsD0chargedRef(0),
  79   fYrangeTrig(0.5),
  80   fEtarangeEleTrig(0.8),
  81   fminDpt(1.),
  82   fArrayDaugh(0x0),
  83   flastdaugh(0),
  84   fArrayAssoc(0x0),
  85   fLastAss(0),
  86   fminPtass(0.3),
  87   fmaxEtaAss(99.),
  88   fhMCcorrel(0x0),
  89   fhMCtrigPart(0x0),
  90   fhMChadroncorrel(0x0),
  91   fhMChadrontrigPart(0x0),
  92   fhDzeroDecay(0x0),
  93   fhDplusDecay(0x0),
  94   fhLambdaCDecay(0x0),
  95   fhAllBDecay(0x0),
  96   fnpionPlus(0),
  97   fnpionMinus(0),
  98   fnpionZero(0),
  99   fnkaonPlus(0),
 100   fnkaonMinus(0),
 101   fnkaonZeroShort(0),
 102   fnProton(0),
 103   fnAntiProton(0),
 104   fnElePlus(0),
 105   fnEleMinus(0),
 106   fnMuonPlus(0),
 107   fnMuonMinus(0),
 108   fnNeutrino(0),
 109   fnJpsi(0),
 110   fnPhysPrim(0),
 111   fpxtotdaugh(0),
 112   fpytotdaugh(0),
 113   fpztotdaugh(0),
 114   fGeneratorString(0)
 115 {// standard constructor
 116
 117 }
 118
 119
 120
 121 //________________________________________________________________________
 122 AliAnalysisTaskSEmcCorr::AliAnalysisTaskSEmcCorr(const char *name)
 123   : AliAnalysisTaskSE(name),
 124     fReadMC(kTRUE),
 125     fFastSimul(kFALSE),
 126     fCheckDecay(kTRUE),
 127     fDoHFCorrelations(kTRUE),
 128     fOnlyKine(kTRUE),
 129     fDoHadronHadron(kFALSE),
 130     fNentries(0),
 131     fhNprongsD0(0),
 132     fhNprongsD0chargedOnly(0),
 133     fhNprongsD0chargedRef(0),
 134     fYrangeTrig(0.5),
 135     fEtarangeEleTrig(0.8),
 136     fminDpt(1.),
 137     fArrayDaugh(0x0),
 138     flastdaugh(0),
 139     fArrayAssoc(0x0),
 140     fLastAss(0),
 141     fminPtass(0.3),
 142     fmaxEtaAss(99.),
 143     fhMCcorrel(0x0),
 144     fhMCtrigPart(0x0),
 145     fhMChadroncorrel(0x0),
 146     fhMChadrontrigPart(0x0),
 147     fhDzeroDecay(0x0),
 148     fhDplusDecay(0x0),
 149     fhLambdaCDecay(0x0),
 150     fhAllBDecay(0x0),
 151     fnpionPlus(0),
 152     fnpionMinus(0),
 153     fnpionZero(0),
 154     fnkaonPlus(0),
 155     fnkaonMinus(0),
 156     fnkaonZeroShort(0),
 157     fnProton(0),
 158     fnAntiProton(0),
 159     fnElePlus(0),
 160     fnEleMinus(0),
 161     fnMuonPlus(0),
 162     fnMuonMinus(0),
 163     fnNeutrino(0),
 164     fnJpsi(0),
 165     fnPhysPrim(0),
 166     fpxtotdaugh(0),
 167     fpytotdaugh(0),
 168     fpztotdaugh(0),
 169     fGeneratorString(0)
 170 { // default constructor
 171
 172   DefineOutput(1, TH1F::Class());
 173   DefineOutput(2, TH1F::Class());
 174   DefineOutput(3, TH1F::Class());
 175   DefineOutput(4, TH1F::Class());
 176   DefineOutput(5, THnSparseF::Class());
 177   DefineOutput(6, THnSparseF::Class());
 178   DefineOutput(7, THnSparseF::Class());
 179   DefineOutput(8, THnSparseF::Class());
 180   DefineOutput(9,TH1F::Class());
 181   DefineOutput(10,TH1F::Class());
 182   DefineOutput(11,TH1F::Class());
 183   DefineOutput(12,TH1F::Class());
 184 }
 185
 186 //________________________________________________________________________
 187 AliAnalysisTaskSEmcCorr::~AliAnalysisTaskSEmcCorr(){
 188   // destructor
 189   delete fNentries;
 190   delete fhNprongsD0;
 191   delete fhNprongsD0chargedOnly;
 192   delete fhNprongsD0chargedRef;
 193   delete fhMCcorrel;
 194   delete fhMCtrigPart;
 195   delete fhMChadroncorrel;
 196   delete fhMChadrontrigPart;
 197   delete fArrayDaugh;
 198   delete fArrayAssoc;
 199   delete fhDzeroDecay;
 200   delete fhDplusDecay;
 201   delete fhLambdaCDecay;
 202   delete fhAllBDecay;
 203 }
 204
 205 //________________________________________________________________________
 206 void AliAnalysisTaskSEmcCorr::Init()
 207 {
 208   // Initialization
 209   fArrayDaugh=new TArrayI(20);
 210
 211   fArrayAssoc=new TArrayI(200);
 212
 213 }
 214 //_______________________________________________
 215 void AliAnalysisTaskSEmcCorr::GetFinalStateParticles(AliAODMCParticle *part,TClonesArray *clarr){
 216
 217   Int_t fdaugh=part->GetDaughter(0);
 218   Int_t ldaugh=part->GetDaughter(1);
 219   if(fdaugh<0||ldaugh<0||(fdaugh>ldaugh))return;
 220   for(Int_t j=fdaugh;j<=ldaugh;j++){
 221     AliAODMCParticle *pd=(AliAODMCParticle*)clarr->At(j);
 222     Int_t pdg=pd->GetPdgCode();
 223     fnPhysPrim++;// increment it now and decrement it after all if, in case the daughter is not a phys prim
 224     fpxtotdaugh+=pd->Px(); // increment it now and decrement it after all if, in case the daughter is not a phys prim
 225     fpytotdaugh+=pd->Py();
 226     fpztotdaugh+=pd->Pz();
 227
 228     if(pdg==211)fnpionPlus++;
 229     else if(pdg==-211)fnpionMinus++;
 230     else if(pdg==111)fnpionZero++;
 231     else if(pdg==321)fnkaonPlus++;
 232     else if(pdg==-321)fnkaonMinus++;
 233     else if(pdg==310)fnkaonZeroShort++;
 234     else if(pdg==2212)fnProton++;
 235     else if(pdg==-2212)fnAntiProton++;
 236     else if(pdg==-11)fnElePlus++;
 237     else if(pdg==11)fnEleMinus++;
 238     else if(pdg==-13)fnMuonPlus++;
 239     else if(pdg==13)fnMuonMinus++;
 240     else if(pdg==12||pdg==14||pdg==-12||pdg==-14)fnNeutrino++;
 241     else if(pdg==443)fnJpsi++;
 242     else if(pdg!=130&&pdg!=22){
 243       fnPhysPrim--;
 244       fpxtotdaugh-=pd->Px();
 245       fpytotdaugh-=pd->Py();
 246       fpztotdaugh-=pd->Pz();
 247       GetFinalStateParticles(pd,clarr);
 248     }
 249   }
 250 }
 251
 252 //_______________________________________________
 253 void AliAnalysisTaskSEmcCorr::CheckDzeroChannels(AliAODMCParticle *part,TClonesArray *clarr){
 254   Int_t pdg=part->GetPdgCode();
 255   if(!(TMath::Abs(pdg)==421))return;
 256   fhDzeroDecay->Fill(0);
 257   fnpionPlus=0;
 258   fnpionMinus=0;
 259   fnpionZero=0;
 260   fnkaonPlus=0;
 261   fnkaonMinus=0;
 262   fnkaonZeroShort=0;
 263   fnProton=0;
 264   fnAntiProton=0;
 265   fnElePlus=0;
 266   fnEleMinus=0;
 267   fnMuonPlus=0;
 268   fnMuonMinus=0;
 269   fnNeutrino=0;
 270   fnJpsi=0;
 271   fnPhysPrim=0;
 272   fpxtotdaugh=0.;
 273   fpytotdaugh=0.;
 274   fpztotdaugh=0.;
 275
 276   GetFinalStateParticles(part,clarr);
 277
 278   Double_t px=part->Px();
 279   Double_t py=part->Py();
 280   Double_t pz=part->Pz();
 281   // check exclusive channels in the list
 282   if(TMath::Abs(px-fpxtotdaugh)<1.e-6&&TMath::Abs(py-fpytotdaugh)<1.e-6&&TMath::Abs(pz-fpztotdaugh)<1.e-6){
 283     if(fnPhysPrim==2&&((pdg==421&&fnpionPlus==1&&fnkaonMinus==1)||(pdg==-421&&fnpionMinus==1&&fnkaonPlus==1))){
 284       fhDzeroDecay->Fill(6);
 285     }
 286     if(fnPhysPrim==3&&((pdg==421&&fnpionPlus==1&&fnkaonMinus==1&&fnpionZero==1)||(pdg==-421&&fnpionMinus==1&&fnkaonPlus==1&&fnpionZero==1))){
 287       fhDzeroDecay->Fill(7);
 288     }
 289     if(fnPhysPrim==4&&((pdg==421&&fnpionPlus==2&&fnkaonMinus==1&&fnpionMinus==1)||(pdg==-421&&fnpionMinus==2&&fnkaonPlus==1&&fnpionPlus==1))){
 290       fhDzeroDecay->Fill(8);
 291     }
 292     if(fnPhysPrim==3&&((pdg==421&&fnElePlus==1&&fnkaonMinus==1&&fnNeutrino==1)||(pdg==-421&&fnEleMinus==1&&fnkaonPlus==1&&fnNeutrino==1))){
 293       fhDzeroDecay->Fill(9);
 294     }
 295     if(fnPhysPrim==3&&((pdg==421&&fnMuonPlus==1&&fnkaonMinus==1&&fnNeutrino==1)||(pdg==-421&&fnMuonMinus==1&&fnkaonPlus==1&&fnNeutrino==1))){
 296       fhDzeroDecay->Fill(10);
 297     }
 298   }
 299
 300   //now check decay to electron + X
 301   if(fnElePlus>=1||fnEleMinus>=1){
 302     fhDzeroDecay->Fill(11);
 303   }
 304   // no check number of prong
 305   Int_t nprong=fnPhysPrim;
 306   nprong-=fnpionZero;
 307   nprong-=fnNeutrino;
 308   nprong-=fnkaonZeroShort;
 309   if(nprong==0){
 310     fhDzeroDecay->Fill(1);
 311   }
 312   else if(nprong==1){// should never happen for D0
 313     fhDzeroDecay->Fill(2);
 314   }
 315   else if(nprong==2){// should never happen for D0
 316     fhDzeroDecay->Fill(3);
 317   }
 318   else if(nprong==3){
 319     fhDzeroDecay->Fill(4);
 320   }
 321   else if(nprong==4){
 322     fhDzeroDecay->Fill(5);
 323   }
 324
 325   return;
 326 }
 327
 328
 329 //_______________________________________________
 330 void AliAnalysisTaskSEmcCorr::CheckDplusChannels(AliAODMCParticle *part,TClonesArray *clarr){
 331   Int_t pdg=part->GetPdgCode();
 332   if(!(TMath::Abs(pdg)==411))return;
 333   fhDplusDecay->Fill(0);
 334   fnpionPlus=0;
 335   fnpionMinus=0;
 336   fnpionZero=0;
 337   fnkaonPlus=0;
 338   fnkaonMinus=0;
 339   fnkaonZeroShort=0;
 340   fnProton=0;
 341   fnAntiProton=0;
 342   fnElePlus=0;
 343   fnEleMinus=0;
 344   fnMuonPlus=0;
 345   fnMuonMinus=0;
 346   fnNeutrino=0;
 347   fnJpsi=0;
 348   fnPhysPrim=0;
 349   fpxtotdaugh=0.;
 350   fpytotdaugh=0.;
 351   fpztotdaugh=0.;
 352   GetFinalStateParticles(part,clarr);
 353
 354   Double_t px=part->Px();
 355   Double_t py=part->Py();
 356   Double_t pz=part->Pz();
 357   // check exclusive channels in the list
 358   if(TMath::Abs(px-fpxtotdaugh)<1.e-6&&TMath::Abs(py-fpytotdaugh)<1.e-6&&TMath::Abs(pz-fpztotdaugh)<1.e-6){
 359     if(fnPhysPrim==3&&((pdg==411&&fnpionPlus==2&&fnkaonMinus==1)||(pdg==-411&&fnpionMinus==2&&fnkaonPlus==1))){
 360       fhDplusDecay->Fill(6);
 361     }
 362     if(fnPhysPrim==3&&((pdg==411&&fnkaonZeroShort==1&&fnpionPlus==1&&fnpionZero==1)||(pdg==-411&&fnkaonZeroShort==1&&fnpionMinus==1&&fnpionZero==1))){
 363       fhDplusDecay->Fill(7);
 364     }
 365     if(fnPhysPrim==4&&((pdg==411&&fnpionPlus==2&&fnkaonMinus==1&&fnpionZero==1)||(pdg==-411&&fnpionMinus==2&&fnkaonPlus==1&&fnpionZero==1))){
 366       fhDzeroDecay->Fill(8);
 367     }
 368     if(fnPhysPrim==3&&((pdg==411&&fnElePlus==1&&fnkaonZeroShort==1&&fnNeutrino==1)||(pdg==-411&&fnEleMinus==1&&fnkaonZeroShort==1&&fnNeutrino==1))){
 369       fhDplusDecay->Fill(9);
 370     }
 371     if(fnPhysPrim==3&&((pdg==411&&fnMuonPlus==1&&fnkaonZeroShort==1&&fnNeutrino==1)||(pdg==-411&&fnMuonMinus==1&&fnkaonZeroShort==1&&fnNeutrino==1))){
 372       fhDplusDecay->Fill(10);
 373     }
 374   }
 375
 376   //now check decay to electron + X
 377   if(fnElePlus>=1||fnEleMinus>=1){
 378     fhDplusDecay->Fill(11);
 379   }
 380   // no check number of prong
 381   Int_t nprong=fnPhysPrim;
 382   nprong-=fnpionZero;
 383   nprong-=fnNeutrino;
 384   nprong-=fnkaonZeroShort;
 385   if(nprong==0){
 386     fhDplusDecay->Fill(1);
 387   }
 388   else if(nprong==1){// should never happen for D0
 389     fhDplusDecay->Fill(2);
 390   }
 391   else if(nprong==2){// should never happen for D0
 392     fhDplusDecay->Fill(3);
 393   }
 394   else if(nprong==3){
 395     fhDplusDecay->Fill(4);
 396   }
 397   else if(nprong==4){
 398     fhDplusDecay->Fill(5);
 399   }
 400
 401   return;
 402 }
 403 //________________________________________________________________________
 404 // void AliAnalysisTaskSEmcCorr::CheckLambdaChannels(AliAODMCParticle *part,TClonesArray *clarr){
 405 //   // not implemented yet
 406 // }
 407
 408 //________________________________________________________________________
 409 // void AliAnalysisTaskSEmcCorr::CheckBallChannels(AliAODMCParticle *part,TClonesArray *clarr){
 410 //   // not implemented yet
 411 // }
 412
 413
 414 //________________________________________________________________________
 415 void AliAnalysisTaskSEmcCorr::UserCreateOutputObjects(){
 416
 417   if(!fArrayDaugh)  fArrayDaugh=new TArrayI(20);
 418
 419   if(!fArrayAssoc)fArrayAssoc=new TArrayI(200);
 420
 421
 422   fNentries=new TH1F("nentriesChFr", "Analyzed sample properties", 21,0.5,21.5);
 423   fNentries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
 424   fNentries->GetXaxis()->SetBinLabel(2,"nEvSel");
 425   fNentries->GetXaxis()->SetBinLabel(3,"nEvPile-up Rej");
 426   fNentries->GetXaxis()->SetBinLabel(4,"nEvGoodVtxS");
 427   fNentries->GetXaxis()->SetBinLabel(5,"nEvRejVtxZ");
 428   fNentries->GetXaxis()->SetBinLabel(6,"nD0");
 429   fNentries->GetXaxis()->SetBinLabel(7,"nEventsSelGenName");
 430
 431   fhNprongsD0=new TH1F("fhNprongsD0","fhNprongsD0",20,-0.5,19.5);
 432   fhNprongsD0chargedOnly=new TH1F("fhNprongsD0chargedOnly","fhNprongsD0chargedOnly",20,-0.5,19.5);
 433   fhNprongsD0chargedRef=new TH1F("fhNprongsD0chargedRef","fhNprongsD0chargedRef",20,-0.5,19.5);
 434
 435
 436   Int_t nbinsCorrMC[8]={15,20,20,20,50,63,20,11};
 437   Double_t binlowCorrMC[8]={-7.5,0.,-1.,0.,-5.,-0.5*TMath::Pi(),-2,-1.5};
 438   Double_t binupCorrMC[8]={7.5,20.,1.,2.,5.,1.5*TMath::Pi(),2.,9.5};
 439   fhMCcorrel=new THnSparseF("fhMCcorrel","fhMCcorrel;pdg;ptTrig;etaTrig;ptAss;etaAss;deltaPhi;deltaEta;pdgAss;",8,nbinsCorrMC,binlowCorrMC,binupCorrMC);
 440
 441   Int_t nbinsTrigMC[3]={15,20,20};
 442   Double_t binlowTrigMC[3]={-7.5,0.,-1.};
 443   Double_t binupTrigMC[3]={7.5,20.,1.};
 444   fhMCtrigPart=new THnSparseF("fhMCtrigPart","fhMCcorrel;pdg;ptTrig;etaTrig;",3,nbinsTrigMC,binlowTrigMC,binupTrigMC);
 445
 446   fhDzeroDecay=new TH1F("fhDzeroDecay","Dzero Decay channels",14,-0.5,13.5);
 447   fhDzeroDecay->GetXaxis()->SetBinLabel(1,"All D0");
 448   fhDzeroDecay->GetXaxis()->SetBinLabel(2,"0 prong");
 449   fhDzeroDecay->GetXaxis()->SetBinLabel(3,"1 prong");
 450   fhDzeroDecay->GetXaxis()->SetBinLabel(4,"2 prong");
 451   fhDzeroDecay->GetXaxis()->SetBinLabel(5,"3 prong");
 452   fhDzeroDecay->GetXaxis()->SetBinLabel(6,"4 prong");
 453   fhDzeroDecay->GetXaxis()->SetBinLabel(7,"K pi (3.88)");
 454   fhDzeroDecay->GetXaxis()->SetBinLabel(8,"K pi pi0 (13.9)");
 455   fhDzeroDecay->GetXaxis()->SetBinLabel(9,"K pi pi pi (8.09)");
 456   fhDzeroDecay->GetXaxis()->SetBinLabel(10,"e K nu (3.55)");
 457   fhDzeroDecay->GetXaxis()->SetBinLabel(11,"mu K nu (3.31)");
 458   fhDzeroDecay->GetXaxis()->SetBinLabel(12,"e X (6.49)");
 459
 460
 461   fhDplusDecay=new TH1F("fhDplusDecay","Dplus Decay channels",14,-0.5,13.5);
 462   fhDplusDecay->GetXaxis()->SetBinLabel(1,"All D+");
 463   fhDplusDecay->GetXaxis()->SetBinLabel(2,"0 prong");
 464   fhDplusDecay->GetXaxis()->SetBinLabel(3,"1 prong");
 465   fhDplusDecay->GetXaxis()->SetBinLabel(4,"2 prong");
 466   fhDplusDecay->GetXaxis()->SetBinLabel(5,"3 prong");
 467   fhDplusDecay->GetXaxis()->SetBinLabel(6,"4 prong");
 468   fhDplusDecay->GetXaxis()->SetBinLabel(7,"K pi pi (9.4)");
 469   fhDplusDecay->GetXaxis()->SetBinLabel(8,"K0s pi pi0 (6.9)");
 470   fhDplusDecay->GetXaxis()->SetBinLabel(9,"K pi pi pi0 (6.08)");
 471   fhDplusDecay->GetXaxis()->SetBinLabel(10,"e K0bar nu (8.83)");
 472   fhDplusDecay->GetXaxis()->SetBinLabel(11,"mu K0bar nu (9.4)");
 473   fhDplusDecay->GetXaxis()->SetBinLabel(12,"e X (16.07");
 474
 475   fhLambdaCDecay=new TH1F("fhLambdaCDecay","LambdaC Decay channels",10,-0.5,9.5);
 476   fhLambdaCDecay->GetXaxis()->SetBinLabel(1,"All LambdaC");
 477   fhLambdaCDecay->GetXaxis()->SetBinLabel(2,"p K0s");
 478   fhLambdaCDecay->GetXaxis()->SetBinLabel(3,"p K pi");
 479
 480   fhAllBDecay=new TH1F("fhAllBDecay","All bhadron (also LambdaB) Decay channels",10,-0.5,9.5);
 481   fhAllBDecay->GetXaxis()->SetBinLabel(1,"All b-hadrons");
 482   fhAllBDecay->GetXaxis()->SetBinLabel(2,"D0bar X (59.6)");
 483   fhAllBDecay->GetXaxis()->SetBinLabel(3,"D0bar D0 X (5.1) ");
 484   fhAllBDecay->GetXaxis()->SetBinLabel(4,"D+- D0 X (2.7)");
 485   fhAllBDecay->GetXaxis()->SetBinLabel(5,"D- X (22.7)");
 486   fhAllBDecay->GetXaxis()->SetBinLabel(6,"J/Psi X (1.16)");
 487
 488
 489
 490   Int_t nbinsCorrMChadron[8]={6,20,20,20,50,63,20,7};
 491   Double_t binlowCorrMChadron[8]={-1.5,0.,-1.,0.,-5.,-0.5*TMath::Pi(),-2,-1.5};
 492   Double_t binupCorrMChadron[8]={4.5,20.,1.,2.,5.,1.5*TMath::Pi(),2.,5.5};
 493   fhMChadroncorrel=new THnSparseF("fhMChadroncorrel","fhMChadroncorrel;pdg;ptTrig;etaTrig;ptAss;etaAss;deltaPhi;deltaEta;pdgAss;",8,nbinsCorrMChadron,binlowCorrMChadron,binupCorrMChadron);
 494
 495   Int_t nbinsTrigMChadron[3]={6,20,20};
 496   Double_t binlowTrigMChadron[3]={-1.5,0.,-1.};
 497   Double_t binupTrigMChadron[3]={4.5,20.,1.};
 498   fhMChadrontrigPart=new THnSparseF("fhMChadrontrigPart","fhMChadroncorrel;pdg;ptTrig;etaTrig;",3,nbinsTrigMChadron,binlowTrigMChadron,binupTrigMChadron);
 499
 500
 501   PostData(1,fNentries);
 502   PostData(2,fhNprongsD0);
 503   PostData(3,fhNprongsD0chargedOnly);
 504   PostData(4,fhNprongsD0chargedRef);
 505   PostData(5,fhMCcorrel);
 506   PostData(6,fhMCtrigPart);
 507   PostData(7,fhMChadroncorrel);
 508   PostData(8,fhMChadrontrigPart);
 509   PostData(9,fhDzeroDecay);
 510   PostData(10,fhDplusDecay);
 511   PostData(11,fhLambdaCDecay);
 512   PostData(12,fhAllBDecay);
 513
 514
 515 }
 516
 517
 518
 519 //________________________________________________________________________
 520 void AliAnalysisTaskSEmcCorr::UserExec(Option_t */*option*/){
 521
 522
 523
 524   // Execute analysis for current event:
 525   // heavy flavor candidates association to MC truth
 526
 527   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 528   if (!aod) {
 529     Printf("ERROR: aod not available");
 530     return;
 531   }
 532
 533
 534
 535   fNentries->Fill(1);
 536
 537   // AOD primary vertex
 538
 539
 540   if(!fOnlyKine){
 541     AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 542     TString primTitle = vtx1->GetTitle();
 543     if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
 544
 545       fNentries->Fill(3);
 546
 547     }
 548     else {
 549       PostData(1,fNentries);
 550       return;
 551
 552     }
 553   }
 554
 555
 556   TClonesArray *arrayMC=0x0;
 557   AliAODMCHeader *aodmcHeader=0x0;
 558   Double_t vtxTrue[3];
 559
 560
 561   if(fReadMC){
 562     //    aod->GetList()->Print();
 563     // load MC particles
 564     if(fFastSimul){
 565       arrayMC =
 566         (TClonesArray*)aod->GetList()->At(23);//FindObject("mcparticles");
 567     }
 568     else{
 569       arrayMC =
 570         (TClonesArray*)aod->GetList()->FindObject("mcparticles");
 571     }
 572     if(!arrayMC) {
 573       Printf("AliAnalysisTaskSEmcCorr::UserExec: MC particles branch not found!\n");
 574       fNentries->Fill(2);
 575       PostData(1,fNentries);
 576       return;
 577     }
 578     // load MC header
 579     aodmcHeader =
 580       (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 581     if(!aodmcHeader) {
 582       Printf("AliAnalysisTaskSEmcCorr::UserExec: MC header branch not found!\n");
 583       fNentries->Fill(3);
 584       PostData(1,fNentries);
 585       return;
 586     }
 587
 588     if(!(fGeneratorString.IsNull())){
 589       TString strGenTitle=((AliGenEventHeader*)aodmcHeader->GetCocktailHeader(0))->GetTitle();
 590       // Printf("Gen Title: %s",strGenTitle.Data());
 591       if(!fGeneratorString.Contains(strGenTitle.Data())){
 592         //      Printf("Event rejected from generator title");
 593         return;
 594       }
 595     }
 596     fNentries->Fill(4);
 597     //    printf("N MC entries: %d \n",arrayMC->GetEntriesFast());
 598     // MC primary vertex
 599     aodmcHeader->GetVertex(vtxTrue);
 600
 601
 602     // FILL HISTOS FOR D0 mesons, c quarks and D0 from B properties
 603     if(fDoHFCorrelations)      SelectAssociatedParticles(arrayMC);
 604     for(Int_t jp=0;jp<arrayMC->GetEntriesFast();jp++){
 605       AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->At(jp);
 606       Int_t pdg=TMath::Abs(part->GetPdgCode());
 607       if(fDoHFCorrelations){
 608         if(pdg==421||pdg==411||pdg==413){
 609
 610           if(TMath::Abs(part->Y())>fYrangeTrig)continue;
 611           if(TMath::Abs(part->Pt())<fminDpt)continue;
 612           FillSkipParticleArray(part,arrayMC,jp);
 613           FillCorrelationPlots(part,arrayMC);
 614         }
 615         else if(pdg==11){
 616           if(TMath::Abs(part->Eta())>fEtarangeEleTrig)continue;
 617           if(TMath::Abs(part->Pt())<fminDpt)continue;
 618           FillSkipParticleArray(part,arrayMC,jp);
 619           FillCorrelationPlots(part,arrayMC);
 620         }
 621       }
 622       if(fCheckDecay){
 623         CheckDzeroChannels(part,arrayMC);
 624         CheckDplusChannels(part,arrayMC);
 625         //        CheckLambdaChannels(part,arrayMC);
 626         //        CheckBallChannels(part,arrayMC);
 627
 628         // NOW STUDY OF N-PRONGS (OLD PART OF THE TASK)
 629         if(pdg==421){
 630           Bool_t semilept=kFALSE;
 631           Int_t nprongsD0=0;
 632           Int_t nprongsD0charged=0;
 633
 634           Int_t lab1=part->GetDaughter(0);
 635           Int_t lab2=part->GetDaughter(1);
 636           if(lab1>=0&&lab2>=0){
 637             for(Int_t jd=lab1;jd<=lab2;jd++){
 638               AliAODMCParticle *d=(AliAODMCParticle*)arrayMC->At(jd);
 639               Int_t pdgd=TMath::Abs(d->GetPdgCode());
 640               if(pdgd>=11&&pdgd<=16){
 641                 semilept=kTRUE;
 642                 break;
 643               }
 644               if(pdgd==211||pdgd==111||pdgd==321||pdgd==311||pdgd==310||pdgd==130){
 645                 nprongsD0++;
 646               }
 647               if(pdgd==211||pdgd==321){
 648                 nprongsD0charged++;
 649               }
 650               else{
 651                 Int_t lab1d=d->GetDaughter(0);
 652                 Int_t lab2d=d->GetDaughter(1);
 653                 if(lab1d>=0&&lab2d>=0){
 654                   for(Int_t jdd=lab1d;jdd<=lab2d;jdd++){
 655                     AliAODMCParticle *dd=(AliAODMCParticle*)arrayMC->At(jdd);
 656                     Int_t pdgdd=TMath::Abs(dd->GetPdgCode());
 657                     if(pdgdd==211||pdgdd==111||pdgdd==321||pdgdd==311||pdgdd==310||pdgdd==130){
 658                       nprongsD0++;
 659                     }
 660                     if(pdgd==211||pdgd==321){
 661                       nprongsD0charged++;
 662                     }
 663                   }
 664                 }
 665               }
 666             }
 667           }
 668           if(!semilept){
 669             fhNprongsD0->Fill(nprongsD0);
 670             fhNprongsD0chargedOnly->Fill(nprongsD0charged);
 671             fhNprongsD0chargedRef->Fill(nprongsD0charged);
 672           }
 673         }
 674       }
 675     }
 676
 677     if(fDoHadronHadron){
 678       for(Int_t jp=0;jp<fLastAss;jp++){
 679
 680         AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->At(fArrayAssoc->At(jp));
 681         Int_t pdg=TMath::Abs(part->GetPdgCode());
 682         if(pdg==11||pdg==13||pdg==211||pdg==2212||pdg==321){
 683           if(TMath::Abs(part->Eta())>fEtarangeEleTrig)continue;
 684           if(TMath::Abs(part->Pt())<fminDpt)continue;
 685           FillSkipParticleArray(part,arrayMC,fArrayAssoc->At(jp));
 686           FillCorrelationPlotsHadrons(part,arrayMC);
 687         }
 688       }
 689
 690     }
 691
 692   }
 693
 694
 695   PostData(1,fNentries);
 696   PostData(2,fhNprongsD0);
 697   PostData(3,fhNprongsD0chargedOnly);
 698   PostData(4,fhNprongsD0chargedRef);
 699   PostData(5,fhMCcorrel);
 700   PostData(6,fhMCtrigPart);
 701   PostData(7,fhMChadroncorrel);
 702   PostData(8,fhMChadrontrigPart);
 703   PostData(9,fhDzeroDecay);
 704   PostData(10,fhDplusDecay);
 705   PostData(11,fhLambdaCDecay);
 706   PostData(12,fhAllBDecay);
 707 }
 708
 709
 710
 711 void AliAnalysisTaskSEmcCorr::FillSkipParticleArray(AliAODMCParticle *part,TClonesArray *arrayMC,Int_t partPos){
 712   fArrayDaugh->Reset(0);
 713   fArrayDaugh->AddAt(partPos,0);
 714   flastdaugh=1;
 715
 716   // Loop on daughters and nephews
 717   if(part->GetNDaughters()>0){
 718     Int_t pos1=part->GetDaughter(0);
 719     Int_t pos2=part->GetDaughter(1);
 720     if(pos2<0)pos2=pos1;
 721     if(pos1>=0){
 722       for(Int_t jd=pos1;jd<=pos2;jd++){
 723         AliAODMCParticle *pd=(AliAODMCParticle*)arrayMC->At(jd);
 724         Int_t pdgd=TMath::Abs(pd->GetPdgCode());
 725         if(pdgd==11||pdgd==13||pdgd==211||pdgd==321||pdgd==2212){
 726           fArrayDaugh->AddAt(jd,flastdaugh);
 727           flastdaugh++;
 728         }
 729         else if(pdgd!=12&&pdgd!=14&&pdgd!=16&&pdgd!=310&&pdgd!=111){// CHECK NEPHEWS (FOR RESONANT CHANNELS)
 730           Int_t nephw=pd->GetNDaughters();
 731           if(nephw>0){
 732             Int_t posN1=pd->GetDaughter(0);
 733             Int_t posN2=pd->GetDaughter(1);
 734             if(posN2<0)posN2=posN1;
 735             if(posN1>=0){
 736               for(Int_t jN=posN1;jN<=posN2;jN++){
 737                 AliAODMCParticle *pN=(AliAODMCParticle*)arrayMC->At(jN);
 738                 Int_t pdgN=TMath::Abs(pN->GetPdgCode());
 739                 if(pdgN==11||pdgN==13||pdgN==211||pdgN==321||pdgN==2212){
 740                   fArrayDaugh->AddAt(jN,flastdaugh);
 741                   flastdaugh++;
 742                 }
 743                 else if(pdgN!=12&&pdgN!=14&&pdgN!=16&&pdgN!=310&&pdgN!=111){// CHECK one step more (FOR RESONANT CHANNELS)
 744
 745                   Int_t gnephw=pN->GetNDaughters();
 746                   if(gnephw>0){
 747                     Int_t posGN1=pN->GetDaughter(0);
 748                     Int_t posGN2=pN->GetDaughter(1);
 749                     if(posGN2<0)posGN2=posGN1;
 750                     if(posGN1>=0){
 751                       for(Int_t jGN=posGN1;jGN<=posGN2;jGN++){
 752                         AliAODMCParticle *pGN=(AliAODMCParticle*)arrayMC->At(jGN);
 753                         Int_t pdgGN=TMath::Abs(pGN->GetPdgCode());
 754                         if(pdgGN==11||pdgGN==13||pdgGN==211||pdgGN==321||pdgGN==2212){
 755                           fArrayDaugh->AddAt(jGN,flastdaugh);
 756                           flastdaugh++;
 757                         }
 758                       }
 759                     }
 760                   }
 761                 }
 762               }
 763             }
 764           }
 765         }
 766       }
 767     }
 768   }
 769   //  printf("particles to be skipped %d\n",flastdaugh++);
 770 }
 771
 772 void AliAnalysisTaskSEmcCorr::SelectAssociatedParticles(TClonesArray *arrayMC){
 773   fLastAss=0;
 774   fArrayAssoc->Reset(0);
 775   Double_t pt,eta;
 776   for(Int_t jAss=0;jAss<arrayMC->GetEntriesFast();jAss++){
 777     AliAODMCParticle *partAss=(AliAODMCParticle*)arrayMC->At(jAss);
 778     Int_t pdgAss=TMath::Abs(partAss->GetPdgCode());
 779     if(pdgAss==11||pdgAss==13||pdgAss==211||pdgAss==321||pdgAss==2212){// check type: keep only e,mu,pi,K+,protons
 780       pt=partAss->Pt();
 781       eta=TMath::Abs(partAss->Eta());
 782       if(pt<fminPtass)continue;
 783       if(eta>fmaxEtaAss)continue;
 784       if(!partAss->IsPhysicalPrimary()){
 785         if(pdgAss!=11)  continue;// keep ele
 786       }
 787       fArrayAssoc->AddAt(jAss,fLastAss);
 788       fLastAss++;
 789       if(fLastAss==fArrayAssoc->GetSize()-1){
 790         fArrayAssoc->Set(fArrayAssoc->GetSize()+200);
 791       }
 792     }
 793   }
 794   //  printf("n ass part: %d\n",fLastAss++);
 795 }
 796
 797 void AliAnalysisTaskSEmcCorr::FillCorrelationPlots(AliAODMCParticle *part,TClonesArray *arrayMC){
 798
 799
 800   Bool_t hasToSkip=kFALSE;
 801   Int_t index=0,softpi=-1;
 802   Int_t pdgTrig=TMath::Abs(part->GetPdgCode());
 803   if(pdgTrig==421){
 804     pdgTrig=1;
 805   }
 806   else if(pdgTrig==411){
 807     pdgTrig=2;
 808   }
 809   else if(pdgTrig==413){
 810     pdgTrig=3;
 811   }
 812   else if(pdgTrig==11){
 813
 814     pdgTrig=4;
 815   }
 816   else pdgTrig=0;
 817   // check if it is prompt or from B and, for electrons the origin
 818   Int_t labmum=part->GetMother();
 819   if(labmum>=0){
 820     AliAODMCParticle *moth=(AliAODMCParticle*)arrayMC->At(labmum);
 821     Int_t pdgmum=TMath::Abs(moth->GetPdgCode());
 822     if(pdgTrig==4){// for electrons go to the parent hadron
 823       //      labmum=part->GetMother();
 824       //      moth=(AliAODMCParticle*)arrayMC->At(labmum);
 825       //      pdgmum=TMath::Abs(moth->GetPdgCode());
 826       if(!(pdgmum==5||pdgmum==4||(400<pdgmum&&pdgmum<600)||(4000<pdgmum&&pdgmum<6000))){// NON HF electron
 827         if(pdgmum==22){// GAMMA CONV
 828           pdgTrig=5;
 829         }
 830         else if(pdgmum==111||pdgmum==113||pdgmum==221||pdgmum==223){ //DALITZ DECAY
 831           pdgTrig=6;
 832         }
 833         else pdgTrig=7;
 834       }
 835     }
 836     if(pdgmum==413){ // D from Dstar: check soft pion  and go to the grandmother
 837       if(pdgTrig==1){
 838         for(Int_t isp=moth->GetDaughter(0);isp<=moth->GetDaughter(1);isp++){
 839           AliAODMCParticle *sfp=(AliAODMCParticle*)arrayMC->At(isp);
 840           Int_t pdgsp=TMath::Abs(sfp->GetPdgCode());
 841           if(pdgsp==211)softpi=isp;
 842         }
 843       }
 844       labmum=moth->GetMother();
 845       if(labmum>=0){
 846         moth=(AliAODMCParticle*)arrayMC->At(labmum);
 847         pdgmum=TMath::Abs(moth->GetPdgCode());
 848       }
 849       else pdgmum=-1;
 850     }
 851     else if(pdgmum==423){// D*0 -> go to the mother
 852       labmum=moth->GetMother();
 853       if(labmum>=0){
 854         moth=(AliAODMCParticle*)arrayMC->At(labmum);
 855         pdgmum=TMath::Abs(moth->GetPdgCode());
 856       }
 857       else pdgmum=-1;
 858     }
 859     if(pdgmum==5||(500<pdgmum&&pdgmum<600)||(5000<pdgmum&&pdgmum<6000)){
 860       pdgTrig*=-1;
 861     }
 862     else if(pdgmum==4){
 863       labmum=moth->GetMother();
 864       if(labmum>=0){
 865         moth=(AliAODMCParticle*)arrayMC->At(labmum);
 866         pdgmum=TMath::Abs(moth->GetPdgCode());
 867         if(pdgmum==5||(500<pdgmum&&pdgmum<600)||(5000<pdgmum&&pdgmum<6000)){
 868           //      Printf("c quark from b quark");
 869           pdgTrig*=-1;
 870         }
 871       }
 872     }
 873     else {
 874       //      printf("Charm not from charm or beauty \n");
 875     }
 876   }
 877   Double_t phiTrig=part->Phi();
 878   Double_t etaTrig=part->Eta();
 879   Double_t deltaPhi;
 880   Double_t point[8]={ static_cast<Double_t>(pdgTrig),part->Pt(),etaTrig,0,0,0,0,0};//pdg,ptTrig,etaTrig,ptAss,etaAss,deltaPhi,deltaEta,pdgAss
 881   fhMCtrigPart->Fill(point);
 882   for(Int_t jAss=0;jAss<fLastAss;jAss++){
 883     index=fArrayAssoc->At(jAss);
 884     hasToSkip=kFALSE;
 885       for(Int_t jd=0;jd<flastdaugh;jd++){
 886         if(index==fArrayDaugh->At(jd)){
 887           hasToSkip=kTRUE;
 888           break;
 889         }
 890       }
 891       if(hasToSkip)continue;
 892       AliAODMCParticle *partAss=(AliAODMCParticle*)arrayMC->At(index);
 893       point[3]=partAss->Pt();
 894       point[4]=partAss->Eta();
 895       deltaPhi=partAss->Phi()-phiTrig;
 896       if(deltaPhi<-0.5*TMath::Pi())deltaPhi+=2.*TMath::Pi();
 897       else if(deltaPhi>1.5*TMath::Pi())deltaPhi-=2.*TMath::Pi();
 898       point[5]=deltaPhi;
 899       point[6]=partAss->Eta()-etaTrig;
 900       UInt_t pdgassmum=TMath::Abs(partAss->GetPdgCode());
 901       if(index==softpi){
 902         point[7]=-1;
 903       }
 904       else if(pdgassmum==11){
 905         if(partAss->IsPhysicalPrimary())point[7]=1;
 906         Int_t elemum=partAss->GetMother();
 907         if(elemum>0){
 908           AliAODMCParticle *partEleMum=(AliAODMCParticle*)arrayMC->At(elemum);
 909           UInt_t pdgelemum=TMath::Abs(partEleMum->GetPdgCode());
 910           if(pdgelemum==111||pdgelemum==113||pdgelemum==221||pdgelemum==223){// DALITZ DECAY
 911             point[7]=6;
 912           }
 913           else if(pdgelemum==22){// GAMMA CONV
 914             point[7]=7;
 915           }
 916           else if(pdgelemum==4||pdgelemum==5||(pdgelemum>400&&pdgelemum<600)||(pdgelemum>4000&&pdgelemum<6000)){// HF e
 917             point[7]=8;
 918           }
 919           else point[7]=9;
 920         }
 921       }
 922       else if(pdgassmum==13){
 923         point[7]=2;
 924       }
 925       else if(pdgassmum==211){
 926         point[7]=3;
 927       }
 928       else if(pdgassmum==321){
 929         point[7]=4;
 930       }
 931       else if(pdgassmum==2212){
 932         point[7]=5;
 933       }
 934       else point[7]=0;
 935
 936       fhMCcorrel->Fill(point);
 937   }
 938
 939 }
 940
 941
 942
 943 //______________________________________________________________
 944 void AliAnalysisTaskSEmcCorr::FillCorrelationPlotsHadrons(AliAODMCParticle *part,TClonesArray *arrayMC){
 945
 946
 947   Bool_t hasToSkip=kFALSE;
 948   Int_t index=0;
 949   Int_t pdgTrig=TMath::Abs(part->GetPdgCode());
 950   if(pdgTrig==11){
 951     pdgTrig=0;
 952   }
 953   else if(pdgTrig==13){
 954     pdgTrig=1;
 955   }
 956   else if(pdgTrig==211){
 957     pdgTrig=2;
 958   }
 959   else if(pdgTrig==321){
 960     pdgTrig=3;
 961   }
 962   else if(pdgTrig==2212){
 963     pdgTrig=4;
 964   }
 965   else pdgTrig=-1;
 966
 967
 968   Double_t phiTrig=part->Phi();
 969   Double_t etaTrig=part->Eta();
 970   Double_t deltaPhi;
 971   Double_t point[8]={ static_cast<Double_t>(pdgTrig),part->Pt(),etaTrig,0,0,0,0,0};//pdg,ptTrig,etaTrig,ptAss,etaAss,deltaPhi,deltaEta,pdgAss
 972   fhMChadrontrigPart->Fill(point);
 973   for(Int_t jAss=0;jAss<fLastAss;jAss++){
 974     index=fArrayAssoc->At(jAss);
 975     hasToSkip=kFALSE;
 976     for(Int_t jd=0;jd<flastdaugh;jd++){
 977       if(index==fArrayDaugh->At(jd)){
 978         hasToSkip=kTRUE;
 979         break;
 980       }
 981     }
 982     if(hasToSkip)continue;
 983     AliAODMCParticle *partAss=(AliAODMCParticle*)arrayMC->At(index);
 984     point[3]=partAss->Pt();
 985     if(point[3]>point[1])continue;
 986     point[4]=partAss->Eta();
 987     deltaPhi=partAss->Phi()-phiTrig;
 988     if(deltaPhi<-0.5*TMath::Pi())deltaPhi+=2.*TMath::Pi();
 989     else if(deltaPhi>1.5*TMath::Pi())deltaPhi-=2.*TMath::Pi();
 990     point[5]=deltaPhi;
 991     point[6]=partAss->Eta()-etaTrig;
 992     point[7]=TMath::Abs(partAss->GetPdgCode());
 993     if(point[7]==11){
 994       point[7]=0;
 995     }
 996     else if(point[7]==13){
 997       point[7]=1;
 998     }
 999     else if(point[7]==211){
1000       point[7]=2;
1001     }
1002     else if(point[7]==321){
1003       point[7]=3;
1004     }
1005     else if(point[7]==2212){
1006       point[7]=4;
1007     }
1008     else point[7]=-1;
1009
1010     fhMChadroncorrel->Fill(point);
1011   }
1012
1013 }
1014
1015 //_______________________________________________________________
1016 void AliAnalysisTaskSEmcCorr::Terminate(const Option_t*){
1017   //TERMINATE METHOD: NOTHING TO DO
1018
1019
1020
1021 }
1022