[u/mrichter/AliRoot.git] / PWGHF / correlationHF / AliAnalysisTaskSEmcCorr.cxx

/**************************************************************************
 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/////////////////////////////////////////////////////////////
//
// Class AliAnalysisTaskSEmcCorr
// AliAnalysisTaskSE for studying HF-(hadron,electrons) and hadron-hadron correlations
//     at MC level
// Author: Andrea Rossi, andrea.rossi@cern.ch
/////////////////////////////////////////////////////////////


#include <TH1F.h>
#include <TH2F.h>
#include <TAxis.h>
#include <THnSparse.h>
#include <TDatabasePDG.h>
#include <TMath.h>
#include <TROOT.h>
#include <TArrayI.h>
#include "AliAODEvent.h"
#include "AliAODRecoDecayHF2Prong.h"
#include "AliAODRecoDecayHF.h"
#include "AliAODRecoDecay.h"
#include "AliAnalysisDataSlot.h"
#include "AliAnalysisDataContainer.h"
#include "AliAODTrack.h"
#include "AliAODHandler.h"
#include "AliESDtrack.h"
#include "AliAODVertex.h"
#include "AliESDVertex.h"
#include "AliVertexerTracks.h"
#include "AliAODMCParticle.h"
#include "AliAODPid.h"
#include "AliTPCPIDResponse.h"
#include "AliAODMCHeader.h"
#include "AliAnalysisVertexingHF.h"
#include "AliRDHFCutsD0toKpi.h"
#include "AliAODInputHandler.h"
#include "AliAnalysisManager.h"
#include "AliNormalizationCounter.h"
//#include "/Users/administrator/ALICE/CHARM/HFCJ/TestsForProduction/AliAnalysisTaskSEmcCorr.h"
#include "AliAnalysisTaskSEmcCorr.h"

class TCanvas;
class TTree;
class TChain;
class AliAnalysisTaskSE;


ClassImp(AliAnalysisTaskSEmcCorr)

AliAnalysisTaskSEmcCorr::AliAnalysisTaskSEmcCorr() 
: AliAnalysisTaskSE(),
  fReadMC(kTRUE),
  fOnlyKine(kTRUE),
  fDoHadronHadron(kFALSE),
  fNentries(0),
  fhNprongsD0(0),
  fhNprongsD0chargedOnly(0),
  fhNprongsD0chargedRef(0),
  fYrangeTrig(0.5),
  fEtarangeEleTrig(0.8),
  fminDpt(1.),
  fArrayDaugh(0x0),
  flastdaugh(0),
  fArrayAssoc(0x0),
  fLastAss(0),
  fminPtass(0.3),
  fmaxEtaAss(99.),
  fhMCcorrel(0x0),
  fhMCtrigPart(0x0),
  fhMChadroncorrel(0x0),
  fhMChadrontrigPart(0x0)
{// standard constructor
  
}


//________________________________________________________________________
AliAnalysisTaskSEmcCorr::AliAnalysisTaskSEmcCorr(const char *name) 
  : AliAnalysisTaskSE(name),
    fReadMC(kTRUE),
    fOnlyKine(kTRUE),
    fDoHadronHadron(kFALSE),
    fNentries(0),
    fhNprongsD0(0),
    fhNprongsD0chargedOnly(0),
    fhNprongsD0chargedRef(0),
    fYrangeTrig(0.5),
    fEtarangeEleTrig(0.8),
    fminDpt(1.),
    fArrayDaugh(0x0),
    flastdaugh(0),
    fArrayAssoc(0x0),
    fLastAss(0),
    fminPtass(0.3),
    fmaxEtaAss(99.),
    fhMCcorrel(0x0),
    fhMCtrigPart(0x0),
    fhMChadroncorrel(0x0),
    fhMChadrontrigPart(0x0)
{ // default constructor
  
  DefineOutput(1, TH1F::Class());
  DefineOutput(2, TH1F::Class());
  DefineOutput(3, TH1F::Class());
  DefineOutput(4, TH1F::Class());
  DefineOutput(5, THnSparseF::Class());
  DefineOutput(6, THnSparseF::Class());
  DefineOutput(7, THnSparseF::Class());
  DefineOutput(8, THnSparseF::Class());

}

//________________________________________________________________________
AliAnalysisTaskSEmcCorr::~AliAnalysisTaskSEmcCorr(){
  // destructor
  delete fNentries;
  delete fhNprongsD0;
  delete fhNprongsD0chargedOnly;
  delete fhNprongsD0chargedRef;
  delete fhMCcorrel;
  delete fhMCtrigPart;
  delete fhMChadroncorrel;
  delete fhMChadrontrigPart;
  delete fArrayDaugh;
  delete fArrayAssoc;
}

//________________________________________________________________________
void AliAnalysisTaskSEmcCorr::Init()
{
  // Initialization
  fArrayDaugh=new TArrayI(20);

  fArrayAssoc=new TArrayI(200);

}


//________________________________________________________________________
void AliAnalysisTaskSEmcCorr::UserCreateOutputObjects(){
 
  if(!fArrayDaugh)  fArrayDaugh=new TArrayI(20);

  if(!fArrayAssoc)fArrayAssoc=new TArrayI(200);

  fNentries=new TH1F("nentriesChFr", "Analyzed sample properties", 21,0.5,21.5);
  fNentries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
  fNentries->GetXaxis()->SetBinLabel(2,"nEvSel");
  fNentries->GetXaxis()->SetBinLabel(3,"nEvPile-up Rej");
  fNentries->GetXaxis()->SetBinLabel(4,"nEvGoodVtxS");
  fNentries->GetXaxis()->SetBinLabel(5,"nEvRejVtxZ");
  fNentries->GetXaxis()->SetBinLabel(6,"nD0");

  fhNprongsD0=new TH1F("fhNprongsD0","fhNprongsD0",20,-0.5,19.5);
  fhNprongsD0chargedOnly=new TH1F("fhNprongsD0chargedOnly","fhNprongsD0chargedOnly",20,-0.5,19.5);
  fhNprongsD0chargedRef=new TH1F("fhNprongsD0chargedRef","fhNprongsD0chargedRef",20,-0.5,19.5);


  Int_t nbinsCorrMC[8]={15,20,20,20,50,63,20,11};
  Double_t binlowCorrMC[8]={-7.5,0.,-1.,0.,-5.,-0.5*TMath::Pi(),-2,-1.5};
  Double_t binupCorrMC[8]={7.5,20.,1.,2.,5.,1.5*TMath::Pi(),2.,9.5};
  fhMCcorrel=new THnSparseF("fhMCcorrel","fhMCcorrel;pdg;ptTrig;etaTrig;ptAss;etaAss;deltaPhi;deltaEta;pdgAss;",8,nbinsCorrMC,binlowCorrMC,binupCorrMC);

  Int_t nbinsTrigMC[3]={15,20,20};
  Double_t binlowTrigMC[3]={-7.5,0.,-1.};
  Double_t binupTrigMC[3]={7.5,20.,1.};
  fhMCtrigPart=new THnSparseF("fhMCtrigPart","fhMCcorrel;pdg;ptTrig;etaTrig;",3,nbinsTrigMC,binlowTrigMC,binupTrigMC);


  Int_t nbinsCorrMChadron[8]={6,20,20,20,50,63,20,7};
  Double_t binlowCorrMChadron[8]={-1.5,0.,-1.,0.,-5.,-0.5*TMath::Pi(),-2,-1.5};
  Double_t binupCorrMChadron[8]={4.5,20.,1.,2.,5.,1.5*TMath::Pi(),2.,5.5};
  fhMChadroncorrel=new THnSparseF("fhMChadroncorrel","fhMChadroncorrel;pdg;ptTrig;etaTrig;ptAss;etaAss;deltaPhi;deltaEta;pdgAss;",8,nbinsCorrMChadron,binlowCorrMChadron,binupCorrMChadron);

  Int_t nbinsTrigMChadron[3]={6,20,20};
  Double_t binlowTrigMChadron[3]={-1.5,0.,-1.};
  Double_t binupTrigMChadron[3]={4.5,20.,1.};
  fhMChadrontrigPart=new THnSparseF("fhMChadrontrigPart","fhMChadroncorrel;pdg;ptTrig;etaTrig;",3,nbinsTrigMChadron,binlowTrigMChadron,binupTrigMChadron);
  

  PostData(1,fNentries);
  PostData(2,fhNprongsD0);
  PostData(3,fhNprongsD0chargedOnly);
  PostData(4,fhNprongsD0chargedRef);
  PostData(5,fhMCcorrel);
  PostData(6,fhMCtrigPart);
  PostData(7,fhMChadroncorrel);
  PostData(8,fhMChadrontrigPart);
}


//________________________________________________________________________
void AliAnalysisTaskSEmcCorr::UserExec(Option_t */*option*/){


  // Execute analysis for current event:
  // heavy flavor candidates association to MC truth
  
  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
  if (!aod) {
    Printf("ERROR: aod not available");
    return;
  }


  fNentries->Fill(1);

  // AOD primary vertex
  
  
  if(!fOnlyKine){
    AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
    TString primTitle = vtx1->GetTitle();
    if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) { 
      
      fNentries->Fill(3);
      
    }
    else {
      PostData(1,fNentries);
      return;
      
    }
  }


  TClonesArray *arrayMC=0x0;
  AliAODMCHeader *aodmcHeader=0x0;
  Double_t vtxTrue[3];
 
  
  if(fReadMC){
    //    aod->GetList()->Print();
    // load MC particles
    arrayMC = 
      (TClonesArray*)aod->GetList()->At(23);//FindObject("mcparticles");
    if(!arrayMC) {
      Printf("AliAnalysisTaskSEmcCorr::UserExec: MC particles branch not found!\n");
      fNentries->Fill(2);
      PostData(1,fNentries);
      return;
    }
    // load MC header
    aodmcHeader = 
      (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
    if(!aodmcHeader) {
      Printf("AliAnalysisTaskSEmcCorr::UserExec: MC header branch not found!\n");
            fNentries->Fill(3);
      PostData(1,fNentries);
      return;
    }
    fNentries->Fill(4);
    //    printf("N MC entries: %d \n",arrayMC->GetEntriesFast());
    // MC primary vertex
    aodmcHeader->GetVertex(vtxTrue);

    // FILL HISTOS FOR D0 mesons, c quarks and D0 from B properties
    SelectAssociatedParticles(arrayMC);
    for(Int_t jp=0;jp<arrayMC->GetEntriesFast();jp++){

      AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->At(jp);
      Int_t pdg=TMath::Abs(part->GetPdgCode());

      if(pdg==421||pdg==411||pdg==413){

	if(TMath::Abs(part->Y())>fYrangeTrig)continue;
	if(TMath::Abs(part->Pt())<fminDpt)continue;
	FillSkipParticleArray(part,arrayMC,jp);
	FillCorrelationPlots(part,arrayMC);
	

      }
      else if(pdg==11){

	if(TMath::Abs(part->Eta())>fEtarangeEleTrig)continue;
	if(TMath::Abs(part->Pt())<fminDpt)continue;
	FillSkipParticleArray(part,arrayMC,jp);
	FillCorrelationPlots(part,arrayMC);


      }       

      // NOW STUDY OF N-PRONGS
      if(pdg==421){
	Bool_t semilept=kFALSE;
	Int_t nprongsD0=0;
	Int_t nprongsD0charged=0;

	Int_t lab1=part->GetDaughter(0);
	Int_t lab2=part->GetDaughter(1);
	if(lab1>=0&&lab2>=0){
	  for(Int_t jd=lab1;jd<=lab2;jd++){
	    AliAODMCParticle *d=(AliAODMCParticle*)arrayMC->At(jd);
	    Int_t pdgd=TMath::Abs(d->GetPdgCode());
	    if(pdgd>=11&&pdgd<=16){
	      semilept=kTRUE;	   
	      break;
	    }
	    if(pdgd==211||pdgd==111||pdgd==321||pdgd==311||pdgd==310||pdgd==130){
	      nprongsD0++;
	    }
	    if(pdgd==211||pdgd==321){
	      nprongsD0charged++;
	    }	    
	    else{	    
	      Int_t lab1d=d->GetDaughter(0);
	      Int_t lab2d=d->GetDaughter(1);
	      if(lab1d>=0&&lab2d>=0){
		for(Int_t jdd=lab1d;jdd<=lab2d;jdd++){
		  AliAODMCParticle *dd=(AliAODMCParticle*)arrayMC->At(jdd);
		  Int_t pdgdd=TMath::Abs(dd->GetPdgCode());
		  if(pdgdd==211||pdgdd==111||pdgdd==321||pdgdd==311||pdgdd==310||pdgdd==130){
		    nprongsD0++;
		  }
		  if(pdgd==211||pdgd==321){
		    nprongsD0charged++;
		  } 
		}
	      }
	    }
	  }
	}         
	if(!semilept){
	  fhNprongsD0->Fill(nprongsD0);
	  fhNprongsD0chargedOnly->Fill(nprongsD0charged);
	  fhNprongsD0chargedRef->Fill(nprongsD0charged);
	}
      }

    }
    
    if(fDoHadronHadron){
      for(Int_t jp=0;jp<fLastAss;jp++){
	
	AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->At(fArrayAssoc->At(jp));
	Int_t pdg=TMath::Abs(part->GetPdgCode());      
	if(pdg==11||pdg==13||pdg==211||pdg==2212||pdg==321){
	  if(TMath::Abs(part->Eta())>fEtarangeEleTrig)continue;
	  if(TMath::Abs(part->Pt())<fminDpt)continue;
	  FillSkipParticleArray(part,arrayMC,fArrayAssoc->At(jp));
	  FillCorrelationPlotsHadrons(part,arrayMC);
	}
      }
      
    }
    
  }
  

  PostData(1,fNentries);
  PostData(2,fhNprongsD0);
  PostData(3,fhNprongsD0chargedOnly);
  PostData(4,fhNprongsD0chargedRef);
  PostData(5,fhMCcorrel);
  PostData(6,fhMCtrigPart);
  PostData(7,fhMChadroncorrel);
  PostData(8,fhMChadrontrigPart);
}


