[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliAnalysisTaskDielectronReadAODBranch.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 to retrieve the branch of the dielectron candidates stored in   //
//  filtered AODs file (AliAOD.Dielectron.root). It is possible to        //
//  apply tighter cuts to candidates stored in the branch and do the      //
//  matching with MC truth.                                               //
////////////////////////////////////////////////////////////////////////////

#include "TChain.h"
#include "TNtuple.h"
#include "TList.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TDatabasePDG.h"
#include "TF1.h"
#include "AliAnalysisManager.h"
#include "AliAODHandler.h"
#include "AliAODEvent.h"
#include "AliAODVertex.h"
#include "AliAODTrack.h"
#include "AliAODMCParticle.h"
#include "AliAnalysisTaskDielectronReadAODBranch.h"
#include "AliDielectronPair.h"
#include "AliDielectronMC.h"

ClassImp(AliAnalysisTaskDielectronReadAODBranch)
	//________________________________________________________________________
	AliAnalysisTaskDielectronReadAODBranch::AliAnalysisTaskDielectronReadAODBranch():
		AliAnalysisTaskSE(),
		fOutput(0),
                fNtupleJPSI(0),
                fNentries(0),
		fInvMass(0),
		fInvMassNoCuts(0),
                fpsproperSignal(0),
		fpsproperSidebands(0),            
		fpsproperAll(0),                  
		fpsproperUnder(0),
                fpsproperUpper(0),
                fLxyVsPtleg1(0),
		fLxyVsPtleg2(0),
		fLxyVsPt(0),
		fMeeVsPt(0),
		fMeeVsLxy(0),
		fprimvtxZ(0),  
		fsecvtxZ(0),  
		fprimvtxX(0),  
		fsecvtxX(0),  
		fprimvtxY(0),  
		fsecvtxY(0),  
		fPt(0),
		fPtLeg1(0),
		fPtLeg2(0),
		fdEdxP(0),
		fHasMC(0),
		fobj(0),
		fobjMC(0),
                fPtCut(1.),
                fSpdFirstRequired(kFALSE),
                fClsTPC(90),
                fPairType(1),
                fPtJpsi(1.3),
                fInvMassSignalLimits(0),
                fInvMassSideBandsLimits(0),
		fSecondary(0)
{
	// Default constructor
}

//________________________________________________________________________
AliAnalysisTaskDielectronReadAODBranch::AliAnalysisTaskDielectronReadAODBranch(const char *name):
	AliAnalysisTaskSE(name),
	fOutput(0),
        fNtupleJPSI(0),
        fNentries(0),
	fInvMass(0),
        fInvMassNoCuts(0),
        fpsproperSignal(0),
	fpsproperSidebands(0),            
	fpsproperAll(0),                  
	fpsproperUnder(0),
        fpsproperUpper(0),
        fLxyVsPtleg1(0),
	fLxyVsPtleg2(0),
	fLxyVsPt(0),
	fMeeVsPt(0),
	fMeeVsLxy(0),
	fprimvtxZ(0), 
	fsecvtxZ(0),  
	fprimvtxX(0), 
	fsecvtxX(0),  
	fprimvtxY(0), 
	fsecvtxY(0),  
	fPt(0),
	fPtLeg1(0),
	fPtLeg2(0),
	fdEdxP(0),
        fHasMC(0),
	fobj(0),
	fobjMC(0),
        fPtCut(1.),
        fSpdFirstRequired(kFALSE),
        fClsTPC(90),
        fPairType(1),
        fPtJpsi(1.3),
        fInvMassSignalLimits(0),
        fInvMassSideBandsLimits(0), 
	fSecondary(0)
{
	// Default constructor
	DefineInput(0,TChain::Class());
	DefineOutput(1,TList::Class());  //My private output

        fInvMassSignalLimits = new Double_t[2]; 
        fInvMassSignalLimits[0] = 0.;  fInvMassSignalLimits[1] = 0.; 
        fInvMassSideBandsLimits = new Double_t[2]; 
        fInvMassSideBandsLimits[0] = 0.;  fInvMassSideBandsLimits[1] = 0.;
}

//________________________________________________________________________
AliAnalysisTaskDielectronReadAODBranch::~AliAnalysisTaskDielectronReadAODBranch()
{
	if(fOutput)   delete fOutput;
	if(fobj)      delete fobj;
	if(fobjMC)    delete fobjMC;
}

//________________________________________________________________________
void AliAnalysisTaskDielectronReadAODBranch::Init()
{
	// Initialization
	if(fDebug > 1) printf("AnalysisTaskReadAOD::Init() \n");
	return;
}

//________________________________________________________________________
void AliAnalysisTaskDielectronReadAODBranch::UserCreateOutputObjects()
{
	//
        // Create the output container
	//
	if(fDebug > 1) printf("AnalysisTaskReadAOD::UserCreateOutputObjects() \n");
	fOutput = new TList();
	fOutput->SetOwner();
	//// Histogram booking
        
        // invariant mass
        fInvMass = new TH1F("fInvMass","fInvMass",300,2.0,2.0+300*.04); // step 40MeV
        fInvMassNoCuts = new TH1F("fInvMass_no_cuts","fInvMass_no_cuts",125,0,125*0.04); //step 40MeV
        fNentries = new TH1F("numberOfevent","numbOfEvent",1,-0.5,0.5);

        //pseudoproper 
	fpsproperSignal = new TH1F("psproper_decay_length",Form("psproper_decay_length_distrib(AODcuts+#clustTPC>%d, pT(leg)>%fGeV, pT(J/#psi)>%fGeV/c, %f < M < %f);X [#mu m];Entries/40#mu m",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]),150,-3000.,3000.);
	fpsproperSidebands = new TH1F("psproper_decay_length_sidebands",Form("psproper_decay_length_distrib_sidebands(AODcuts+#clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,M> %f && M < %f );  X [#mu m]; Entries/40#mu m",fClsTPC,fPtCut,fPtJpsi,fInvMassSideBandsLimits[1],fInvMassSideBandsLimits[0]),150,-3000.,3000.);
	fpsproperAll = new TH1F("psproper_decay_length_all",Form("psproper_decay_length_distrib_all(AODcuts+#clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c);X [#mu m];Entries/15#mu m",fClsTPC,fPtCut,fPtJpsi),400,-3000.,3000.);
	fpsproperUnder = new TH1F("psproper_decay_length_under",Form("psproper_decay_length_distrib_under(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c);X [#mu m];Entries/15#mu m",fClsTPC,fPtCut,fPtJpsi),400,-3000.,3000.); 
        fpsproperUpper = new TH1F("psproper_decay_length_upper",Form("psproper_decay_length_distrib_upper(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c);X [#mu m];Entries/15#mu m",fClsTPC,fPtCut,fPtJpsi),400,-3000.,3000.);  
        //
	fLxyVsPtleg1 = new TH2F("Lxy_vs_Pt_leg1",Form("Lxy_vs_Pt_leg1_distrib(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,%f<M<%f);p_{T}[GeV/c]; X",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]), 700, 0., 7.,400, -3000.,3000.);
	fLxyVsPtleg2 = new TH2F("Lxy_vs_Pt_leg2",Form("Lxy_vs_Pt_leg2_distrib(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,%f<M<%f);p_{T}[GeV/c]; X",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]), 700, 0., 7.,400, -3000.,3000.);
	fLxyVsPt = new TH2F("Lxy_vs_Pt_jpsi",Form("Lxy_vs_Pt_jpsi_distrib(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,%f<M<%f); p_{T}[GeV/c]; X",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]), 700, 0., 7.,400, -3000.,3000.);
	fMeeVsPt = new TH2F("Mee_vs_Pt_jpsi",Form("Mee_vs_Pt_jpsi_distrib(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,%f<M<%f); p_{T}[GeV/c]; M_{ee}",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]), 700, 0., 7.,200, 1.99,4.1);
	fMeeVsLxy = new TH2F("Mee_vs_Lxy_jpsi",Form("Mee_vs_Lxy_jpsi_distrib(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV),pT(J/#psi)>%f GeV/c; X; M_{ee}",fClsTPC,fPtCut,fPtJpsi), 400, -3000, 3000.,200, 1.99,4.1);

	// QA plots
	fprimvtxZ = new TH1F("prim_vtx_Z",Form("prim_vtx_Z_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV, pT(J/#psi)>%f);  Z[cm]; Entries/20 mm",fClsTPC,fPtCut,fPtJpsi),1000,-10.,10.);
	fsecvtxZ = new TH1F("sec_vtx_Z",Form("sec_vtx_Z_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV, pT(J/#psi)>%f);  Z[cm]; Entries/20 mm",fClsTPC,fPtCut,fPtJpsi),1000,-10.,10.);
	fprimvtxX = new TH1F("prim_vtx_X",Form("prim_vtx_X_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f);  X[cm]; Entries/mm",fClsTPC,fPtCut,fPtJpsi),1000,-0.5,0.5);
	fsecvtxX = new TH1F("sec_vtx_X",Form("sec_vtx_X_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f);  X[cm]; Entries/20 mm",fClsTPC,fPtCut,fPtJpsi),100,-1.,1.);
	fprimvtxY = new TH1F("prim_vtx_Y",Form("prim_vtx_Y_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f);  Y[cm]; Entries/mm",fClsTPC,fPtCut,fPtJpsi),1000,-0.5,0.5);
	fsecvtxY = new TH1F("sec_vtx_Y",Form("sec_vtx_Y_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f);  Y[cm]; Entries/20 mm",fClsTPC,fPtCut,fPtJpsi),100,-1.,1.);
	//
	fPt = new TH1F("Pt(J/psi)",Form("Pt_Jpsi_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV, pT(J/#psi)>%f); p_{T}(J/#psi)[GeV/c]; Entries/25 MeV",fClsTPC,fPtCut,fPtJpsi),400,0.,10.);
	fPtLeg1 = new TH1F("Pt_leg1",Form("Pt_leg1_distrib(AODcuts + #clustTPC>%d, pT(e+e-) > %f GeV, pT(J/#psi)>%f);  p_{T}[GeV/c]; Entries/25 MeV",fClsTPC,fPtCut,fPtJpsi),400,0.,10.);
	fPtLeg2 = new TH1F("Pt_leg2",Form("Pt_leg2_distrib(AODcuts + #clustTPC>%d, pT(e+e-) > %f GeV, pT(J/#psi)>%f);  p_{T}[GeV/c]; Entries/25 MeV",fClsTPC,fPtCut,fPtJpsi),400,0.,10.);
	//dE/dx plots
	fdEdxP = new TH2F("dE/dx_vs_Ptpc",Form("dE/dx_vs_Ptpc(AODcuts + #clustTPC>%d, pT(e+e-) > %f GeV, pT(J/#psi)>%f)",fClsTPC,fPtCut,fPtJpsi),400,0.2,20.,200,0.,200.);
        
        //J/psi n-tuple (X and M values)
        fNtupleJPSI = new TNtuple("fNtupleJPSI","J/#psi pseudo-proper decay time & invariant mass","Xdecaytime:Mass");
 
	// add histograms to list
	fOutput->Add(fNentries);
	fOutput->Add(fInvMass);
        fOutput->Add(fInvMassNoCuts); 
	fOutput->Add(fpsproperSignal);
	fOutput->Add(fpsproperSidebands);
	fOutput->Add(fpsproperAll);
	fOutput->Add(fpsproperUnder);
        fOutput->Add(fpsproperUpper);
        fOutput->Add(fLxyVsPtleg1);
	fOutput->Add(fLxyVsPtleg2);
	fOutput->Add(fLxyVsPt);
	fOutput->Add(fMeeVsPt);
	fOutput->Add(fMeeVsLxy);
	fOutput->Add(fprimvtxZ);
	fOutput->Add(fsecvtxZ);
	fOutput->Add(fprimvtxX);
	fOutput->Add(fsecvtxX);
	fOutput->Add(fprimvtxY);
	fOutput->Add(fsecvtxY);
	fOutput->Add(fPt);
	fOutput->Add(fPtLeg1);
	fOutput->Add(fPtLeg2);
        fOutput->Add(fNtupleJPSI);
        fOutput->Add(fdEdxP);
	return;
}

