[u/mrichter/AliRoot.git] / PWGGA / GammaConv / AliAnalysisTaskMaterial.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *																								*
 * Authors: Friederike Bock															*
 * Version 1.0																				*
 *																								*
 * 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.						*
 **************************************************************************/

////////////////////////////////////////////////
//---------------------------------------------
// QA Task for V0 Reader V1
//---------------------------------------------
////////////////////////////////////////////////

#include "AliAnalysisTaskMaterial.h"
#include "TChain.h"
#include "AliAnalysisManager.h"
#include "TParticle.h"
#include "TVectorF.h"
#include "AliPIDResponse.h"
#include "AliESDtrackCuts.h"
#include "TFile.h"

class iostream;

using namespace std;

ClassImp(AliAnalysisTaskMaterial)

AliAnalysisTaskMaterial::AliAnalysisTaskMaterial() : AliAnalysisTaskSE(),
   fV0Reader(NULL),
   fConversionGammas(NULL),
   fConversionCuts(NULL),
   fOutputList(NULL),
   fEventList(NULL),
   fRecGammaList(NULL),
   fAllMCGammaList(NULL),
   fAllMCConvGammaList(NULL),
   fTreeEvent(NULL),
   fTreeMaterialRec(NULL),
   fTreeMaterialAllGamma(NULL),
   fTreeMaterialConvGamma(NULL),
   fPrimVtxZ(0.),
   fNContrVtx(0),
   fNESDtracksEta09(0),
   fNESDtracksEta0914(0),
   fNESDtracksEta14(0),
   fGammaMCPt(0.),
   fGammaMCTheta(0.),
   fGammaMCConvPt(0.),
   fGammaMCConvTheta(0.),
   fMCConvCords(5),
   fMCConvDaughterProp(4),
   fGammaPt(0.),
   fGammaTheta(0.),
   fGammaChi2NDF(0.),
   fRecCords(5),
   fDaughterProp(4),
   fKind(0),      
   fIsHeavyIon(kFALSE),
   fIsMC(kFALSE),
   fESDEvent(NULL),
   fMCEvent(NULL)
{

   
}


//________________________________________________________________________
AliAnalysisTaskMaterial::AliAnalysisTaskMaterial(const char *name) : AliAnalysisTaskSE(name),
   fV0Reader(NULL),
   fConversionGammas(NULL),
   fConversionCuts(NULL),
   fOutputList(NULL),
   fEventList(NULL),
   fRecGammaList(NULL),
   fAllMCGammaList(NULL),
   fAllMCConvGammaList(NULL),
   fTreeEvent(NULL),
   fTreeMaterialRec(NULL),
   fTreeMaterialAllGamma(NULL),
   fTreeMaterialConvGamma(NULL),
   fPrimVtxZ(0.),
   fNContrVtx(0),
   fNESDtracksEta09(0),
   fNESDtracksEta0914(0),
   fNESDtracksEta14(0),
   fGammaMCPt(0.),
   fGammaMCTheta(0.),
   fGammaMCConvPt(0.),
   fGammaMCConvTheta(0.),
   fMCConvCords(5),
   fMCConvDaughterProp(4),
   fGammaPt(0.),
   fGammaTheta(0.),
   fGammaChi2NDF(0.),
   fRecCords(5),
   fDaughterProp(4),
   fKind(0),      
   fIsHeavyIon(kFALSE),
   fIsMC(kFALSE),
   fESDEvent(NULL),
   fMCEvent(NULL)
{
   // Default constructor

   
   DefineInput(0, TChain::Class());
   DefineOutput(1, TList::Class());
}

//________________________________________________________________________
AliAnalysisTaskMaterial::~AliAnalysisTaskMaterial()
{
   // default deconstructor
}
//________________________________________________________________________
void AliAnalysisTaskMaterial::UserCreateOutputObjects()
{
   // Create User Output Objects

   if(fOutputList != NULL){
      delete fOutputList;
      fOutputList = NULL;
   }
   if(fOutputList == NULL){
      fOutputList = new TList();
      fOutputList->SetOwner(kTRUE);
   }
   
   fEventList = new TList();
   fEventList->SetName("EventList");
   fEventList->SetOwner(kTRUE);
   fOutputList->Add(fEventList);
   
   fTreeEvent = new TTree("Event","Event");   
   fTreeEvent->Branch("primVtxZ",&fPrimVtxZ,"fPrimVtxZ/F");
   fTreeEvent->Branch("nContrVtx",&fNContrVtx,"fNContrVtx/I");
   fTreeEvent->Branch("nGoodTracksEta09",&fNESDtracksEta09,"fNESDtracksEta09/I");
   fTreeEvent->Branch("nGoodTracksEta0914",&fNESDtracksEta0914,"fNESDtracksEta0914/I");
   fTreeEvent->Branch("nGoodTracksEta14",&fNESDtracksEta14,"fNESDtracksEta14/I");
   fEventList->Add(fTreeEvent);
   
   fRecGammaList= new TList();
   fRecGammaList->SetName("RecGammaList");
   fRecGammaList->SetOwner(kTRUE);
   fOutputList->Add(fRecGammaList);
   
   
   fTreeMaterialRec = new TTree("ConvPointRec","ConvPointRec");   
   fTreeMaterialRec->Branch("recCords",&fRecCords);
   fTreeMaterialRec->Branch("daughterProp",&fDaughterProp);
   fTreeMaterialRec->Branch("pt",&fGammaPt,"fGammaPt/F");
   fTreeMaterialRec->Branch("theta",&fGammaTheta,"fGammaTheta/F");
   fTreeMaterialRec->Branch("chi2ndf",&fGammaChi2NDF,"fGammaChi2NDF/F");
   if (fIsMC) {
      fTreeMaterialRec->Branch("kind",&fKind,"fKind/b");
   }   
   fRecGammaList->Add(fTreeMaterialRec);
   
   if (fIsMC) {
      fAllMCGammaList = new TList();
      fAllMCGammaList->SetName("AllMCGammaList");
      fAllMCGammaList->SetOwner(kTRUE);
      fOutputList->Add(fAllMCGammaList);
      
      fTreeMaterialAllGamma = new TTree("AllGamma","AllGamma");   
      fTreeMaterialAllGamma->Branch("pt",&fGammaMCPt,"fGammaMCPt/F");
      fTreeMaterialAllGamma->Branch("theta",&fGammaMCTheta,"fGammaMCTheta/I");
      fAllMCGammaList->Add(fTreeMaterialAllGamma);
      
      fAllMCConvGammaList = new TList();
      fAllMCConvGammaList->SetName("AllMCGammaConvList");
      fAllMCConvGammaList->SetOwner(kTRUE);
      fOutputList->Add(fAllMCConvGammaList);

//       fMCConvCords = new Float_t[5];
//       fMCConvDaughterProp = new Float_t[4];

      
      fTreeMaterialConvGamma = new TTree("ConvGammaMC","ConvGammaMC");   
      fTreeMaterialConvGamma->Branch("Cords",&fMCConvCords);
      fTreeMaterialConvGamma->Branch("daughterProp",&fMCConvDaughterProp);
      fTreeMaterialConvGamma->Branch("Pt",&fGammaMCConvPt,"fGammaMCConvPt/F");
      fTreeMaterialConvGamma->Branch("Theta",&fGammaMCConvTheta,"fGammaMCConvTheta/F");   
      fAllMCConvGammaList->Add(fTreeMaterialConvGamma);
   }
   
   // V0 Reader Cuts
   TString cutnumber = fConversionCuts->GetCutNumber();
   PostData(1, fOutputList);
   
}