void AliAnalysisTaskSEmcCorr::FillSkipParticleArray(AliAODMCParticle *part,TClonesArray *arrayMC,Int_t partPos){
  fArrayDaugh->Reset(0);
  fArrayDaugh->AddAt(partPos,0);
  flastdaugh=1;
  
  // Loop on daughters and nephews
  if(part->GetNDaughters()>0){
    Int_t pos1=part->GetDaughter(0);
    Int_t pos2=part->GetDaughter(1);
    if(pos2<0)pos2=pos1;
    if(pos1>=0){
      for(Int_t jd=pos1;jd<=pos2;jd++){ 
	AliAODMCParticle *pd=(AliAODMCParticle*)arrayMC->At(jd);
	Int_t pdgd=TMath::Abs(pd->GetPdgCode());
	if(pdgd==11||pdgd==13||pdgd==211||pdgd==321||pdgd==2212){
	  fArrayDaugh->AddAt(jd,flastdaugh);
	  flastdaugh++;
	}
	else if(pdgd!=12&&pdgd!=14&&pdgd!=16&&pdgd!=310&&pdgd!=111){// CHECK NEPHEWS (FOR RESONANT CHANNELS)
	  Int_t nephw=pd->GetNDaughters();
	  if(nephw>0){
	    Int_t posN1=pd->GetDaughter(0);
	    Int_t posN2=pd->GetDaughter(1);
	    if(posN2<0)posN2=posN1;
	    if(posN1>=0){
	      for(Int_t jN=posN1;jN<=posN2;jN++){ 
		AliAODMCParticle *pN=(AliAODMCParticle*)arrayMC->At(jN);
		Int_t pdgN=TMath::Abs(pN->GetPdgCode());
		if(pdgN==11||pdgN==13||pdgN==211||pdgN==321||pdgN==2212){
		  fArrayDaugh->AddAt(jN,flastdaugh);
		  flastdaugh++;
		}
		else if(pdgN!=12&&pdgN!=14&&pdgN!=16&&pdgN!=310&&pdgN!=111){// CHECK one step more (FOR RESONANT CHANNELS)	 
		  
		  Int_t gnephw=pN->GetNDaughters();
		  if(gnephw>0){
		    Int_t posGN1=pN->GetDaughter(0);
		    Int_t posGN2=pN->GetDaughter(1);
		    if(posGN2<0)posGN2=posGN1;
		    if(posGN1>=0){
		      for(Int_t jGN=posGN1;jGN<=posGN2;jGN++){ 
			AliAODMCParticle *pGN=(AliAODMCParticle*)arrayMC->At(jGN);
			Int_t pdgGN=TMath::Abs(pGN->GetPdgCode());
			if(pdgGN==11||pdgGN==13||pdgGN==211||pdgGN==321||pdgGN==2212){
			  fArrayDaugh->AddAt(jGN,flastdaugh);
			  flastdaugh++;			 
			}
		      }
		    }
		  }
		}
	      }
	    }
	  }
	}
      }
    } 
  }
  //  printf("particles to be skipped %d\n",flastdaugh++);
}