//________________________________________________________________________
void AliAnalysisTaskDielectronReadAODBranch::UserExec(Option_t */*option*/)
{
	// Execute analysis for current event:
	AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
	if (!aod) return;

	Double_t vtxPrim[3] = {0.,0.,0.};

	AliAODVertex* primvtx = aod->GetPrimaryVertex();
	vtxPrim[0] = primvtx->GetX();
	vtxPrim[1] = primvtx->GetY();
	vtxPrim[2] = primvtx->GetZ();

	AliAODHandler* aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());

	TTree *aodTree = aodHandler->GetTree();
	if(!fobj) aodTree->SetBranchAddress("dielectrons",&fobj);

	if(fHasMC){ aodTree->SetBranchAddress("mcparticles",&fobjMC);
	if(!fobjMC) printf("AliAnalysisTaskDielectronReadAODBranch::UserExec: MC particles branch not found!\n"); }

	fNentries->Fill(0);
	aodTree->GetEvent(Entry());

	// loop over candidates
	if(fobj) {
		TObjArray *objArr = (TObjArray*)fobj->UncheckedAt(fPairType);
		for(int j=0;j<objArr->GetEntriesFast();j++)
		{
			AliDielectronPair *pairObj = (AliDielectronPair*)objArr->UncheckedAt(j);
			fInvMassNoCuts->Fill(pairObj->M()); 
                        Double_t vtxSec[3] = {0.,0.,0.};
			fprimvtxX->Fill(vtxPrim[0]);
			fprimvtxY->Fill(vtxPrim[1]);
			fprimvtxZ->Fill(vtxPrim[2]);
			pairObj->XvYvZv(vtxSec);
			fsecvtxX->Fill(vtxSec[0]);
			fsecvtxY->Fill(vtxSec[1]);
			fsecvtxZ->Fill(vtxSec[2]);

			Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(pairObj->Px()) + (vtxSec[1]-vtxPrim[1])*(pairObj->Py()))/pairObj->Pt();
			Double_t psProperDecayLength = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/pairObj->Pt();

			AliAODTrack *trk = (AliAODTrack*)pairObj->GetFirstDaughterP();
			AliAODTrack *trk1 = (AliAODTrack*)pairObj->GetSecondDaughterP();
			if(!trk || !trk1) {printf("ERROR: daughter tracks not available\n"); continue;}

			//check in case of MC analysis if candidate is a true J/psi->ee
                          if(fHasMC && ((AliDielectronMC::Instance()->IsJpsiPrimary(pairObj))!=fSecondary)) continue;
			   AliAODPid *pid  = trk->GetDetPid();
                                AliAODPid *pid1 = trk1->GetDetPid();
                                if(!pid || !pid1){
                                        if(fDebug>1) printf("No AliAODPid found\n");
                                        continue;
                                }
                                 // ptLeg cut
                                if((trk->Pt()<fPtCut) || (trk1->Pt()<fPtCut)) continue;
                           
                                // pt jpsi cut
                                if((pairObj->Pt())<fPtJpsi) continue; 
 
                              //spd first required
                              if(fSpdFirstRequired)
                              {
                               if((!trk->HasPointOnITSLayer(0)) || (!trk1->HasPointOnITSLayer(0))) continue;
                              }
                           //
                          if((trk->GetTPCNcls()) <= fClsTPC || (trk1->GetTPCNcls()) <= fClsTPC) continue;
                         	if(trk->Charge()==-1){
					fPtLeg1->Fill(trk->Pt());
				} else {fPtLeg2->Fill(trk->Pt());}
				if(trk1->Charge()==-1){
					fPtLeg1->Fill(trk1->Pt());
				} else {fPtLeg2->Fill(trk1->Pt());}
				//Fill dE/dx related plots (before PID cuts)
				fdEdxP->Fill(pid->GetTPCmomentum(),pid->GetTPCsignal());
				fdEdxP->Fill(pid1->GetTPCmomentum(),pid1->GetTPCsignal());
					fInvMass->Fill(pairObj->M());
				        fPt->Fill(pairObj->Pt());
                                 	fMeeVsLxy->Fill(10000*psProperDecayLength,pairObj->M());
					fMeeVsPt->Fill(pairObj->Pt(),pairObj->M());
					fpsproperAll->Fill(10000*psProperDecayLength);
                                           
                                           //psproper in the signal region 
					   if((pairObj->M())<fInvMassSignalLimits[1] && (pairObj->M())>fInvMassSignalLimits[0]){
                                         	fpsproperSignal->Fill(10000*psProperDecayLength);
						fLxyVsPt->Fill(pairObj->Pt(),10000*psProperDecayLength);  //jpsi 
						fLxyVsPtleg1->Fill(trk->Pt(),10000*psProperDecayLength);  //leg1 (warning: mixture of pos and neg tracks)
						fLxyVsPtleg2->Fill(trk1->Pt(),10000*psProperDecayLength); //leg2 (warning: mixture of pos and neg tracks)
					        fNtupleJPSI->Fill(10000*psProperDecayLength,pairObj->M()); // fill the N-tuple (imput of the minimization algorithm)
                                                 } 

                                        // psproper in the sidebands 
					if((pairObj->M())> fInvMassSideBandsLimits[1] || (pairObj->M())< fInvMassSideBandsLimits[0])
                                        fpsproperSidebands->Fill(10000*psProperDecayLength);
                                          
                                        // psproper in the lower sideband
				        if(pairObj->M()< fInvMassSideBandsLimits[0]) fpsproperUnder->Fill(10000*psProperDecayLength);
                                          
                                        // psproper in the upper sideband
                                        if(pairObj->M()> fInvMassSideBandsLimits[1]) fpsproperUpper->Fill(10000*psProperDecayLength);
                        } // end loop over candidates
		}// end loop over pair types
        // Post the data
	PostData(1,fOutput);
	return;
} 