//________________________________________________________________________
void AliAnalysisTaskMaterial::UserExec(Option_t *){

   fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1");

   Int_t eventQuality = ((AliConversionCuts*)fV0Reader->GetConversionCuts())->GetEventQuality();
	if(eventQuality != 0){// Event Not Accepted
      return;
   }
   fESDEvent = (AliESDEvent*) InputEvent();
   if (fESDEvent==NULL) return;
   if(fIsHeavyIon && !fConversionCuts->IsCentralitySelected(fESDEvent)) return;
	fNESDtracksEta09 = CountTracks09(); // Estimate Event Multiplicity
	fNESDtracksEta0914 = CountTracks0914(); // Estimate Event Multiplicity
	fNESDtracksEta14 = fNESDtracksEta09 + fNESDtracksEta0914;
	if(fESDEvent){
		if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()>0) {
			fNContrVtx = fESDEvent->GetPrimaryVertexTracks()->GetNContributors();
		} else if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()<1) {
			if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
				fNContrVtx = fESDEvent->GetPrimaryVertexSPD()->GetNContributors();
			} else if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()<1) {
				fNContrVtx = 0;
			}
		}
   }
	fPrimVtxZ = fESDEvent->GetPrimaryVertex()->GetZ();
	
	if (fTreeEvent){
      fTreeEvent->Fill();
	}
	
   fConversionGammas=fV0Reader->GetReconstructedGammas();
	if(MCEvent()){
		fMCEvent = MCEvent();
	}
   ProcessPhotons();
	if(MCEvent()){
		ProcessMCPhotons();
	}
   PostData(1, fOutputList);
}

///________________________________________________________________________
void AliAnalysisTaskMaterial::FillMCTree(Int_t stackPos){
	AliStack *MCStack = fMCEvent->Stack();
	TParticle* candidate = (TParticle *)MCStack->Particle(stackPos);
	if(fConversionCuts->PhotonIsSelectedMC(candidate,MCStack,kFALSE)){
		fGammaMCPt = candidate->Pt();
		fGammaMCTheta = candidate->Theta();
		if (fTreeMaterialAllGamma){
			fTreeMaterialAllGamma->Fill();
		}	
	}
	if(fConversionCuts->PhotonIsSelectedMC(candidate,MCStack,kTRUE)){
		fGammaMCConvPt = candidate->Pt();
		fGammaMCConvTheta = candidate->Theta();
		TParticle* daughter1 = (TParticle *)MCStack->Particle(candidate->GetFirstDaughter()); 
		TParticle* daughter2 = (TParticle *)MCStack->Particle(candidate->GetLastDaughter()); 
      fMCConvCords(0) = (Float_t)daughter1->Vx();
      fMCConvCords(1) = (Float_t)daughter1->Vy();
      fMCConvCords(2) = (Float_t)daughter1->Vz();
      fMCConvCords(3) = (Float_t)daughter1->R();
      fMCConvCords(4) = (Float_t)daughter1->Phi();
      
      fMCConvDaughterProp(0) = (Float_t)daughter1->Pt();
      fMCConvDaughterProp(1) = (Float_t)daughter1->Theta();
      fMCConvDaughterProp(2) = (Float_t)daughter2->Pt();
      fMCConvDaughterProp(3) = (Float_t)daughter2->Theta();      
		
		if (fTreeMaterialConvGamma){
			fTreeMaterialConvGamma->Fill();
		}
	} // Converted MC Gamma		
}