void AliAnalysisTaskSEmcCorr::SelectAssociatedParticles(TClonesArray *arrayMC){
  fLastAss=0;
  fArrayAssoc->Reset(0);
  Double_t pt,eta;
  for(Int_t jAss=0;jAss<arrayMC->GetEntriesFast();jAss++){
    AliAODMCParticle *partAss=(AliAODMCParticle*)arrayMC->At(jAss);
    Int_t pdgAss=TMath::Abs(partAss->GetPdgCode());
    if(pdgAss==11||pdgAss==13||pdgAss==211||pdgAss==321||pdgAss==2212){// check type: keep only e,mu,pi,K+,protons
      pt=partAss->Pt();
      eta=TMath::Abs(partAss->Eta());
      if(pt<fminPtass)continue;
      if(eta>fmaxEtaAss)continue;
      if(!partAss->IsPhysicalPrimary()){
	if(pdgAss!=11)	continue;// keep ele
      }
      fArrayAssoc->AddAt(jAss,fLastAss);
      fLastAss++;
    }
  }  
  //  printf("n ass part: %d\n",fLastAss++);
}

void AliAnalysisTaskSEmcCorr::FillCorrelationPlots(AliAODMCParticle *part,TClonesArray *arrayMC){
  

  Bool_t hasToSkip=kFALSE;  
  Int_t index=0,softpi=-1;
  Int_t pdgTrig=TMath::Abs(part->GetPdgCode());
  if(pdgTrig==421){
    pdgTrig=1;
  }
  else if(pdgTrig==411){
    pdgTrig=2;
  }
  else if(pdgTrig==413){
    pdgTrig=3;
  }
  else if(pdgTrig==11){
    
    pdgTrig=4;
  }
  else pdgTrig=0;
  // check if it is prompt or from B and, for electrons the origin
  Int_t labmum=part->GetMother();
  if(labmum>=0){
    AliAODMCParticle *moth=(AliAODMCParticle*)arrayMC->At(labmum);
    Int_t pdgmum=TMath::Abs(moth->GetPdgCode());
    if(pdgTrig==4){// for electrons go to the parent hadron
      //      labmum=part->GetMother();
      //      moth=(AliAODMCParticle*)arrayMC->At(labmum);
      //      pdgmum=TMath::Abs(moth->GetPdgCode());
      if(!(pdgmum==5||pdgmum==4||(400<pdgmum&&pdgmum<600)||(4000<pdgmum&&pdgmum<6000))){// NON HF electron
	if(pdgmum==22){// GAMMA CONV
	  pdgTrig=5;
	}
	else if(pdgmum==111||pdgmum==113||pdgmum==221||pdgmum==223){ //DALITZ DECAY
	  pdgTrig=6;
	}
	else pdgTrig=7;
      }
    }
    if(pdgmum==413){ // D from Dstar: check soft pion  and go to the grandmother
      if(pdgTrig==1){
	for(Int_t isp=moth->GetDaughter(0);isp<=moth->GetDaughter(1);isp++){
	  AliAODMCParticle *sfp=(AliAODMCParticle*)arrayMC->At(isp);
	  Int_t pdgsp=TMath::Abs(sfp->GetPdgCode());
	  if(pdgsp==211)softpi=isp;
	}
      }
      labmum=moth->GetMother();
      moth=(AliAODMCParticle*)arrayMC->At(labmum);
      pdgmum=TMath::Abs(moth->GetPdgCode());
    }
    else if(pdgmum==423){// D*0 -> go to the mother
      labmum=moth->GetMother();
      moth=(AliAODMCParticle*)arrayMC->At(labmum);
      pdgmum=TMath::Abs(moth->GetPdgCode());
    }
    if(pdgmum==5||(500<pdgmum&&pdgmum<600)||(5000<pdgmum&&pdgmum<6000)){
      pdgTrig*=-1;
    }
    else if(pdgmum==4){
      labmum=moth->GetMother();
      if(labmum>=0){
	moth=(AliAODMCParticle*)arrayMC->At(labmum);
	pdgmum=TMath::Abs(moth->GetPdgCode());
	if(pdgmum==5||(500<pdgmum&&pdgmum<600)||(5000<pdgmum&&pdgmum<6000)){
	  //	  Printf("c quark from b quark");
	  pdgTrig*=-1;
	}   
      }
    }
    else {
      //      printf("Charm not from charm or beauty \n");
    }
  }
  Double_t phiTrig=part->Phi();
  Double_t etaTrig=part->Eta();
  Double_t deltaPhi;
  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
  fhMCtrigPart->Fill(point);
  for(Int_t jAss=0;jAss<fLastAss;jAss++){
    index=fArrayAssoc->At(jAss);
    hasToSkip=kFALSE;  
      for(Int_t jd=0;jd<flastdaugh;jd++){
	if(index==fArrayDaugh->At(jd)){
	  hasToSkip=kTRUE;  
	  break;
	}
      }
      if(hasToSkip)continue;
      AliAODMCParticle *partAss=(AliAODMCParticle*)arrayMC->At(index);
      point[3]=partAss->Pt();
      point[4]=partAss->Eta();
      deltaPhi=partAss->Phi()-phiTrig;
      if(deltaPhi<-0.5*TMath::Pi())deltaPhi+=2.*TMath::Pi();
      else if(deltaPhi>1.5*TMath::Pi())deltaPhi-=2.*TMath::Pi();
      point[5]=deltaPhi;
      point[6]=partAss->Eta()-etaTrig;      
      UInt_t pdgassmum=TMath::Abs(partAss->GetPdgCode());
      if(index==softpi){
	point[7]=-1;
      }
      else if(pdgassmum==11){
	if(partAss->IsPhysicalPrimary())point[7]=1;
	Int_t elemum=partAss->GetMother();
	if(elemum>0){
	  AliAODMCParticle *partEleMum=(AliAODMCParticle*)arrayMC->At(elemum);
	  UInt_t pdgelemum=TMath::Abs(partEleMum->GetPdgCode());
	  if(pdgelemum==111||pdgelemum==113||pdgelemum==221||pdgelemum==223){// DALITZ DECAY
	    point[7]=6;
	  }
	  else if(pdgelemum==22){// GAMMA CONV
	    point[7]=7;
	  }
	  else if(pdgelemum==4||pdgelemum==5||(pdgelemum>400&&pdgelemum<600)||(pdgelemum>4000&&pdgelemum<6000)){// HF e
	    point[7]=8;
	  }
	  else point[7]=9;
	}
      }
      else if(pdgassmum==13){
	point[7]=2;
      }
      else if(pdgassmum==211){
	point[7]=3;
      }
      else if(pdgassmum==321){
	point[7]=4;
      }
      else if(pdgassmum==2212){
	point[7]=5;
      }
      else point[7]=0;

      fhMCcorrel->Fill(point);
  }

}