//________________________________________________________________________
void AliAnalysisTaskDielectronReadAODBranch::Terminate(Option_t */*option*/)
{
	if(fDebug > 1) printf("AliAnalysisTaskDielectronReadAODBranch::Terminate() \n");
	return;
}
Commit	Line	Data
ba15fdfb	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	// Class to retrieve the branch of the dielectron candidates stored in //
	18	// filtered AODs file (AliAOD.Dielectron.root). It is possible to //
	19	// apply tighter cuts to candidates stored in the branch and do the //
	20	// matching with MC truth. //
	21	////////////////////////////////////////////////////////////////////////////
	22
	23	#include "TChain.h"
	24	#include "TNtuple.h"
bc7a3220	25	#include "TList.h"
ba15fdfb	26	#include "TH1F.h"
	27	#include "TH2F.h"
	28	#include "TDatabasePDG.h"
	29	#include "TF1.h"
	30	#include "AliAnalysisManager.h"
	31	#include "AliAODHandler.h"
	32	#include "AliAODEvent.h"
	33	#include "AliAODVertex.h"
	34	#include "AliAODTrack.h"
	35	#include "AliAODMCParticle.h"
	36	#include "AliAnalysisTaskDielectronReadAODBranch.h"
	37	#include "AliDielectronPair.h"