///________________________________________________________________________
void AliAnalysisTaskMaterial::ProcessMCPhotons(){
	// Loop over all primary MC particle
	AliStack *ffMCStack = fMCEvent->Stack();
	for(Int_t i = 0; i < ffMCStack->GetNprimary(); i++) {
		TParticle* particle = (TParticle *)ffMCStack->Particle(i);
		if (!particle) continue;		
		if (particle->GetPdgCode() == 111 && particle->GetFirstDaughter() >= ffMCStack->GetNprimary()){
// 			cout << "Undecayed pi0 found with mother: " << particle->GetMother(0) << endl;
			for (Int_t j = 0; j < 2 ; j++){
 				FillMCTree(particle->GetDaughter(j));
			}
		} else {
			FillMCTree(i);
		}	
	}	
}

///________________________________________________________________________
void AliAnalysisTaskMaterial::ProcessPhotons(){

	// Fill Histograms for QA and MC
   for(Int_t firstGammaIndex=0;firstGammaIndex<fConversionGammas->GetEntriesFast();firstGammaIndex++){
      AliAODConversionPhoton *gamma=dynamic_cast<AliAODConversionPhoton*>(fConversionGammas->At(firstGammaIndex));
      if (gamma ==NULL) continue;
      if(!fConversionCuts->PhotonIsSelected(gamma,fESDEvent)) continue;

      fGammaPt = gamma->GetPhotonPt();
		fGammaTheta = gamma->GetPhotonTheta();
		fGammaChi2NDF = gamma->GetChi2perNDF();
      fRecCords(0) = (Float_t)gamma->GetConversionX();
      fRecCords(1) = (Float_t)gamma->GetConversionY();
      fRecCords(2) = (Float_t)gamma->GetConversionZ();
      fRecCords(3) = (Float_t)gamma->GetConversionRadius();
      fRecCords(4) = (Float_t)gamma->GetPhotonPhi();
      
		AliESDtrack * negTrack = fConversionCuts->GetESDTrack(fESDEvent, gamma->GetTrackLabelNegative());
      AliESDtrack * posTrack = fConversionCuts->GetESDTrack(fESDEvent, gamma->GetTrackLabelPositive());
      fDaughterProp(0) = (Float_t)posTrack->Pt();
      fDaughterProp(1) = (Float_t)posTrack->Theta();
      fDaughterProp(2) = (Float_t)negTrack->Pt();
      fDaughterProp(3) = (Float_t)negTrack->Theta();

		fKind = 9;	
		
		if(fMCEvent){
// 			cout << "generating MC stack"<< endl;
			AliStack *fMCStack = fMCEvent->Stack();
			if (!fMCStack) continue;
			TParticle *posDaughter = gamma->GetPositiveMCDaughter(fMCStack);
			TParticle *negDaughter = gamma->GetNegativeMCDaughter(fMCStack);
// 			cout << "generate Daughters: "<<posDaughter << "\t" << negDaughter << endl;
			
			if(posDaughter == NULL || negDaughter == NULL){ 
				fKind = 9; // garbage
// 				cout << "one of the daughters not available" << endl;
			} else if(posDaughter->GetMother(0) != negDaughter->GetMother(0) || (posDaughter->GetMother(0) == negDaughter->GetMother(0) && posDaughter->GetMother(0) ==-1)){ 
				// Not Same Mother == Combinatorial Bck
				fKind = 1;
// 				cout << "not the same mother" << endl;
				Int_t pdgCodePos; 
				if (posDaughter->GetPdgCode()) pdgCodePos = posDaughter->GetPdgCode(); else continue;
				Int_t pdgCodeNeg; 
				if (negDaughter->GetPdgCode()) pdgCodeNeg = negDaughter->GetPdgCode(); else continue;
// 				cout << "PDG codes daughters: " << pdgCodePos << "\t" << pdgCodeNeg << endl;
				if(TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==11)
					fKind = 10; //Electron Combinatorial
				if(TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==11 && (posDaughter->GetMother(0) == negDaughter->GetMother(0) && posDaughter->GetMother(0) ==-1))
					fKind = 15; //direct Electron Combinatorial
				
				if(TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==211)
					fKind = 11; //Pion Combinatorial
				if((TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==2212) ||
					(TMath::Abs(pdgCodePos)==2212 && TMath::Abs(pdgCodeNeg)==211))
					fKind = 12; //Pion, Proton Combinatorics
				if((TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==11) ||
					(TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==211))
					fKind = 13; //Pion, Electron Combinatorics
				if (TMath::Abs(pdgCodePos)==321 || TMath::Abs(pdgCodeNeg)==321)	
					fKind = 14; //Kaon combinatorics

			} else {
// 				cout << "same mother" << endl;
				Int_t pdgCodePos; 
				if (posDaughter->GetPdgCode()) pdgCodePos = posDaughter->GetPdgCode(); else continue;
				Int_t pdgCodeNeg; 
				if (negDaughter->GetPdgCode()) pdgCodeNeg = negDaughter->GetPdgCode(); else continue;
// 				cout << "PDG codes daughters: " << pdgCodePos << "\t" << pdgCodeNeg << endl;
				Int_t pdgCode; 
				if (gamma->GetMCParticle(fMCStack)->GetPdgCode()) pdgCode = gamma->GetMCParticle(fMCStack)->GetPdgCode(); else continue;
// 				cout << "PDG code: " << pdgCode << endl;
				if(TMath::Abs(pdgCodePos)!=11 || TMath::Abs(pdgCodeNeg)!=11)
					fKind = 2; // combinatorics from hadronic decays
				else if ( !(pdgCodeNeg==pdgCodePos)){
					TParticle *truePhotonCanditate = gamma->GetMCParticle(fMCStack);
					Int_t motherLabelPhoton = truePhotonCanditate->GetMother(0);
					if(pdgCode == 111) 
						fKind = 3; // pi0 Dalitz
					else if (pdgCode == 221) 
						fKind = 4; // eta Dalitz
					else if (!(negDaughter->GetUniqueID() != 5 || posDaughter->GetUniqueID() !=5)){
						if(pdgCode == 22 && motherLabelPhoton < fMCStack->GetNprimary()){
							fKind = 0; // primary photons
						} else if (pdgCode == 22){
							fKind = 5; //secondary photons
						}		
					} else 	fKind = 9; //garbage
				} else fKind = 9; //garbage
			}					
		}
		if (fTreeMaterialRec){
         fTreeMaterialRec->Fill();
      }
	}
}