//______________________________________________________________
void AliAnalysisTaskSEmcCorr::FillCorrelationPlotsHadrons(AliAODMCParticle *part,TClonesArray *arrayMC){
  

  Bool_t hasToSkip=kFALSE;  
  Int_t index=0;
  Int_t pdgTrig=TMath::Abs(part->GetPdgCode());
  if(pdgTrig==11){
    pdgTrig=0;
  }
  else if(pdgTrig==13){
    pdgTrig=1;
  }
  else if(pdgTrig==211){
    pdgTrig=2;
  }
  else if(pdgTrig==321){
    pdgTrig=3;
  }
  else if(pdgTrig==2212){
    pdgTrig=4;
  }
  else pdgTrig=-1;
  
  
  Double_t phiTrig=part->Phi();
  Double_t etaTrig=part->Eta();
  Double_t deltaPhi;
  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
  fhMChadrontrigPart->Fill(point);
  for(Int_t jAss=0;jAss<fLastAss;jAss++){
    index=fArrayAssoc->At(jAss);
    hasToSkip=kFALSE;  
    for(Int_t jd=0;jd<flastdaugh;jd++){
      if(index==fArrayDaugh->At(jd)){
	hasToSkip=kTRUE;  
	break;
      }
    }
    if(hasToSkip)continue;
    AliAODMCParticle *partAss=(AliAODMCParticle*)arrayMC->At(index);
    point[3]=partAss->Pt();
    if(point[3]>point[1])continue;
    point[4]=partAss->Eta();
    deltaPhi=partAss->Phi()-phiTrig;
    if(deltaPhi<-0.5*TMath::Pi())deltaPhi+=2.*TMath::Pi();
    else if(deltaPhi>1.5*TMath::Pi())deltaPhi-=2.*TMath::Pi();
    point[5]=deltaPhi;
    point[6]=partAss->Eta()-etaTrig;
    point[7]=TMath::Abs(partAss->GetPdgCode());
    if(point[7]==11){
      point[7]=0;
    }
    else if(point[7]==13){
      point[7]=1;
    }
    else if(point[7]==211){
      point[7]=2;
    }
    else if(point[7]==321){
      point[7]=3;
    }
    else if(point[7]==2212){
      point[7]=4;
    }
    else point[7]=-1;
    
    fhMChadroncorrel->Fill(point);
  }
 
}

//_______________________________________________________________
void AliAnalysisTaskSEmcCorr::Terminate(const Option_t*){
  //TERMINATE METHOD: NOTHING TO DO


}
Commit	Line	Data
d33db6d9 AR	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	//
bb236562	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
d33db6d9 AR	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 "AliAODEvent.h"
	34	#include "AliAODRecoDecayHF2Prong.h"
	35	#include "AliAODRecoDecayHF.h"
	36	#include "AliAODRecoDecay.h"
	37	#include "AliAnalysisDataSlot.h"
	38	#include "AliAnalysisDataContainer.h"
	39	#include "AliAODTrack.h"
	40	#include "AliAODHandler.h"
	41	#include "AliESDtrack.h"
	42	#include "AliAODVertex.h"
	43	#include "AliESDVertex.h"
	44	#include "AliVertexerTracks.h"
	45	#include "AliAODMCParticle.h"
	46	#include "AliAODPid.h"
	47	#include "AliTPCPIDResponse.h"
	48	#include "AliAODMCHeader.h"
	49	#include "AliAnalysisVertexingHF.h"
	50	#include "AliRDHFCutsD0toKpi.h"
	51	#include "AliAODInputHandler.h"
	52	#include "AliAnalysisManager.h"
	53	#include "AliNormalizationCounter.h"
	54	//#include "/Users/administrator/ALICE/CHARM/HFCJ/TestsForProduction/AliAnalysisTaskSEmcCorr.h"
	55	#include "AliAnalysisTaskSEmcCorr.h"
	56
	57	class TCanvas;
	58	class TTree;
	59	class TChain;
	60	class AliAnalysisTaskSE;
	61
	62
	63	ClassImp(AliAnalysisTaskSEmcCorr)
	64
	65	AliAnalysisTaskSEmcCorr::AliAnalysisTaskSEmcCorr()
	66	: AliAnalysisTaskSE(),
	67	fReadMC(kTRUE),
	68	fOnlyKine(kTRUE),
	69	fDoHadronHadron(kFALSE),
	70	fNentries(0),
	71	fhNprongsD0(0),
	72	fhNprongsD0chargedOnly(0),
	73	fhNprongsD0chargedRef(0),
	74	fYrangeTrig(0.5),
	75	fEtarangeEleTrig(0.8),
	76	fminDpt(1.),
	77	fArrayDaugh(0x0),
	78	flastdaugh(0),
	79	fArrayAssoc(0x0),
	80	fLastAss(0),
	81	fminPtass(0.3),
	82	fmaxEtaAss(99.),
	83	fhMCcorrel(0x0),
	84	fhMCtrigPart(0x0),
	85	fhMChadroncorrel(0x0),