bc7a3220	38	#include "AliDielectronMC.h"
ba15fdfb	39
	40	ClassImp(AliAnalysisTaskDielectronReadAODBranch)
	41	//________________________________________________________________________
	42	AliAnalysisTaskDielectronReadAODBranch::AliAnalysisTaskDielectronReadAODBranch():
	43	AliAnalysisTaskSE(),
	44	fOutput(0),
	45	fNtupleJPSI(0),
	46	fNentries(0),
	47	fInvMass(0),
	48	fInvMassNoCuts(0),
	49	fpsproperSignal(0),
	50	fpsproperSidebands(0),
	51	fpsproperAll(0),
	52	fpsproperUnder(0),
	53	fpsproperUpper(0),
	54	fLxyVsPtleg1(0),
	55	fLxyVsPtleg2(0),
	56	fLxyVsPt(0),
	57	fMeeVsPt(0),
	58	fMeeVsLxy(0),
	59	fprimvtxZ(0),
	60	fsecvtxZ(0),
	61	fprimvtxX(0),
	62	fsecvtxX(0),
	63	fprimvtxY(0),
	64	fsecvtxY(0),
	65	fPt(0),
	66	fPtLeg1(0),
	67	fPtLeg2(0),
	68	fdEdxP(0),
	69	fHasMC(0),
	70	fobj(0),
	71	fobjMC(0),
	72	fPtCut(1.),
	73	fSpdFirstRequired(kFALSE),
	74	fClsTPC(90),
	75	fPairType(1),
	76	fPtJpsi(1.3),
	77	fInvMassSignalLimits(0),