//________________________________________________________________________
Int_t AliAnalysisTaskMaterial::CountTracks09(){
   Int_t fNumberOfESDTracks = 0;
   if(fInputEvent->IsA()==AliESDEvent::Class()){
   // Using standard function for setting Cuts
      Bool_t selectPrimaries=kTRUE;
      AliESDtrackCuts *EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
      EsdTrackCuts->SetMaxDCAToVertexZ(2);
      EsdTrackCuts->SetEtaRange(-0.9, 0.9);
      EsdTrackCuts->SetPtRange(0.15);
      
      for(Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++){
         AliESDtrack* curTrack = (AliESDtrack*) fInputEvent->GetTrack(iTracks);
         if(!curTrack) continue;
         if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
      }
      delete EsdTrackCuts;
      EsdTrackCuts=0x0;
   }
   else if(fInputEvent->IsA()==AliAODEvent::Class()){
      for(Int_t iTracks = 0; iTracks<fInputEvent->GetNumberOfTracks(); iTracks++){
         AliAODTrack* curTrack = (AliAODTrack*) fInputEvent->GetTrack(iTracks);
         if(!curTrack->IsPrimaryCandidate()) continue;
         if(abs(curTrack->Eta())>0.9) continue;
         if(curTrack->Pt()<0.15) continue;
         if(abs(curTrack->ZAtDCA())>2) continue;
         fNumberOfESDTracks++;
      }
   }

   return fNumberOfESDTracks;
}

Int_t AliAnalysisTaskMaterial::CountTracks0914(){

    Int_t fNumberOfESDTracks = 0;
    if(fInputEvent->IsA()==AliESDEvent::Class()){
      // Using standard function for setting Cuts
      AliESDtrackCuts *EsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
      EsdTrackCuts->SetMaxDCAToVertexZ(5);
      EsdTrackCuts->SetEtaRange(0.9, 1.4);
      EsdTrackCuts->SetPtRange(0.15);
      
      for(Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++){
         AliESDtrack* curTrack = (AliESDtrack*) fInputEvent->GetTrack(iTracks);
         if(!curTrack) continue;
         if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
      }
      EsdTrackCuts->SetEtaRange(-1.4, -0.9);
      for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
         AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
         if(!curTrack) continue;
         if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
      }
      delete EsdTrackCuts;
      EsdTrackCuts=0x0;
   }
   else if(fInputEvent->IsA()==AliAODEvent::Class()){
      for(Int_t iTracks = 0; iTracks<fInputEvent->GetNumberOfTracks(); iTracks++){
         AliAODTrack* curTrack = (AliAODTrack*) fInputEvent->GetTrack(iTracks);
         if(abs(curTrack->Eta())<0.9 || abs(curTrack->Eta())>1.4 ) continue;
         if(curTrack->Pt()<0.15) continue;
         if(abs(curTrack->ZAtDCA())>5) continue;
         fNumberOfESDTracks++;
      }
   }
   
   return fNumberOfESDTracks;
}

//________________________________________________________________________
void AliAnalysisTaskMaterial::Terminate(Option_t *)
{
//    if (fStreamMaterial){
//       fStreamMaterial->GetFile()->Write();
//    }
}
Commit	Line	Data
ca91a3e1	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Authors: Friederike Bock *
	5	* Version 1.0 *
	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	// QA Task for V0 Reader V1
	19	//---------------------------------------------
	20	////////////////////////////////////////////////
	21
	22	#include "AliAnalysisTaskMaterial.h"
	23	#include "TChain.h"
	24	#include "AliAnalysisManager.h"
	25	#include "TParticle.h"
	26	#include "TVectorF.h"
	27	#include "AliPIDResponse.h"
	28	#include "AliESDtrackCuts.h"
	29	#include "TFile.h"
	30
	31	class iostream;
	32
	33	using namespace std;
	34
	35	ClassImp(AliAnalysisTaskMaterial)
	36
72395bd9	37	AliAnalysisTaskMaterial::AliAnalysisTaskMaterial() : AliAnalysisTaskSE(),
	38	fV0Reader(NULL),
	39	fConversionGammas(NULL),
	40	fConversionCuts(NULL),
	41	fOutputList(NULL),
	42	fEventList(NULL),
	43	fRecGammaList(NULL),
	44	fAllMCGammaList(NULL),
	45	fAllMCConvGammaList(NULL),
	46	fTreeEvent(NULL),
	47	fTreeMaterialRec(NULL),
	48	fTreeMaterialAllGamma(NULL),
	49	fTreeMaterialConvGamma(NULL),
	50	fPrimVtxZ(0.),
	51	fNContrVtx(0),
	52	fNESDtracksEta09(0),
	53	fNESDtracksEta0914(0),
	54	fNESDtracksEta14(0),
	55	fGammaMCPt(0.),
	56	fGammaMCTheta(0.),
	57	fGammaMCConvPt(0.),
	58	fGammaMCConvTheta(0.),
31216736	59	fMCConvCords(5),
31216736	60	fMCConvDaughterProp(4),
72395bd9	61	fGammaPt(0.),
	62	fGammaTheta(0.),
	63	fGammaChi2NDF(0.),