86	fhMChadrontrigPart(0x0)
87	{// standard constructor
88
89	}
90
91
92
93	//________________________________________________________________________
94	AliAnalysisTaskSEmcCorr::AliAnalysisTaskSEmcCorr(const char *name)
95	: AliAnalysisTaskSE(name),
96	fReadMC(kTRUE),
97	fOnlyKine(kTRUE),
98	fDoHadronHadron(kFALSE),
99	fNentries(0),
100	fhNprongsD0(0),
101	fhNprongsD0chargedOnly(0),
102	fhNprongsD0chargedRef(0),
103	fYrangeTrig(0.5),
104	fEtarangeEleTrig(0.8),
105	fminDpt(1.),
106	fArrayDaugh(0x0),
107	flastdaugh(0),
108	fArrayAssoc(0x0),
109	fLastAss(0),
110	fminPtass(0.3),
111	fmaxEtaAss(99.),
112	fhMCcorrel(0x0),
113	fhMCtrigPart(0x0),
114	fhMChadroncorrel(0x0),
115	fhMChadrontrigPart(0x0)
116	{ // default constructor
117
118	DefineOutput(1, TH1F::Class());
119	DefineOutput(2, TH1F::Class());
120	DefineOutput(3, TH1F::Class());
121	DefineOutput(4, TH1F::Class());
122	DefineOutput(5, THnSparseF::Class());
123	DefineOutput(6, THnSparseF::Class());
124	DefineOutput(7, THnSparseF::Class());
125	DefineOutput(8, THnSparseF::Class());
126
127	}
128
129	//________________________________________________________________________
130	AliAnalysisTaskSEmcCorr::~AliAnalysisTaskSEmcCorr(){
131	// destructor
132	delete fNentries;
133	delete fhNprongsD0;
134	delete fhNprongsD0chargedOnly;
135	delete fhNprongsD0chargedRef;
136	delete fhMCcorrel;
137	delete fhMCtrigPart;
138	delete fhMChadroncorrel;
139	delete fhMChadrontrigPart;
140	delete fArrayDaugh;
141	delete fArrayAssoc;
142	}
143
144	//________________________________________________________________________
145	void AliAnalysisTaskSEmcCorr::Init()
146	{
147	// Initialization
148	fArrayDaugh=new TArrayI(20);
149
150	fArrayAssoc=new TArrayI(200);
151
152	}
153
154
155	//________________________________________________________________________
156	void AliAnalysisTaskSEmcCorr::UserCreateOutputObjects(){
157
158	if(!fArrayDaugh) fArrayDaugh=new TArrayI(20);
159
160	if(!fArrayAssoc)fArrayAssoc=new TArrayI(200);
161
162	fNentries=new TH1F("nentriesChFr", "Analyzed sample properties", 21,0.5,21.5);
163	fNentries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
164	fNentries->GetXaxis()->SetBinLabel(2,"nEvSel");
165	fNentries->GetXaxis()->SetBinLabel(3,"nEvPile-up Rej");
166	fNentries->GetXaxis()->SetBinLabel(4,"nEvGoodVtxS");
167	fNentries->GetXaxis()->SetBinLabel(5,"nEvRejVtxZ");
168	fNentries->GetXaxis()->SetBinLabel(6,"nD0");
169
170	fhNprongsD0=new TH1F("fhNprongsD0","fhNprongsD0",20,-0.5,19.5);
171	fhNprongsD0chargedOnly=new TH1F("fhNprongsD0chargedOnly","fhNprongsD0chargedOnly",20,-0.5,19.5);
172	fhNprongsD0chargedRef=new TH1F("fhNprongsD0chargedRef","fhNprongsD0chargedRef",20,-0.5,19.5);
173
174
2942f542	175	Int_t nbinsCorrMC[8]={15,20,20,20,50,63,20,11};
d33db6d9 AR	176	Double_t binlowCorrMC[8]={-7.5,0.,-1.,0.,-5.,-0.5*TMath::Pi(),-2,-1.5};
	177	Double_t binupCorrMC[8]={7.5,20.,1.,2.,5.,1.5*TMath::Pi(),2.,9.5};
	178	fhMCcorrel=new THnSparseF("fhMCcorrel","fhMCcorrel;pdg;ptTrig;etaTrig;ptAss;etaAss;deltaPhi;deltaEta;pdgAss;",8,nbinsCorrMC,binlowCorrMC,binupCorrMC);
	179
2942f542	180	Int_t nbinsTrigMC[3]={15,20,20};
d33db6d9 AR	181	Double_t binlowTrigMC[3]={-7.5,0.,-1.};
	182	Double_t binupTrigMC[3]={7.5,20.,1.};
	183	fhMCtrigPart=new THnSparseF("fhMCtrigPart","fhMCcorrel;pdg;ptTrig;etaTrig;",3,nbinsTrigMC,binlowTrigMC,binupTrigMC);
	184
	185
2942f542	186	Int_t nbinsCorrMChadron[8]={6,20,20,20,50,63,20,7};
d33db6d9 AR	187	Double_t binlowCorrMChadron[8]={-1.5,0.,-1.,0.,-5.,-0.5*TMath::Pi(),-2,-1.5};
	188	Double_t binupCorrMChadron[8]={4.5,20.,1.,2.,5.,1.5*TMath::Pi(),2.,5.5};
	189	fhMChadroncorrel=new THnSparseF("fhMChadroncorrel","fhMChadroncorrel;pdg;ptTrig;etaTrig;ptAss;etaAss;deltaPhi;deltaEta;pdgAss;",8,nbinsCorrMChadron,binlowCorrMChadron,binupCorrMChadron);
	190
2942f542	191	Int_t nbinsTrigMChadron[3]={6,20,20};
d33db6d9 AR	192	Double_t binlowTrigMChadron[3]={-1.5,0.,-1.};
	193	Double_t binupTrigMChadron[3]={4.5,20.,1.};
	194	fhMChadrontrigPart=new THnSparseF("fhMChadrontrigPart","fhMChadroncorrel;pdg;ptTrig;etaTrig;",3,nbinsTrigMChadron,binlowTrigMChadron,binupTrigMChadron);
	195
	196
	197	PostData(1,fNentries);
	198	PostData(2,fhNprongsD0);
	199	PostData(3,fhNprongsD0chargedOnly);
	200	PostData(4,fhNprongsD0chargedRef);
	201	PostData(5,fhMCcorrel);
	202	PostData(6,fhMCtrigPart);
	203	PostData(7,fhMChadroncorrel);
	204	PostData(8,fhMChadrontrigPart);
	205	}
	206
	207
	208
	209	//________________________________________________________________________
	210	void AliAnalysisTaskSEmcCorr::UserExec(Option_t /option*/){
	211
	212
	213
	214	// Execute analysis for current event:
	215	// heavy flavor candidates association to MC truth
	216
	217	AliAODEvent aod = dynamic_cast<AliAODEvent> (InputEvent());
	218	if (!aod) {
	219	Printf("ERROR: aod not available");
	220	return;
	221	}
	222
	223
d33db6d9 AR	224
	225	fNentries->Fill(1);
	226
	227	// AOD primary vertex
	228
	229
	230	if(!fOnlyKine){
	231	AliAODVertex vtx1 = (AliAODVertex)aod->GetPrimaryVertex();
	232	TString primTitle = vtx1->GetTitle();
	233	if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
	234
	235	fNentries->Fill(3);
	236
	237	}
	238	else {
	239	PostData(1,fNentries);
	240	return;
	241
	242	}
	243	}
	244
	245
	246	TClonesArray *arrayMC=0x0;
	247	AliAODMCHeader *aodmcHeader=0x0;
	248	Double_t vtxTrue[3];
	249
	250
	251	if(fReadMC){
	252	// aod->GetList()->Print();
	253	// load MC particles
	254	arrayMC =
	255	(TClonesArray*)aod->GetList()->At(23);//FindObject("mcparticles");
	256	if(!arrayMC) {
	257	Printf("AliAnalysisTaskSEmcCorr::UserExec: MC particles branch not found!\n");
	258	fNentries->Fill(2);
	259	PostData(1,fNentries);
	260	return;
	261	}
	262	// load MC header
	263	aodmcHeader =
	264	(AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
	265	if(!aodmcHeader) {
	266	Printf("AliAnalysisTaskSEmcCorr::UserExec: MC header branch not found!\n");
	267	fNentries->Fill(3);
	268	PostData(1,fNentries);
	269	return;
	270	}
	271	fNentries->Fill(4);
	272	// printf("N MC entries: %d \n",arrayMC->GetEntriesFast());
	273	// MC primary vertex
	274	aodmcHeader->GetVertex(vtxTrue);
	275
	276	// FILL HISTOS FOR D0 mesons, c quarks and D0 from B properties
	277	SelectAssociatedParticles(arrayMC);
	278	for(Int_t jp=0;jp<arrayMC->GetEntriesFast();jp++){
	279
	280	AliAODMCParticle part=(AliAODMCParticle)arrayMC->At(jp);
	281	Int_t pdg=TMath::Abs(part->GetPdgCode());
	282
	283	if(pdg==421\|\|pdg==411\|\|pdg==413){
	284
	285	if(TMath::Abs(part->Y())>fYrangeTrig)continue;
	286	if(TMath::Abs(part->Pt())<fminDpt)continue;
	287	FillSkipParticleArray(part,arrayMC,jp);
288	FillCorrelationPlots(part,arrayMC);
289
290
291	}
292	else if(pdg==11){
293
294	if(TMath::Abs(part->Eta())>fEtarangeEleTrig)continue;
295	if(TMath::Abs(part->Pt())<fminDpt)continue;
296	FillSkipParticleArray(part,arrayMC,jp);
297	FillCorrelationPlots(part,arrayMC);
298
299
300	}
301
302	// NOW STUDY OF N-PRONGS
303	if(pdg==421){
304	Bool_t semilept=kFALSE;
305	Int_t nprongsD0=0;
306	Int_t nprongsD0charged=0;
307
308	Int_t lab1=part->GetDaughter(0);
309	Int_t lab2=part->GetDaughter(1);
310	if(lab1>=0&&lab2>=0){
311	for(Int_t jd=lab1;jd<=lab2;jd++){
312	AliAODMCParticle d=(AliAODMCParticle)arrayMC->At(jd);
313	Int_t pdgd=TMath::Abs(d->GetPdgCode());
314	if(pdgd>=11&&pdgd<=16){
315	semilept=kTRUE;
316	break;
317	}
318	if(pdgd==211\|\|pdgd==111\|\|pdgd==321\|\|pdgd==311\|\|pdgd==310\|\|pdgd==130){
319	nprongsD0++;
320	}
321	if(pdgd==211\|\|pdgd==321){
322	nprongsD0charged++;
323	}
324	else{
325	Int_t lab1d=d->GetDaughter(0);
326	Int_t lab2d=d->GetDaughter(1);
327	if(lab1d>=0&&lab2d>=0){
328	for(Int_t jdd=lab1d;jdd<=lab2d;jdd++){
329	AliAODMCParticle dd=(AliAODMCParticle)arrayMC->At(jdd);
330	Int_t pdgdd=TMath::Abs(dd->GetPdgCode());
331	if(pdgdd==211\|\|pdgdd==111\|\|pdgdd==321\|\|pdgdd==311\|\|pdgdd==310\|\|pdgdd==130){
332	nprongsD0++;
333	}
334	if(pdgd==211\|\|pdgd==321){
335	nprongsD0charged++;
336	}
337	}
338	}
339	}
340	}
341	}
342	if(!semilept){
343	fhNprongsD0->Fill(nprongsD0);
344	fhNprongsD0chargedOnly->Fill(nprongsD0charged);
345	fhNprongsD0chargedRef->Fill(nprongsD0charged);
346	}
347	}
348
349	}
350
351	if(fDoHadronHadron){
352	for(Int_t jp=0;jp<fLastAss;jp++){
353
354	AliAODMCParticle part=(AliAODMCParticle)arrayMC->At(fArrayAssoc->At(jp));
355	Int_t pdg=TMath::Abs(part->GetPdgCode());
356	if(pdg==11\|\|pdg==13\|\|pdg==211\|\|pdg==2212\|\|pdg==321){
357	if(TMath::Abs(part->Eta())>fEtarangeEleTrig)continue;
358	if(TMath::Abs(part->Pt())<fminDpt)continue;
359	FillSkipParticleArray(part,arrayMC,fArrayAssoc->At(jp));
360	FillCorrelationPlotsHadrons(part,arrayMC);
361	}
362	}
363
364	}
365
366	}
367
368
369	PostData(1,fNentries);
370	PostData(2,fhNprongsD0);
371	PostData(3,fhNprongsD0chargedOnly);
372	PostData(4,fhNprongsD0chargedRef);
373	PostData(5,fhMCcorrel);
374	PostData(6,fhMCtrigPart);
375	PostData(7,fhMChadroncorrel);
376	PostData(8,fhMChadrontrigPart);
377	}
378
379
380
381	void AliAnalysisTaskSEmcCorr::FillSkipParticleArray(AliAODMCParticle part,TClonesArray arrayMC,Int_t partPos){
382	fArrayDaugh->Reset(0);
383	fArrayDaugh->AddAt(partPos,0);
384	flastdaugh=1;
385
386	// Loop on daughters and nephews
387	if(part->GetNDaughters()>0){
388	Int_t pos1=part->GetDaughter(0);
389	Int_t pos2=part->GetDaughter(1);
390	if(pos2<0)pos2=pos1;
391	if(pos1>=0){
392	for(Int_t jd=pos1;jd<=pos2;jd++){
393	AliAODMCParticle pd=(AliAODMCParticle)arrayMC->At(jd);
394	Int_t pdgd=TMath::Abs(pd->GetPdgCode());
395	if(pdgd==11\|\|pdgd==13\|\|pdgd==211\|\|pdgd==321\|\|pdgd==2212){
396	fArrayDaugh->AddAt(jd,flastdaugh);
397	flastdaugh++;
398	}
399	else if(pdgd!=12&&pdgd!=14&&pdgd!=16&&pdgd!=310&&pdgd!=111){// CHECK NEPHEWS (FOR RESONANT CHANNELS)
400	Int_t nephw=pd->GetNDaughters();
401	if(nephw>0){
402	Int_t posN1=pd->GetDaughter(0);
403	Int_t posN2=pd->GetDaughter(1);
404	if(posN2<0)posN2=posN1;
405	if(posN1>=0){
406	for(Int_t jN=posN1;jN<=posN2;jN++){
407	AliAODMCParticle pN=(AliAODMCParticle)arrayMC->At(jN);
408	Int_t pdgN=TMath::Abs(pN->GetPdgCode());
409	if(pdgN==11\|\|pdgN==13\|\|pdgN==211\|\|pdgN==321\|\|pdgN==2212){
410	fArrayDaugh->AddAt(jN,flastdaugh);
411	flastdaugh++;
412	}
413	else if(pdgN!=12&&pdgN!=14&&pdgN!=16&&pdgN!=310&&pdgN!=111){// CHECK one step more (FOR RESONANT CHANNELS)
414
415	Int_t gnephw=pN->GetNDaughters();
416	if(gnephw>0){
417	Int_t posGN1=pN->GetDaughter(0);
418	Int_t posGN2=pN->GetDaughter(1);
419	if(posGN2<0)posGN2=posGN1;
420	if(posGN1>=0){
421	for(Int_t jGN=posGN1;jGN<=posGN2;jGN++){
422	AliAODMCParticle pGN=(AliAODMCParticle)arrayMC->At(jGN);
423	Int_t pdgGN=TMath::Abs(pGN->GetPdgCode());
424	if(pdgGN==11\|\|pdgGN==13\|\|pdgGN==211\|\|pdgGN==321\|\|pdgGN==2212){
425	fArrayDaugh->AddAt(jGN,flastdaugh);
426	flastdaugh++;
427	}
428	}
429	}
430	}
431	}
432	}
433	}
434	}
435	}
436	}
437	}
438	}
439	// printf("particles to be skipped %d\n",flastdaugh++);
440	}
441
442	void AliAnalysisTaskSEmcCorr::SelectAssociatedParticles(TClonesArray *arrayMC){
443	fLastAss=0;
444	fArrayAssoc->Reset(0);
445	Double_t pt,eta;
446	for(Int_t jAss=0;jAss<arrayMC->GetEntriesFast();jAss++){
447	AliAODMCParticle partAss=(AliAODMCParticle)arrayMC->At(jAss);
448	Int_t pdgAss=TMath::Abs(partAss->GetPdgCode());
449	if(pdgAss==11\|\|pdgAss==13\|\|pdgAss==211\|\|pdgAss==321\|\|pdgAss==2212){// check type: keep only e,mu,pi,K+,protons
450	pt=partAss->Pt();
451	eta=TMath::Abs(partAss->Eta());
452	if(pt<fminPtass)continue;
453	if(eta>fmaxEtaAss)continue;
454	if(!partAss->IsPhysicalPrimary()){
455	if(pdgAss!=11) continue;// keep ele
456	}
457	fArrayAssoc->AddAt(jAss,fLastAss);
458	fLastAss++;
459	}
460	}
461	// printf("n ass part: %d\n",fLastAss++);
462	}
463
464	void AliAnalysisTaskSEmcCorr::FillCorrelationPlots(AliAODMCParticle part,TClonesArray arrayMC){
465
466
467	Bool_t hasToSkip=kFALSE;
468	Int_t index=0,softpi=-1;
469	Int_t pdgTrig=TMath::Abs(part->GetPdgCode());
470	if(pdgTrig==421){
471	pdgTrig=1;
472	}
473	else if(pdgTrig==411){
474	pdgTrig=2;
475	}
476	else if(pdgTrig==413){
477	pdgTrig=3;
478	}
479	else if(pdgTrig==11){
480
481	pdgTrig=4;
482	}
483	else pdgTrig=0;
484	// check if it is prompt or from B and, for electrons the origin
485	Int_t labmum=part->GetMother();
486	if(labmum>=0){
487	AliAODMCParticle moth=(AliAODMCParticle)arrayMC->At(labmum);
488	Int_t pdgmum=TMath::Abs(moth->GetPdgCode());
489	if(pdgTrig==4){// for electrons go to the parent hadron
490	// labmum=part->GetMother();
491	// moth=(AliAODMCParticle*)arrayMC->At(labmum);
492	// pdgmum=TMath::Abs(moth->GetPdgCode());
493	if(!(pdgmum==5\|\|pdgmum==4\|\|(400<pdgmum&&pdgmum<600)\|\|(4000<pdgmum&&pdgmum<6000))){// NON HF electron
494	if(pdgmum==22){// GAMMA CONV
495	pdgTrig=5;
496	}
497	else if(pdgmum==111\|\|pdgmum==113\|\|pdgmum==221\|\|pdgmum==223){ //DALITZ DECAY
498	pdgTrig=6;
499	}
500	else pdgTrig=7;
501	}
502	}
503	if(pdgmum==413){ // D from Dstar: check soft pion and go to the grandmother
504	if(pdgTrig==1){
505	for(Int_t isp=moth->GetDaughter(0);isp<=moth->GetDaughter(1);isp++){
506	AliAODMCParticle sfp=(AliAODMCParticle)arrayMC->At(isp);
507	Int_t pdgsp=TMath::Abs(sfp->GetPdgCode());
508	if(pdgsp==211)softpi=isp;
509	}
510	}
511	labmum=moth->GetMother();
512	moth=(AliAODMCParticle*)arrayMC->At(labmum);
513	pdgmum=TMath::Abs(moth->GetPdgCode());
514	}
515	else if(pdgmum==423){// D*0 -> go to the mother
516	labmum=moth->GetMother();
517	moth=(AliAODMCParticle*)arrayMC->At(labmum);
518	pdgmum=TMath::Abs(moth->GetPdgCode());
519	}
520	if(pdgmum==5\|\|(500<pdgmum&&pdgmum<600)\|\|(5000<pdgmum&&pdgmum<6000)){
521	pdgTrig*=-1;
522	}
523	else if(pdgmum==4){
524	labmum=moth->GetMother();
525	if(labmum>=0){
526	moth=(AliAODMCParticle*)arrayMC->At(labmum);
527	pdgmum=TMath::Abs(moth->GetPdgCode());
528	if(pdgmum==5\|\|(500<pdgmum&&pdgmum<600)\|\|(5000<pdgmum&&pdgmum<6000)){
529	// Printf("c quark from b quark");
530	pdgTrig*=-1;
531	}
532	}
533	}
534	else {
535	// printf("Charm not from charm or beauty \n");
536	}
537	}
538	Double_t phiTrig=part->Phi();
539	Double_t etaTrig=part->Eta();
540	Double_t deltaPhi;
2942f542	541	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
d33db6d9 AR	542	fhMCtrigPart->Fill(point);
	543	for(Int_t jAss=0;jAss<fLastAss;jAss++){
	544	index=fArrayAssoc->At(jAss);
	545	hasToSkip=kFALSE;
	546	for(Int_t jd=0;jd<flastdaugh;jd++){
	547	if(index==fArrayDaugh->At(jd)){
	548	hasToSkip=kTRUE;
	549	break;
	550	}
	551	}
	552	if(hasToSkip)continue;
	553	AliAODMCParticle partAss=(AliAODMCParticle)arrayMC->At(index);
	554	point[3]=partAss->Pt();
	555	point[4]=partAss->Eta();
	556	deltaPhi=partAss->Phi()-phiTrig;
	557	if(deltaPhi<-0.5TMath::Pi())deltaPhi+=2.TMath::Pi();
	558	else if(deltaPhi>1.5TMath::Pi())deltaPhi-=2.TMath::Pi();
	559	point[5]=deltaPhi;
	560	point[6]=partAss->Eta()-etaTrig;
	561	UInt_t pdgassmum=TMath::Abs(partAss->GetPdgCode());
	562	if(index==softpi){
	563	point[7]=-1;
	564	}
	565	else if(pdgassmum==11){
	566	if(partAss->IsPhysicalPrimary())point[7]=1;
	567	Int_t elemum=partAss->GetMother();
	568	if(elemum>0){
	569	AliAODMCParticle partEleMum=(AliAODMCParticle)arrayMC->At(elemum);
	570	UInt_t pdgelemum=TMath::Abs(partEleMum->GetPdgCode());
	571	if(pdgelemum==111\|\|pdgelemum==113\|\|pdgelemum==221\|\|pdgelemum==223){// DALITZ DECAY
	572	point[7]=6;
	573	}
	574	else if(pdgelemum==22){// GAMMA CONV
	575	point[7]=7;
	576	}
	577	else if(pdgelemum==4\|\|pdgelemum==5\|\|(pdgelemum>400&&pdgelemum<600)\|\|(pdgelemum>4000&&pdgelemum<6000)){// HF e
	578	point[7]=8;
	579	}
	580	else point[7]=9;
	581	}
	582	}
	583	else if(pdgassmum==13){
	584	point[7]=2;
	585	}
	586	else if(pdgassmum==211){
	587	point[7]=3;
	588	}
	589	else if(pdgassmum==321){
	590	point[7]=4;
	591	}
	592	else if(pdgassmum==2212){
	593	point[7]=5;
	594	}
	595	else point[7]=0;
	596
	597	fhMCcorrel->Fill(point);
	598	}
	599
	600	}
	601
	602
	603
	604	//______________________________________________________________
	605	void AliAnalysisTaskSEmcCorr::FillCorrelationPlotsHadrons(AliAODMCParticle part,TClonesArray arrayMC){
606
607
608	Bool_t hasToSkip=kFALSE;
609	Int_t index=0;
610	Int_t pdgTrig=TMath::Abs(part->GetPdgCode());
611	if(pdgTrig==11){
612	pdgTrig=0;
613	}
614	else if(pdgTrig==13){
615	pdgTrig=1;
616	}
617	else if(pdgTrig==211){
618	pdgTrig=2;
619	}
620	else if(pdgTrig==321){
621	pdgTrig=3;
622	}
623	else if(pdgTrig==2212){
624	pdgTrig=4;
625	}
626	else pdgTrig=-1;
627
628
629	Double_t phiTrig=part->Phi();
630	Double_t etaTrig=part->Eta();
631	Double_t deltaPhi;
2942f542	632	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
d33db6d9 AR	633	fhMChadrontrigPart->Fill(point);
	634	for(Int_t jAss=0;jAss<fLastAss;jAss++){
	635	index=fArrayAssoc->At(jAss);
	636	hasToSkip=kFALSE;
	637	for(Int_t jd=0;jd<flastdaugh;jd++){
	638	if(index==fArrayDaugh->At(jd)){
	639	hasToSkip=kTRUE;
	640	break;
	641	}
	642	}
	643	if(hasToSkip)continue;
	644	AliAODMCParticle partAss=(AliAODMCParticle)arrayMC->At(index);
	645	point[3]=partAss->Pt();
	646	if(point[3]>point[1])continue;
	647	point[4]=partAss->Eta();
	648	deltaPhi=partAss->Phi()-phiTrig;
	649	if(deltaPhi<-0.5TMath::Pi())deltaPhi+=2.TMath::Pi();
	650	else if(deltaPhi>1.5TMath::Pi())deltaPhi-=2.TMath::Pi();
	651	point[5]=deltaPhi;
	652	point[6]=partAss->Eta()-etaTrig;
	653	point[7]=TMath::Abs(partAss->GetPdgCode());
	654	if(point[7]==11){
	655	point[7]=0;
	656	}
	657	else if(point[7]==13){
	658	point[7]=1;
	659	}
	660	else if(point[7]==211){
	661	point[7]=2;
	662	}
	663	else if(point[7]==321){
	664	point[7]=3;
	665	}
	666	else if(point[7]==2212){
	667	point[7]=4;
	668	}
	669	else point[7]=-1;
	670
	671	fhMChadroncorrel->Fill(point);
	672	}
	673
	674	}
	675
	676	//_______________________________________________________________
	677	void AliAnalysisTaskSEmcCorr::Terminate(const Option_t*){
	678	//TERMINATE METHOD: NOTHING TO DO
	679
	680
	681
	682	}
	683