bc7a3220	78	fInvMassSideBandsLimits(0),
bc7a3220	79	fSecondary(0)
ba15fdfb	80	{
	81	// Default constructor
	82	}
	83
	84	//________________________________________________________________________
	85	AliAnalysisTaskDielectronReadAODBranch::AliAnalysisTaskDielectronReadAODBranch(const char *name):
	86	AliAnalysisTaskSE(name),
	87	fOutput(0),
	88	fNtupleJPSI(0),
	89	fNentries(0),
	90	fInvMass(0),
	91	fInvMassNoCuts(0),
	92	fpsproperSignal(0),
	93	fpsproperSidebands(0),
	94	fpsproperAll(0),
	95	fpsproperUnder(0),
	96	fpsproperUpper(0),
	97	fLxyVsPtleg1(0),
	98	fLxyVsPtleg2(0),
	99	fLxyVsPt(0),
	100	fMeeVsPt(0),
	101	fMeeVsLxy(0),
	102	fprimvtxZ(0),
	103	fsecvtxZ(0),
	104	fprimvtxX(0),
	105	fsecvtxX(0),
	106	fprimvtxY(0),
	107	fsecvtxY(0),
	108	fPt(0),
	109	fPtLeg1(0),
	110	fPtLeg2(0),
	111	fdEdxP(0),
	112	fHasMC(0),
	113	fobj(0),
	114	fobjMC(0),
	115	fPtCut(1.),
	116	fSpdFirstRequired(kFALSE),
	117	fClsTPC(90),
	118	fPairType(1),
	119	fPtJpsi(1.3),
	120	fInvMassSignalLimits(0),