31216736	64	fRecCords(5),
31216736	65	fDaughterProp(4),
72395bd9	66	fKind(0),
	67	fIsHeavyIon(kFALSE),
	68	fIsMC(kFALSE),
	69	fESDEvent(NULL),
	70	fMCEvent(NULL)
	71	{
	72
	73
	74	}
	75
	76
ca91a3e1	77	//________________________________________________________________________
ca91a3e1	78	AliAnalysisTaskMaterial::AliAnalysisTaskMaterial(const char *name) : AliAnalysisTaskSE(name),
a280ac15	79	fV0Reader(NULL),
ca91a3e1	80	fConversionGammas(NULL),
ca91a3e1	81	fConversionCuts(NULL),
ca91a3e1	82	fOutputList(NULL),
72395bd9	83	fEventList(NULL),
	84	fRecGammaList(NULL),
	85	fAllMCGammaList(NULL),
	86	fAllMCConvGammaList(NULL),
	87	fTreeEvent(NULL),
	88	fTreeMaterialRec(NULL),
	89	fTreeMaterialAllGamma(NULL),
	90	fTreeMaterialConvGamma(NULL),
	91	fPrimVtxZ(0.),
	92	fNContrVtx(0),
	93	fNESDtracksEta09(0),
	94	fNESDtracksEta0914(0),
	95	fNESDtracksEta14(0),
	96	fGammaMCPt(0.),
	97	fGammaMCTheta(0.),
	98	fGammaMCConvPt(0.),
	99	fGammaMCConvTheta(0.),
31216736	100	fMCConvCords(5),
31216736	101	fMCConvDaughterProp(4),
72395bd9	102	fGammaPt(0.),
	103	fGammaTheta(0.),
	104	fGammaChi2NDF(0.),
31216736	105	fRecCords(5),
31216736	106	fDaughterProp(4),
72395bd9	107	fKind(0),
	108	fIsHeavyIon(kFALSE),
	109	fIsMC(kFALSE),
ca91a3e1	110	fESDEvent(NULL),
	111	fMCEvent(NULL)
	112	{
	113	// Default constructor
	114
72395bd9	115
ca91a3e1	116	DefineInput(0, TChain::Class());
	117	DefineOutput(1, TList::Class());
	118	}
	119
	120	//________________________________________________________________________
	121	AliAnalysisTaskMaterial::~AliAnalysisTaskMaterial()
	122	{
	123	// default deconstructor
ca91a3e1	124	}
	125	//________________________________________________________________________
	126	void AliAnalysisTaskMaterial::UserCreateOutputObjects()
	127	{
	128	// Create User Output Objects
	129
	130	if(fOutputList != NULL){
	131	delete fOutputList;
	132	fOutputList = NULL;
	133	}
	134	if(fOutputList == NULL){
	135	fOutputList = new TList();
	136	fOutputList->SetOwner(kTRUE);
	137	}
72395bd9	138
	139	fEventList = new TList();
	140	fEventList->SetName("EventList");
	141	fEventList->SetOwner(kTRUE);
	142	fOutputList->Add(fEventList);
	143
	144	fTreeEvent = new TTree("Event","Event");
	145	fTreeEvent->Branch("primVtxZ",&fPrimVtxZ,"fPrimVtxZ/F");
	146	fTreeEvent->Branch("nContrVtx",&fNContrVtx,"fNContrVtx/I");
	147	fTreeEvent->Branch("nGoodTracksEta09",&fNESDtracksEta09,"fNESDtracksEta09/I");
	148	fTreeEvent->Branch("nGoodTracksEta0914",&fNESDtracksEta0914,"fNESDtracksEta0914/I");
	149	fTreeEvent->Branch("nGoodTracksEta14",&fNESDtracksEta14,"fNESDtracksEta14/I");
	150	fEventList->Add(fTreeEvent);
	151
	152	fRecGammaList= new TList();
	153	fRecGammaList->SetName("RecGammaList");
	154	fRecGammaList->SetOwner(kTRUE);
	155	fOutputList->Add(fRecGammaList);
	156
72395bd9	157
72395bd9	158	fTreeMaterialRec = new TTree("ConvPointRec","ConvPointRec");
31216736	159	fTreeMaterialRec->Branch("recCords",&fRecCords);
31216736	160	fTreeMaterialRec->Branch("daughterProp",&fDaughterProp);
72395bd9	161	fTreeMaterialRec->Branch("pt",&fGammaPt,"fGammaPt/F");
	162	fTreeMaterialRec->Branch("theta",&fGammaTheta,"fGammaTheta/F");
	163	fTreeMaterialRec->Branch("chi2ndf",&fGammaChi2NDF,"fGammaChi2NDF/F");
72395bd9	164	if (fIsMC) {
	165	fTreeMaterialRec->Branch("kind",&fKind,"fKind/b");
	166	}
	167	fRecGammaList->Add(fTreeMaterialRec);
	168
	169	if (fIsMC) {
	170	fAllMCGammaList = new TList();
	171	fAllMCGammaList->SetName("AllMCGammaList");
	172	fAllMCGammaList->SetOwner(kTRUE);
	173	fOutputList->Add(fAllMCGammaList);
	174
	175	fTreeMaterialAllGamma = new TTree("AllGamma","AllGamma");
	176	fTreeMaterialAllGamma->Branch("pt",&fGammaMCPt,"fGammaMCPt/F");
	177	fTreeMaterialAllGamma->Branch("theta",&fGammaMCTheta,"fGammaMCTheta/I");
	178	fAllMCGammaList->Add(fTreeMaterialAllGamma);
	179
	180	fAllMCConvGammaList = new TList();
	181	fAllMCConvGammaList->SetName("AllMCGammaConvList");
	182	fAllMCConvGammaList->SetOwner(kTRUE);
	183	fOutputList->Add(fAllMCConvGammaList);
	184
	185	// fMCConvCords = new Float_t[5];
	186	// fMCConvDaughterProp = new Float_t[4];
ca91a3e1	187
72395bd9	188
72395bd9	189	fTreeMaterialConvGamma = new TTree("ConvGammaMC","ConvGammaMC");
31216736	190	fTreeMaterialConvGamma->Branch("Cords",&fMCConvCords);
31216736	191	fTreeMaterialConvGamma->Branch("daughterProp",&fMCConvDaughterProp);
72395bd9	192	fTreeMaterialConvGamma->Branch("Pt",&fGammaMCConvPt,"fGammaMCConvPt/F");
72395bd9	193	fTreeMaterialConvGamma->Branch("Theta",&fGammaMCConvTheta,"fGammaMCConvTheta/F");
72395bd9	194	fAllMCConvGammaList->Add(fTreeMaterialConvGamma);
	195	}
	196
ca91a3e1	197	// V0 Reader Cuts
ca91a3e1	198	TString cutnumber = fConversionCuts->GetCutNumber();
ca91a3e1	199	PostData(1, fOutputList);
72395bd9	200
ca91a3e1	201	}
	202
	203	//________________________________________________________________________
	204	void AliAnalysisTaskMaterial::UserExec(Option_t *){
	205
	206	fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1");
	207
	208	Int_t eventQuality = ((AliConversionCuts*)fV0Reader->GetConversionCuts())->GetEventQuality();
e5b6e8a6	209	if(eventQuality != 0){// Event Not Accepted
ca91a3e1	210	return;
	211	}
	212	fESDEvent = (AliESDEvent*) InputEvent();
0a2b2b4b	213	if (fESDEvent==NULL) return;
ca91a3e1	214	if(fIsHeavyIon && !fConversionCuts->IsCentralitySelected(fESDEvent)) return;
72395bd9	215	fNESDtracksEta09 = CountTracks09(); // Estimate Event Multiplicity
	216	fNESDtracksEta0914 = CountTracks0914(); // Estimate Event Multiplicity
	217	fNESDtracksEta14 = fNESDtracksEta09 + fNESDtracksEta0914;
ca91a3e1	218	if(fESDEvent){
ca91a3e1	219	if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()>0) {
72395bd9	220	fNContrVtx = fESDEvent->GetPrimaryVertexTracks()->GetNContributors();
ca91a3e1	221	} else if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()<1) {
ca91a3e1	222	if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
72395bd9	223	fNContrVtx = fESDEvent->GetPrimaryVertexSPD()->GetNContributors();
ca91a3e1	224	} else if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()<1) {
72395bd9	225	fNContrVtx = 0;
ca91a3e1	226	}
	227	}
	228	}
72395bd9	229	fPrimVtxZ = fESDEvent->GetPrimaryVertex()->GetZ();
ca91a3e1	230
72395bd9	231	if (fTreeEvent){
72395bd9	232	fTreeEvent->Fill();
ca91a3e1	233	}
	234
	235	fConversionGammas=fV0Reader->GetReconstructedGammas();
	236	if(MCEvent()){
	237	fMCEvent = MCEvent();
	238	}
	239	ProcessPhotons();
	240	if(MCEvent()){
	241	ProcessMCPhotons();
	242	}
	243	PostData(1, fOutputList);
	244	}
	245
e5b6e8a6	246	///________________________________________________________________________
	247	void AliAnalysisTaskMaterial::FillMCTree(Int_t stackPos){
	248	AliStack *MCStack = fMCEvent->Stack();
	249	TParticle* candidate = (TParticle *)MCStack->Particle(stackPos);
	250	if(fConversionCuts->PhotonIsSelectedMC(candidate,MCStack,kFALSE)){
72395bd9	251	fGammaMCPt = candidate->Pt();
	252	fGammaMCTheta = candidate->Theta();
	253	if (fTreeMaterialAllGamma){
	254	fTreeMaterialAllGamma->Fill();
e5b6e8a6	255	}
	256	}
	257	if(fConversionCuts->PhotonIsSelectedMC(candidate,MCStack,kTRUE)){
72395bd9	258	fGammaMCConvPt = candidate->Pt();
72395bd9	259	fGammaMCConvTheta = candidate->Theta();
e5b6e8a6	260	TParticle* daughter1 = (TParticle *)MCStack->Particle(candidate->GetFirstDaughter());
e5b6e8a6	261	TParticle* daughter2 = (TParticle *)MCStack->Particle(candidate->GetLastDaughter());
31216736	262	fMCConvCords(0) = (Float_t)daughter1->Vx();
	263	fMCConvCords(1) = (Float_t)daughter1->Vy();
	264	fMCConvCords(2) = (Float_t)daughter1->Vz();
	265	fMCConvCords(3) = (Float_t)daughter1->R();
	266	fMCConvCords(4) = (Float_t)daughter1->Phi();
72395bd9	267
31216736	268	fMCConvDaughterProp(0) = (Float_t)daughter1->Pt();
	269	fMCConvDaughterProp(1) = (Float_t)daughter1->Theta();
	270	fMCConvDaughterProp(2) = (Float_t)daughter2->Pt();
	271	fMCConvDaughterProp(3) = (Float_t)daughter2->Theta();
72395bd9	272
	273	if (fTreeMaterialConvGamma){
	274	fTreeMaterialConvGamma->Fill();
e5b6e8a6	275	}
	276	} // Converted MC Gamma
	277	}