bc7a3220	121	fInvMassSideBandsLimits(0),
bc7a3220	122	fSecondary(0)
ba15fdfb	123	{
	124	// Default constructor
	125	DefineInput(0,TChain::Class());
	126	DefineOutput(1,TList::Class()); //My private output
	127
	128	fInvMassSignalLimits = new Double_t[2];
	129	fInvMassSignalLimits[0] = 0.; fInvMassSignalLimits[1] = 0.;
	130	fInvMassSideBandsLimits = new Double_t[2];
	131	fInvMassSideBandsLimits[0] = 0.; fInvMassSideBandsLimits[1] = 0.;
	132	}
	133
	134	//________________________________________________________________________
	135	AliAnalysisTaskDielectronReadAODBranch::~AliAnalysisTaskDielectronReadAODBranch()
	136	{
	137	if(fOutput) delete fOutput;
	138	if(fobj) delete fobj;
	139	if(fobjMC) delete fobjMC;
	140	}
	141
	142	//________________________________________________________________________
	143	void AliAnalysisTaskDielectronReadAODBranch::Init()
	144	{
	145	// Initialization
	146	if(fDebug > 1) printf("AnalysisTaskReadAOD::Init() \n");
	147	return;
	148	}
	149
	150	//________________________________________________________________________
	151	void AliAnalysisTaskDielectronReadAODBranch::UserCreateOutputObjects()
	152	{
	153	//
	154	// Create the output container
	155	//
	156	if(fDebug > 1) printf("AnalysisTaskReadAOD::UserCreateOutputObjects() \n");
	157	fOutput = new TList();
	158	fOutput->SetOwner();
	159	//// Histogram booking
	160
	161	// invariant mass
	162	fInvMass = new TH1F("fInvMass","fInvMass",300,2.0,2.0+300*.04); // step 40MeV
	163	fInvMassNoCuts = new TH1F("fInvMass_no_cuts","fInvMass_no_cuts",125,0,125*0.04); //step 40MeV
	164	fNentries = new TH1F("numberOfevent","numbOfEvent",1,-0.5,0.5);
bc7a3220	165
ba15fdfb	166	//pseudoproper
	167	fpsproperSignal = new TH1F("psproper_decay_length",Form("psproper_decay_length_distrib(AODcuts+#clustTPC>%d, pT(leg)>%fGeV, pT(J/#psi)>%fGeV/c, %f < M < %f);X [#mu m];Entries/40#mu m",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]),150,-3000.,3000.);
	168	fpsproperSidebands = new TH1F("psproper_decay_length_sidebands",Form("psproper_decay_length_distrib_sidebands(AODcuts+#clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,M> %f && M < %f ); X [#mu m]; Entries/40#mu m",fClsTPC,fPtCut,fPtJpsi,fInvMassSideBandsLimits[1],fInvMassSideBandsLimits[0]),150,-3000.,3000.);
	169	fpsproperAll = new TH1F("psproper_decay_length_all",Form("psproper_decay_length_distrib_all(AODcuts+#clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c);X [#mu m];Entries/15#mu m",fClsTPC,fPtCut,fPtJpsi),400,-3000.,3000.);
	170	fpsproperUnder = new TH1F("psproper_decay_length_under",Form("psproper_decay_length_distrib_under(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c);X [#mu m];Entries/15#mu m",fClsTPC,fPtCut,fPtJpsi),400,-3000.,3000.);
	171	fpsproperUpper = new TH1F("psproper_decay_length_upper",Form("psproper_decay_length_distrib_upper(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c);X [#mu m];Entries/15#mu m",fClsTPC,fPtCut,fPtJpsi),400,-3000.,3000.);
	172	//
	173	fLxyVsPtleg1 = new TH2F("Lxy_vs_Pt_leg1",Form("Lxy_vs_Pt_leg1_distrib(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,%f<M<%f);p_{T}[GeV/c]; X",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]), 700, 0., 7.,400, -3000.,3000.);
	174	fLxyVsPtleg2 = new TH2F("Lxy_vs_Pt_leg2",Form("Lxy_vs_Pt_leg2_distrib(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,%f<M<%f);p_{T}[GeV/c]; X",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]), 700, 0., 7.,400, -3000.,3000.);
	175	fLxyVsPt = new TH2F("Lxy_vs_Pt_jpsi",Form("Lxy_vs_Pt_jpsi_distrib(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,%f<M<%f); p_{T}[GeV/c]; X",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]), 700, 0., 7.,400, -3000.,3000.);
	176	fMeeVsPt = new TH2F("Mee_vs_Pt_jpsi",Form("Mee_vs_Pt_jpsi_distrib(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,%f<M<%f); p_{T}[GeV/c]; M_{ee}",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]), 700, 0., 7.,200, 1.99,4.1);
	177	fMeeVsLxy = new TH2F("Mee_vs_Lxy_jpsi",Form("Mee_vs_Lxy_jpsi_distrib(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV),pT(J/#psi)>%f GeV/c; X; M_{ee}",fClsTPC,fPtCut,fPtJpsi), 400, -3000, 3000.,200, 1.99,4.1);
	178
	179	// QA plots
	180	fprimvtxZ = new TH1F("prim_vtx_Z",Form("prim_vtx_Z_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV, pT(J/#psi)>%f); Z[cm]; Entries/20 mm",fClsTPC,fPtCut,fPtJpsi),1000,-10.,10.);
	181	fsecvtxZ = new TH1F("sec_vtx_Z",Form("sec_vtx_Z_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV, pT(J/#psi)>%f); Z[cm]; Entries/20 mm",fClsTPC,fPtCut,fPtJpsi),1000,-10.,10.);
	182	fprimvtxX = new TH1F("prim_vtx_X",Form("prim_vtx_X_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f); X[cm]; Entries/mm",fClsTPC,fPtCut,fPtJpsi),1000,-0.5,0.5);
	183	fsecvtxX = new TH1F("sec_vtx_X",Form("sec_vtx_X_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f); X[cm]; Entries/20 mm",fClsTPC,fPtCut,fPtJpsi),100,-1.,1.);
	184	fprimvtxY = new TH1F("prim_vtx_Y",Form("prim_vtx_Y_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f); Y[cm]; Entries/mm",fClsTPC,fPtCut,fPtJpsi),1000,-0.5,0.5);
	185	fsecvtxY = new TH1F("sec_vtx_Y",Form("sec_vtx_Y_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f); Y[cm]; Entries/20 mm",fClsTPC,fPtCut,fPtJpsi),100,-1.,1.);
	186	//
	187	fPt = new TH1F("Pt(J/psi)",Form("Pt_Jpsi_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV, pT(J/#psi)>%f); p_{T}(J/#psi)[GeV/c]; Entries/25 MeV",fClsTPC,fPtCut,fPtJpsi),400,0.,10.);
	188	fPtLeg1 = new TH1F("Pt_leg1",Form("Pt_leg1_distrib(AODcuts + #clustTPC>%d, pT(e+e-) > %f GeV, pT(J/#psi)>%f); p_{T}[GeV/c]; Entries/25 MeV",fClsTPC,fPtCut,fPtJpsi),400,0.,10.);
	189	fPtLeg2 = new TH1F("Pt_leg2",Form("Pt_leg2_distrib(AODcuts + #clustTPC>%d, pT(e+e-) > %f GeV, pT(J/#psi)>%f); p_{T}[GeV/c]; Entries/25 MeV",fClsTPC,fPtCut,fPtJpsi),400,0.,10.);
	190	//dE/dx plots
	191	fdEdxP = new TH2F("dE/dx_vs_Ptpc",Form("dE/dx_vs_Ptpc(AODcuts + #clustTPC>%d, pT(e+e-) > %f GeV, pT(J/#psi)>%f)",fClsTPC,fPtCut,fPtJpsi),400,0.2,20.,200,0.,200.);
	192
	193	//J/psi n-tuple (X and M values)
	194	fNtupleJPSI = new TNtuple("fNtupleJPSI","J/#psi pseudo-proper decay time & invariant mass","Xdecaytime:Mass");
	195
	196	// add histograms to list
	197	fOutput->Add(fNentries);
	198	fOutput->Add(fInvMass);
	199	fOutput->Add(fInvMassNoCuts);
	200	fOutput->Add(fpsproperSignal);
	201	fOutput->Add(fpsproperSidebands);
	202	fOutput->Add(fpsproperAll);
	203	fOutput->Add(fpsproperUnder);
	204	fOutput->Add(fpsproperUpper);
	205	fOutput->Add(fLxyVsPtleg1);
	206	fOutput->Add(fLxyVsPtleg2);
	207	fOutput->Add(fLxyVsPt);
	208	fOutput->Add(fMeeVsPt);
	209	fOutput->Add(fMeeVsLxy);
	210	fOutput->Add(fprimvtxZ);
	211	fOutput->Add(fsecvtxZ);
	212	fOutput->Add(fprimvtxX);
	213	fOutput->Add(fsecvtxX);
	214	fOutput->Add(fprimvtxY);
	215	fOutput->Add(fsecvtxY);
	216	fOutput->Add(fPt);
	217	fOutput->Add(fPtLeg1);
	218	fOutput->Add(fPtLeg2);
	219	fOutput->Add(fNtupleJPSI);
	220	fOutput->Add(fdEdxP);
	221	return;
	222	}
	223
	224	//________________________________________________________________________
	225	void AliAnalysisTaskDielectronReadAODBranch::UserExec(Option_t /option*/)
	226	{
	227	// Execute analysis for current event:
	228	AliAODEvent aod = dynamic_cast<AliAODEvent> (InputEvent());
c8f0f810	229	if (!aod) return;
bc7a3220	230
ba15fdfb	231	Double_t vtxPrim[3] = {0.,0.,0.};
	232
	233	AliAODVertex* primvtx = aod->GetPrimaryVertex();
	234	vtxPrim[0] = primvtx->GetX();
	235	vtxPrim[1] = primvtx->GetY();
	236	vtxPrim[2] = primvtx->GetZ();
	237
	238	AliAODHandler* aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
	239
	240	TTree *aodTree = aodHandler->GetTree();
	241	if(!fobj) aodTree->SetBranchAddress("dielectrons",&fobj);
	242
	243	if(fHasMC){ aodTree->SetBranchAddress("mcparticles",&fobjMC);
	244	if(!fobjMC) printf("AliAnalysisTaskDielectronReadAODBranch::UserExec: MC particles branch not found!\n"); }
	245
	246	fNentries->Fill(0);
	247	aodTree->GetEvent(Entry());
ba15fdfb	248
	249	// loop over candidates
	250	if(fobj) {
	251	TObjArray objArr = (TObjArray)fobj->UncheckedAt(fPairType);
	252	for(int j=0;j<objArr->GetEntriesFast();j++)
	253	{
	254	AliDielectronPair pairObj = (AliDielectronPair)objArr->UncheckedAt(j);
	255	fInvMassNoCuts->Fill(pairObj->M());
	256	Double_t vtxSec[3] = {0.,0.,0.};
	257	fprimvtxX->Fill(vtxPrim[0]);
	258	fprimvtxY->Fill(vtxPrim[1]);
	259	fprimvtxZ->Fill(vtxPrim[2]);
	260	pairObj->XvYvZv(vtxSec);
	261	fsecvtxX->Fill(vtxSec[0]);
	262	fsecvtxY->Fill(vtxSec[1]);
	263	fsecvtxZ->Fill(vtxSec[2]);
	264
	265	Double_t lxy = ((vtxSec[0]-vtxPrim[0])(pairObj->Px()) + (vtxSec[1]-vtxPrim[1])(pairObj->Py()))/pairObj->Pt();
	266	Double_t psProperDecayLength = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/pairObj->Pt();
	267
54ce9dc4	268	AliAODTrack trk = (AliAODTrack)pairObj->GetFirstDaughterP();
54ce9dc4	269	AliAODTrack trk1 = (AliAODTrack)pairObj->GetSecondDaughterP();
ba15fdfb	270	if(!trk \|\| !trk1) {printf("ERROR: daughter tracks not available\n"); continue;}
ba15fdfb	271
ba15fdfb	272	//check in case of MC analysis if candidate is a true J/psi->ee
bc7a3220	273	if(fHasMC && ((AliDielectronMC::Instance()->IsJpsiPrimary(pairObj))!=fSecondary)) continue;
bc7a3220	274	AliAODPid *pid = trk->GetDetPid();
ba15fdfb	275	AliAODPid *pid1 = trk1->GetDetPid();
	276	if(!pid \|\| !pid1){
	277	if(fDebug>1) printf("No AliAODPid found\n");
	278	continue;
	279	}
	280	// ptLeg cut
	281	if((trk->Pt()<fPtCut) \|\| (trk1->Pt()<fPtCut)) continue;
	282
	283	// pt jpsi cut
	284	if((pairObj->Pt())<fPtJpsi) continue;
	285
	286	//spd first required
	287	if(fSpdFirstRequired)
	288	{
	289	if((!trk->HasPointOnITSLayer(0)) \|\| (!trk1->HasPointOnITSLayer(0))) continue;
	290	}
	291	//
	292	if((trk->GetTPCNcls()) <= fClsTPC \|\| (trk1->GetTPCNcls()) <= fClsTPC) continue;
	293	if(trk->Charge()==-1){
	294	fPtLeg1->Fill(trk->Pt());
	295	} else {fPtLeg2->Fill(trk->Pt());}
	296	if(trk1->Charge()==-1){
	297	fPtLeg1->Fill(trk1->Pt());
	298	} else {fPtLeg2->Fill(trk1->Pt());}
	299	//Fill dE/dx related plots (before PID cuts)
	300	fdEdxP->Fill(pid->GetTPCmomentum(),pid->GetTPCsignal());
	301	fdEdxP->Fill(pid1->GetTPCmomentum(),pid1->GetTPCsignal());
	302	fInvMass->Fill(pairObj->M());
	303	fPt->Fill(pairObj->Pt());
	304	fMeeVsLxy->Fill(10000*psProperDecayLength,pairObj->M());
	305	fMeeVsPt->Fill(pairObj->Pt(),pairObj->M());
	306	fpsproperAll->Fill(10000*psProperDecayLength);
	307
	308	//psproper in the signal region
	309	if((pairObj->M())<fInvMassSignalLimits[1] && (pairObj->M())>fInvMassSignalLimits[0]){
	310	fpsproperSignal->Fill(10000*psProperDecayLength);
	311	fLxyVsPt->Fill(pairObj->Pt(),10000*psProperDecayLength); //jpsi
	312	fLxyVsPtleg1->Fill(trk->Pt(),10000*psProperDecayLength); //leg1 (warning: mixture of pos and neg tracks)
	313	fLxyVsPtleg2->Fill(trk1->Pt(),10000*psProperDecayLength); //leg2 (warning: mixture of pos and neg tracks)
	314	fNtupleJPSI->Fill(10000*psProperDecayLength,pairObj->M()); // fill the N-tuple (imput of the minimization algorithm)
	315	}
	316
	317	// psproper in the sidebands
	318	if((pairObj->M())> fInvMassSideBandsLimits[1] \|\| (pairObj->M())< fInvMassSideBandsLimits[0])
	319	fpsproperSidebands->Fill(10000*psProperDecayLength);
	320
	321	// psproper in the lower sideband
	322	if(pairObj->M()< fInvMassSideBandsLimits[0]) fpsproperUnder->Fill(10000*psProperDecayLength);
	323
	324	// psproper in the upper sideband
	325	if(pairObj->M()> fInvMassSideBandsLimits[1]) fpsproperUpper->Fill(10000*psProperDecayLength);
	326	} // end loop over candidates
	327	}// end loop over pair types
	328	// Post the data
	329	PostData(1,fOutput);
	330	return;
	331	}
	332
	333	//________________________________________________________________________
	334	void AliAnalysisTaskDielectronReadAODBranch::Terminate(Option_t /option*/)
	335	{
	336	if(fDebug > 1) printf("AliAnalysisTaskDielectronReadAODBranch::Terminate() \n");
	337	return;
	338	}
339