ca91a3e1	278
	279	///________________________________________________________________________
	280	void AliAnalysisTaskMaterial::ProcessMCPhotons(){
	281	// Loop over all primary MC particle
	282	AliStack *ffMCStack = fMCEvent->Stack();
	283	for(Int_t i = 0; i < ffMCStack->GetNprimary(); i++) {
	284	TParticle* particle = (TParticle *)ffMCStack->Particle(i);
	285	if (!particle) continue;
e5b6e8a6	286	if (particle->GetPdgCode() == 111 && particle->GetFirstDaughter() >= ffMCStack->GetNprimary()){
	287	// cout << "Undecayed pi0 found with mother: " << particle->GetMother(0) << endl;
	288	for (Int_t j = 0; j < 2 ; j++){
	289	FillMCTree(particle->GetDaughter(j));
ca91a3e1	290	}
e5b6e8a6	291	} else {
	292	FillMCTree(i);
	293	}
ca91a3e1	294	}
	295	}
	296
	297	///________________________________________________________________________
	298	void AliAnalysisTaskMaterial::ProcessPhotons(){
	299
	300	// Fill Histograms for QA and MC
	301	for(Int_t firstGammaIndex=0;firstGammaIndex<fConversionGammas->GetEntriesFast();firstGammaIndex++){
	302	AliAODConversionPhoton gamma=dynamic_cast<AliAODConversionPhoton>(fConversionGammas->At(firstGammaIndex));
0a2b2b4b	303	if (gamma ==NULL) continue;
ca91a3e1	304	if(!fConversionCuts->PhotonIsSelected(gamma,fESDEvent)) continue;
72395bd9	305
	306	fGammaPt = gamma->GetPhotonPt();
	307	fGammaTheta = gamma->GetPhotonTheta();
	308	fGammaChi2NDF = gamma->GetChi2perNDF();
31216736	309	fRecCords(0) = (Float_t)gamma->GetConversionX();
	310	fRecCords(1) = (Float_t)gamma->GetConversionY();
	311	fRecCords(2) = (Float_t)gamma->GetConversionZ();
	312	fRecCords(3) = (Float_t)gamma->GetConversionRadius();
	313	fRecCords(4) = (Float_t)gamma->GetPhotonPhi();
72395bd9	314
e5b6e8a6	315	AliESDtrack * negTrack = fConversionCuts->GetESDTrack(fESDEvent, gamma->GetTrackLabelNegative());
e5b6e8a6	316	AliESDtrack * posTrack = fConversionCuts->GetESDTrack(fESDEvent, gamma->GetTrackLabelPositive());
31216736	317	fDaughterProp(0) = (Float_t)posTrack->Pt();
	318	fDaughterProp(1) = (Float_t)posTrack->Theta();
	319	fDaughterProp(2) = (Float_t)negTrack->Pt();
	320	fDaughterProp(3) = (Float_t)negTrack->Theta();
72395bd9	321
	322	fKind = 9;
	323
ca91a3e1	324	if(fMCEvent){
	325	// cout << "generating MC stack"<< endl;
	326	AliStack *fMCStack = fMCEvent->Stack();
	327	if (!fMCStack) continue;
	328	TParticle *posDaughter = gamma->GetPositiveMCDaughter(fMCStack);
	329	TParticle *negDaughter = gamma->GetNegativeMCDaughter(fMCStack);
	330	// cout << "generate Daughters: "<<posDaughter << "\t" << negDaughter << endl;
	331
	332	if(posDaughter == NULL \|\| negDaughter == NULL){
72395bd9	333	fKind = 9; // garbage
ca91a3e1	334	// cout << "one of the daughters not available" << endl;
	335	} else if(posDaughter->GetMother(0) != negDaughter->GetMother(0) \|\| (posDaughter->GetMother(0) == negDaughter->GetMother(0) && posDaughter->GetMother(0) ==-1)){
	336	// Not Same Mother == Combinatorial Bck
72395bd9	337	fKind = 1;
ca91a3e1	338	// cout << "not the same mother" << endl;
	339	Int_t pdgCodePos;
	340	if (posDaughter->GetPdgCode()) pdgCodePos = posDaughter->GetPdgCode(); else continue;
	341	Int_t pdgCodeNeg;
	342	if (negDaughter->GetPdgCode()) pdgCodeNeg = negDaughter->GetPdgCode(); else continue;
	343	// cout << "PDG codes daughters: " << pdgCodePos << "\t" << pdgCodeNeg << endl;
	344	if(TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==11)
72395bd9	345	fKind = 10; //Electron Combinatorial
ca91a3e1	346	if(TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==11 && (posDaughter->GetMother(0) == negDaughter->GetMother(0) && posDaughter->GetMother(0) ==-1))
72395bd9	347	fKind = 15; //direct Electron Combinatorial
ca91a3e1	348
ca91a3e1	349	if(TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==211)
72395bd9	350	fKind = 11; //Pion Combinatorial
ca91a3e1	351	if((TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==2212) \|\|
ca91a3e1	352	(TMath::Abs(pdgCodePos)==2212 && TMath::Abs(pdgCodeNeg)==211))
72395bd9	353	fKind = 12; //Pion, Proton Combinatorics
ca91a3e1	354	if((TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==11) \|\|
ca91a3e1	355	(TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==211))
72395bd9	356	fKind = 13; //Pion, Electron Combinatorics
ca91a3e1	357	if (TMath::Abs(pdgCodePos)==321 \|\| TMath::Abs(pdgCodeNeg)==321)
72395bd9	358	fKind = 14; //Kaon combinatorics
ca91a3e1	359
	360	} else {
	361	// cout << "same mother" << endl;
	362	Int_t pdgCodePos;
	363	if (posDaughter->GetPdgCode()) pdgCodePos = posDaughter->GetPdgCode(); else continue;
	364	Int_t pdgCodeNeg;
	365	if (negDaughter->GetPdgCode()) pdgCodeNeg = negDaughter->GetPdgCode(); else continue;
	366	// cout << "PDG codes daughters: " << pdgCodePos << "\t" << pdgCodeNeg << endl;
	367	Int_t pdgCode;
	368	if (gamma->GetMCParticle(fMCStack)->GetPdgCode()) pdgCode = gamma->GetMCParticle(fMCStack)->GetPdgCode(); else continue;
	369	// cout << "PDG code: " << pdgCode << endl;
	370	if(TMath::Abs(pdgCodePos)!=11 \|\| TMath::Abs(pdgCodeNeg)!=11)
72395bd9	371	fKind = 2; // combinatorics from hadronic decays
ca91a3e1	372	else if ( !(pdgCodeNeg==pdgCodePos)){
ca91a3e1	373	TParticle *truePhotonCanditate = gamma->GetMCParticle(fMCStack);
e5b6e8a6	374	Int_t motherLabelPhoton = truePhotonCanditate->GetMother(0);
ca91a3e1	375	if(pdgCode == 111)
72395bd9	376	fKind = 3; // pi0 Dalitz
ca91a3e1	377	else if (pdgCode == 221)
72395bd9	378	fKind = 4; // eta Dalitz
ca91a3e1	379	else if (!(negDaughter->GetUniqueID() != 5 \|\| posDaughter->GetUniqueID() !=5)){
e5b6e8a6	380	if(pdgCode == 22 && motherLabelPhoton < fMCStack->GetNprimary()){
72395bd9	381	fKind = 0; // primary photons
ca91a3e1	382	} else if (pdgCode == 22){
72395bd9	383	fKind = 5; //secondary photons
	384	}
	385	} else fKind = 9; //garbage
	386	} else fKind = 9; //garbage
	387	}
ca91a3e1	388	}
72395bd9	389	if (fTreeMaterialRec){
	390	fTreeMaterialRec->Fill();
	391	}
ca91a3e1	392	}
	393	}
	394
	395	//________________________________________________________________________
ae947965	396	Int_t AliAnalysisTaskMaterial::CountTracks09(){
e5b6e8a6	397	Int_t fNumberOfESDTracks = 0;
ae947965	398	if(fInputEvent->IsA()==AliESDEvent::Class()){
	399	// Using standard function for setting Cuts
	400	Bool_t selectPrimaries=kTRUE;
	401	AliESDtrackCuts *EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
	402	EsdTrackCuts->SetMaxDCAToVertexZ(2);
	403	EsdTrackCuts->SetEtaRange(-0.9, 0.9);
	404	EsdTrackCuts->SetPtRange(0.15);
	405
	406	for(Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++){
	407	AliESDtrack* curTrack = (AliESDtrack*) fInputEvent->GetTrack(iTracks);
	408	if(!curTrack) continue;
ae947965	409	if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
	410	}
	411	delete EsdTrackCuts;
	412	EsdTrackCuts=0x0;
	413	}
	414	else if(fInputEvent->IsA()==AliAODEvent::Class()){
	415	for(Int_t iTracks = 0; iTracks<fInputEvent->GetNumberOfTracks(); iTracks++){
	416	AliAODTrack* curTrack = (AliAODTrack*) fInputEvent->GetTrack(iTracks);
	417	if(!curTrack->IsPrimaryCandidate()) continue;
	418	if(abs(curTrack->Eta())>0.9) continue;
	419	if(curTrack->Pt()<0.15) continue;
	420	if(abs(curTrack->ZAtDCA())>2) continue;
	421	fNumberOfESDTracks++;
	422	}
e5b6e8a6	423	}
e5b6e8a6	424
	425	return fNumberOfESDTracks;
	426	}
	427
ae947965	428	Int_t AliAnalysisTaskMaterial::CountTracks0914(){
e5b6e8a6	429
ae947965	430	Int_t fNumberOfESDTracks = 0;
	431	if(fInputEvent->IsA()==AliESDEvent::Class()){
	432	// Using standard function for setting Cuts
	433	AliESDtrackCuts *EsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
	434	EsdTrackCuts->SetMaxDCAToVertexZ(5);
	435	EsdTrackCuts->SetEtaRange(0.9, 1.4);
	436	EsdTrackCuts->SetPtRange(0.15);
	437
	438	for(Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++){
	439	AliESDtrack* curTrack = (AliESDtrack*) fInputEvent->GetTrack(iTracks);
	440	if(!curTrack) continue;
ae947965	441	if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
	442	}
	443	EsdTrackCuts->SetEtaRange(-1.4, -0.9);
	444	for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
	445	AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
	446	if(!curTrack) continue;
	447	if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
	448	}
	449	delete EsdTrackCuts;
	450	EsdTrackCuts=0x0;
e5b6e8a6	451	}
ae947965	452	else if(fInputEvent->IsA()==AliAODEvent::Class()){
	453	for(Int_t iTracks = 0; iTracks<fInputEvent->GetNumberOfTracks(); iTracks++){
	454	AliAODTrack* curTrack = (AliAODTrack*) fInputEvent->GetTrack(iTracks);
ae947965	455	if(abs(curTrack->Eta())<0.9 \|\| abs(curTrack->Eta())>1.4 ) continue;
	456	if(curTrack->Pt()<0.15) continue;
	457	if(abs(curTrack->ZAtDCA())>5) continue;
	458	fNumberOfESDTracks++;
	459	}
e5b6e8a6	460	}
ae947965	461
ca91a3e1	462	return fNumberOfESDTracks;
	463	}
	464
ca91a3e1	465	//________________________________________________________________________
	466	void AliAnalysisTaskMaterial::Terminate(Option_t *)
	467	{
0a2b2b4b	468	// if (fStreamMaterial){
	469	// fStreamMaterial->GetFile()->Write();
	470	// }
ca91a3e1	471	}