[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(const char *name) : AliAnalysisTaskSE(name),
   fV0Reader(NULL),
   fConversionGammas(NULL),
   fConversionCuts(NULL),
   fStreamMaterial(NULL),
   fStreamResolution(NULL),
   fIsHeavyIon(kFALSE),
   fOutputList(NULL),
   fESDEvent(NULL),
   fMCEvent(NULL)
{
   // Default constructor

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

//________________________________________________________________________
AliAnalysisTaskMaterial::~AliAnalysisTaskMaterial()
{
   // default deconstructor
   if(fStreamMaterial){
      delete fStreamMaterial;
      fStreamMaterial = 0x0;
   }
   if(fStreamResolution){
      delete fStreamResolution;
      fStreamResolution = 0x0;
   }
}
//________________________________________________________________________
void AliAnalysisTaskMaterial::UserCreateOutputObjects()
{
   // Create User Output Objects

   if(fOutputList != NULL){
      delete fOutputList;
      fOutputList = NULL;
   }
   if(fOutputList == NULL){
      fOutputList = new TList();
      fOutputList->SetOwner(kTRUE);
   }

   // V0 Reader Cuts
   TString cutnumber = fConversionCuts->GetCutNumber();

   fStreamMaterial = new TTreeSRedirector(Form("GammaConvV1_Material_%s.root",cutnumber.Data()),"recreate");
	 fStreamResolution = new TTreeSRedirector(Form("GammaConvV1_Resolution_%s.root",cutnumber.Data()),"recreate");
   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;
	Int_t nESDtracksEta09 = CountTracks09(); // Estimate Event Multiplicity
	Int_t nESDtracksEta0914 = CountTracks0914(); // Estimate Event Multiplicity
	Int_t nESDtracksEta14; // Estimate Event Multiplicity
  nESDtracksEta14= nESDtracksEta09 + nESDtracksEta0914;
	Int_t nContrVtx;
	if(fESDEvent){
		if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()>0) {
			nContrVtx = fESDEvent->GetPrimaryVertexTracks()->GetNContributors();
		} else if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()<1) {
			if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
				nContrVtx = fESDEvent->GetPrimaryVertexSPD()->GetNContributors();

			} else if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()<1) {
				nContrVtx = 0;
			}
		}
   }
	Float_t primVtxZ = fESDEvent->GetPrimaryVertex()->GetZ();
	
	if (fStreamMaterial){
		(*fStreamMaterial)<<"Event"
						<< "primVtxZ=" << primVtxZ
						<< "nContrVtx=" << nContrVtx
						<< "nGoodTracksEta09=" << nESDtracksEta09
						<< "nGoodTracksEta0914=" << nESDtracksEta0914
						<< "nGoodTracksEta14=" << nESDtracksEta14
						<< "\n";
	}
	
   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)){
		Float_t gammaPt = candidate->Pt();
		Float_t gammaTheta = candidate->Theta();
		if (fStreamMaterial){
			(*fStreamMaterial)<<"AllGamma"
							<< "pt=" << gammaPt
// 							<< "p=" << gammaP
							<< "theta=" << gammaTheta
							<< "\n";
		}	
	}
	if(fConversionCuts->PhotonIsSelectedMC(candidate,MCStack,kTRUE)){
		Float_t gammaPt = candidate->Pt();
		Float_t gammaTheta = candidate->Theta();
		TParticle* daughter1 = (TParticle *)MCStack->Particle(candidate->GetFirstDaughter()); 
		TParticle* daughter2 = (TParticle *)MCStack->Particle(candidate->GetLastDaughter()); 
		TVectorF coord(5);	
		coord(0) = daughter1->Vx();
		coord(1) = daughter1->Vy();
		coord(2) = daughter1->Vz();
		coord(3) = daughter1->R();
		coord(4) = candidate->Phi();
		TVectorF daughterProp(4);	
		if (daughter1->	GetPdgCode() < 0){
			daughterProp(0) = daughter2->Pt();
			daughterProp(1) = daughter2->Theta();
			daughterProp(2) = daughter1->Pt();
			daughterProp(3) = daughter1->Theta();
		} else {
			daughterProp(0) = daughter1->Pt();
			daughterProp(1) = daughter1->Theta();
			daughterProp(2) = daughter2->Pt();
			daughterProp(3) = daughter2->Theta();
		}
		if (fStreamMaterial){
			(*fStreamMaterial)<<"ConvGammaMC"
							<< "pt=" << gammaPt
// 							<< "p=" << gammaP
 							<< "theta=" << gammaTheta
							<< "coord.=" << &coord
							<< "daughterProp.=" << &daughterProp
							<< "\n";
		}
	} // 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;
// 		cout << "i=  " <<firstGammaIndex << " of "<< fConversionGammas->GetEntriesFast() << endl;
      Float_t gammaPt = gamma->GetPhotonPt();
// 		Float_t gammaP = gamma->GetPhotonP();
		Float_t gammaTheta = gamma->GetPhotonTheta();
		Float_t gammaEta = gamma->GetPhotonEta();
		Float_t gammaChi2NDF = gamma->GetChi2perNDF();
// 		Float_t gammaX = gamma->GetConversionX();
// 		Float_t gammaY = gamma->GetConversionY();
 		Float_t gammaZ = gamma->GetConversionZ();
      Float_t gammaR = gamma->GetConversionRadius();
 		Float_t gammaPhi = gamma->GetPhotonPhi();
		
		TVectorF coord(5);	
		coord(0) = gamma->GetConversionX();
		coord(1) = gamma->GetConversionY();
		coord(2) = gamma->GetConversionZ();
		coord(3) = gamma->GetConversionRadius();
		coord(4) = gamma->GetPhotonPhi();
		TVectorF daughterProp(4);	
		AliESDtrack * negTrack = fConversionCuts->GetESDTrack(fESDEvent, gamma->GetTrackLabelNegative());
      AliESDtrack * posTrack = fConversionCuts->GetESDTrack(fESDEvent, gamma->GetTrackLabelPositive());
		daughterProp(0) = posTrack->Pt();
      daughterProp(1) = posTrack->Theta();
      daughterProp(2) = negTrack->Pt();
      daughterProp(3) = negTrack->Theta();
				
		UInt_t kind = 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){ 
				kind = 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
				kind = 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)
					kind = 10; //Electron Combinatorial
				if(TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==11 && (posDaughter->GetMother(0) == negDaughter->GetMother(0) && posDaughter->GetMother(0) ==-1))
					kind = 15; //direct Electron Combinatorial
				
				if(TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==211)
					kind = 11; //Pion Combinatorial
				if((TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==2212) ||
					(TMath::Abs(pdgCodePos)==2212 && TMath::Abs(pdgCodeNeg)==211))
					kind = 12; //Pion, Proton Combinatorics
				if((TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==11) ||
					(TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==211))
					kind = 13; //Pion, Electron Combinatorics
				if (TMath::Abs(pdgCodePos)==321 || TMath::Abs(pdgCodeNeg)==321)	
					kind = 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)
					kind = 2; // combinatorics from hadronic decays
				else if ( !(pdgCodeNeg==pdgCodePos)){
					TParticle *truePhotonCanditate = gamma->GetMCParticle(fMCStack);
					Int_t motherLabelPhoton = truePhotonCanditate->GetMother(0);
					if(pdgCode == 111) 
						kind = 3; // pi0 Dalitz
					else if (pdgCode == 221) 
						kind = 4; // eta Dalitz
					else if (!(negDaughter->GetUniqueID() != 5 || posDaughter->GetUniqueID() !=5)){
						if(pdgCode == 22 && motherLabelPhoton < fMCStack->GetNprimary()){
							kind = 0; // primary photons
						} else if (pdgCode == 22){
							kind = 5; //secondary photons
						}
						if(pdgCode == 22){
							Float_t mcPt   = truePhotonCanditate->Pt();
							Float_t mcR = gamma->GetNegativeMCDaughter(fMCStack)->R();
							Float_t mcZ = gamma->GetNegativeMCDaughter(fMCStack)->Vz();
							Float_t mcPhi = gamma->GetNegativeMCDaughter(fMCStack)->Phi();
							Float_t mcEta = gamma->GetNegativeMCDaughter(fMCStack)->Eta();
							if (fStreamResolution){
								(*fStreamResolution)<<"Resolution"
								<< "ESDpt=" << gammaPt
								<< "ESDphi=" << gammaPhi
								<< "ESDeta=" << gammaEta
								<< "ESDr="<< gammaR
								<< "ESDz="<< gammaZ
								<< "MCpt=" << mcPt
								<< "MCphi=" << mcPhi
								<< "MCeta=" << mcEta
								<< "MCr="<< mcR
								<< "MCz="<< mcZ
								<< "chi2ndf=" << gammaChi2NDF
								<< "\n";
							}
						}		
					} else 	kind = 9; //garbage
				} else kind = 9; //garbage
			}
// 			cout << gammaPt << "\t" << gammaP<< "\t" << gammaEta<<  "\t" <<gammaChi2NDF << "\t" << gammaX<<  "\t" <<gammaY << "\t" << gammaZ<< "\t" << gammaR<<  "\t" <<  gammaPhi<< "\t" <<kind << endl;
		
			if (fStreamMaterial){
				(*fStreamMaterial)<<"ConvPointRec"
								<< "pt=" << gammaPt
								<< "theta=" << gammaTheta
								<< "chi2ndf=" << gammaChi2NDF
								<< "kind=" << kind
								<< "coord.=" << &coord
								<< "daughterProp.=" << &daughterProp
								<< "\n";
			}		
		} else {
				if (fStreamMaterial){
				(*fStreamMaterial)<<"ConvPointRec"
								<< "pt=" << gammaPt
								<< "theta=" << gammaTheta
								<< "chi2ndf=" << gammaChi2NDF
								<< "coord.=" << &coord
								<< "daughterProp.=" << &daughterProp
								<< "\n";
			}
		}
	}
}

//________________________________________________________________________
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(fMCEvent && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){
         //    if(!((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(abs(curTrack->GetLabel()), fMCStack)) 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(){

   // Using standard function for setting Cuts ; We use TPCOnlyTracks for outer eta region
   //Bool_t selectPrimaries=kTRUE;
	 //   EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
    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(fMCEvent && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){
         //    if(!((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(abs(curTrack->GetLabel()), fMCStack)) 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(!curTrack->IsPrimaryCandidate()) continue;
         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();
//    }
//    if (fStreamResolution){
//       fStreamResolution->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
	37	//________________________________________________________________________
	38	AliAnalysisTaskMaterial::AliAnalysisTaskMaterial(const char *name) : AliAnalysisTaskSE(name),
a280ac15	39	fV0Reader(NULL),
ca91a3e1	40	fConversionGammas(NULL),
	41	fConversionCuts(NULL),
	42	fStreamMaterial(NULL),
	43	fStreamResolution(NULL),
	44	fIsHeavyIon(kFALSE),
	45	fOutputList(NULL),
	46	fESDEvent(NULL),
	47	fMCEvent(NULL)
	48	{
	49	// Default constructor
	50
	51	DefineInput(0, TChain::Class());
	52	DefineOutput(1, TList::Class());
	53	}
	54
	55	//________________________________________________________________________
	56	AliAnalysisTaskMaterial::~AliAnalysisTaskMaterial()
	57	{
	58	// default deconstructor
	59	if(fStreamMaterial){
	60	delete fStreamMaterial;
	61	fStreamMaterial = 0x0;
	62	}
	63	if(fStreamResolution){
	64	delete fStreamResolution;
	65	fStreamResolution = 0x0;
	66	}
	67	}
	68	//________________________________________________________________________
	69	void AliAnalysisTaskMaterial::UserCreateOutputObjects()
	70	{
	71	// Create User Output Objects
	72
	73	if(fOutputList != NULL){
	74	delete fOutputList;
	75	fOutputList = NULL;
	76	}
	77	if(fOutputList == NULL){
	78	fOutputList = new TList();
	79	fOutputList->SetOwner(kTRUE);
	80	}
	81
	82	// V0 Reader Cuts
	83	TString cutnumber = fConversionCuts->GetCutNumber();
	84
a280ac15	85	fStreamMaterial = new TTreeSRedirector(Form("GammaConvV1_Material_%s.root",cutnumber.Data()),"recreate");
a280ac15	86	fStreamResolution = new TTreeSRedirector(Form("GammaConvV1_Resolution_%s.root",cutnumber.Data()),"recreate");
ca91a3e1	87	PostData(1, fOutputList);
	88	}
	89
	90	//________________________________________________________________________
	91	void AliAnalysisTaskMaterial::UserExec(Option_t *){
	92
	93	fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1");
	94
	95	Int_t eventQuality = ((AliConversionCuts*)fV0Reader->GetConversionCuts())->GetEventQuality();
e5b6e8a6	96	if(eventQuality != 0){// Event Not Accepted
ca91a3e1	97	return;
	98	}
	99	fESDEvent = (AliESDEvent*) InputEvent();
0a2b2b4b	100	if (fESDEvent==NULL) return;
ca91a3e1	101	if(fIsHeavyIon && !fConversionCuts->IsCentralitySelected(fESDEvent)) return;
ae947965	102	Int_t nESDtracksEta09 = CountTracks09(); // Estimate Event Multiplicity
ae947965	103	Int_t nESDtracksEta0914 = CountTracks0914(); // Estimate Event Multiplicity
e5b6e8a6	104	Int_t nESDtracksEta14; // Estimate Event Multiplicity
e5b6e8a6	105	nESDtracksEta14= nESDtracksEta09 + nESDtracksEta0914;
ca91a3e1	106	Int_t nContrVtx;
	107	if(fESDEvent){
	108	if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()>0) {
	109	nContrVtx = fESDEvent->GetPrimaryVertexTracks()->GetNContributors();
	110	} else if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()<1) {
	111	if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
	112	nContrVtx = fESDEvent->GetPrimaryVertexSPD()->GetNContributors();
	113
	114	} else if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()<1) {
	115	nContrVtx = 0;
	116	}
	117	}
	118	}
e5b6e8a6	119	Float_t primVtxZ = fESDEvent->GetPrimaryVertex()->GetZ();
ca91a3e1	120
	121	if (fStreamMaterial){
	122	(*fStreamMaterial)<<"Event"
	123	<< "primVtxZ=" << primVtxZ
	124	<< "nContrVtx=" << nContrVtx
e5b6e8a6	125	<< "nGoodTracksEta09=" << nESDtracksEta09
	126	<< "nGoodTracksEta0914=" << nESDtracksEta0914
	127	<< "nGoodTracksEta14=" << nESDtracksEta14
ca91a3e1	128	<< "\n";
	129	}
	130
	131	fConversionGammas=fV0Reader->GetReconstructedGammas();
	132	if(MCEvent()){
	133	fMCEvent = MCEvent();
	134	}
	135	ProcessPhotons();
	136	if(MCEvent()){
	137	ProcessMCPhotons();
	138	}
	139	PostData(1, fOutputList);
	140	}
	141
e5b6e8a6	142	///________________________________________________________________________
	143	void AliAnalysisTaskMaterial::FillMCTree(Int_t stackPos){
	144	AliStack *MCStack = fMCEvent->Stack();
	145	TParticle* candidate = (TParticle *)MCStack->Particle(stackPos);
	146	if(fConversionCuts->PhotonIsSelectedMC(candidate,MCStack,kFALSE)){
	147	Float_t gammaPt = candidate->Pt();
	148	Float_t gammaTheta = candidate->Theta();
	149	if (fStreamMaterial){
	150	(*fStreamMaterial)<<"AllGamma"
	151	<< "pt=" << gammaPt
	152	// << "p=" << gammaP
	153	<< "theta=" << gammaTheta
	154	<< "\n";
	155	}
	156	}
	157	if(fConversionCuts->PhotonIsSelectedMC(candidate,MCStack,kTRUE)){
	158	Float_t gammaPt = candidate->Pt();
	159	Float_t gammaTheta = candidate->Theta();
	160	TParticle* daughter1 = (TParticle *)MCStack->Particle(candidate->GetFirstDaughter());
	161	TParticle* daughter2 = (TParticle *)MCStack->Particle(candidate->GetLastDaughter());
	162	TVectorF coord(5);
	163	coord(0) = daughter1->Vx();
	164	coord(1) = daughter1->Vy();
	165	coord(2) = daughter1->Vz();
	166	coord(3) = daughter1->R();
	167	coord(4) = candidate->Phi();
	168	TVectorF daughterProp(4);
	169	if (daughter1-> GetPdgCode() < 0){
	170	daughterProp(0) = daughter2->Pt();
	171	daughterProp(1) = daughter2->Theta();
	172	daughterProp(2) = daughter1->Pt();
	173	daughterProp(3) = daughter1->Theta();
	174	} else {
	175	daughterProp(0) = daughter1->Pt();
	176	daughterProp(1) = daughter1->Theta();
	177	daughterProp(2) = daughter2->Pt();
	178	daughterProp(3) = daughter2->Theta();
	179	}
	180	if (fStreamMaterial){
	181	(*fStreamMaterial)<<"ConvGammaMC"
	182	<< "pt=" << gammaPt
	183	// << "p=" << gammaP
	184	<< "theta=" << gammaTheta
	185	<< "coord.=" << &coord
	186	<< "daughterProp.=" << &daughterProp
	187	<< "\n";
	188	}
	189	} // Converted MC Gamma
	190	}
ca91a3e1	191
	192	///________________________________________________________________________
	193	void AliAnalysisTaskMaterial::ProcessMCPhotons(){
	194	// Loop over all primary MC particle
	195	AliStack *ffMCStack = fMCEvent->Stack();
	196	for(Int_t i = 0; i < ffMCStack->GetNprimary(); i++) {
	197	TParticle* particle = (TParticle *)ffMCStack->Particle(i);
	198	if (!particle) continue;
e5b6e8a6	199	if (particle->GetPdgCode() == 111 && particle->GetFirstDaughter() >= ffMCStack->GetNprimary()){
	200	// cout << "Undecayed pi0 found with mother: " << particle->GetMother(0) << endl;
	201	for (Int_t j = 0; j < 2 ; j++){
	202	FillMCTree(particle->GetDaughter(j));
ca91a3e1	203	}
e5b6e8a6	204	} else {
	205	FillMCTree(i);
	206	}
ca91a3e1	207	}
	208	}
	209
	210	///________________________________________________________________________
	211	void AliAnalysisTaskMaterial::ProcessPhotons(){
	212
	213	// Fill Histograms for QA and MC
	214	for(Int_t firstGammaIndex=0;firstGammaIndex<fConversionGammas->GetEntriesFast();firstGammaIndex++){
	215	AliAODConversionPhoton gamma=dynamic_cast<AliAODConversionPhoton>(fConversionGammas->At(firstGammaIndex));
0a2b2b4b	216	if (gamma ==NULL) continue;
ca91a3e1	217	if(!fConversionCuts->PhotonIsSelected(gamma,fESDEvent)) continue;
	218	// cout << "i= " <<firstGammaIndex << " of "<< fConversionGammas->GetEntriesFast() << endl;
	219	Float_t gammaPt = gamma->GetPhotonPt();
e5b6e8a6	220	// Float_t gammaP = gamma->GetPhotonP();
e5b6e8a6	221	Float_t gammaTheta = gamma->GetPhotonTheta();
ca91a3e1	222	Float_t gammaEta = gamma->GetPhotonEta();
ca91a3e1	223	Float_t gammaChi2NDF = gamma->GetChi2perNDF();
e5b6e8a6	224	// Float_t gammaX = gamma->GetConversionX();
	225	// Float_t gammaY = gamma->GetConversionY();
	226	Float_t gammaZ = gamma->GetConversionZ();
ca91a3e1	227	Float_t gammaR = gamma->GetConversionRadius();
e5b6e8a6	228	Float_t gammaPhi = gamma->GetPhotonPhi();
	229
	230	TVectorF coord(5);
	231	coord(0) = gamma->GetConversionX();
	232	coord(1) = gamma->GetConversionY();
	233	coord(2) = gamma->GetConversionZ();
	234	coord(3) = gamma->GetConversionRadius();
	235	coord(4) = gamma->GetPhotonPhi();
	236	TVectorF daughterProp(4);
	237	AliESDtrack * negTrack = fConversionCuts->GetESDTrack(fESDEvent, gamma->GetTrackLabelNegative());
	238	AliESDtrack * posTrack = fConversionCuts->GetESDTrack(fESDEvent, gamma->GetTrackLabelPositive());
	239	daughterProp(0) = posTrack->Pt();
	240	daughterProp(1) = posTrack->Theta();
	241	daughterProp(2) = negTrack->Pt();
	242	daughterProp(3) = negTrack->Theta();
	243
ca91a3e1	244	UInt_t kind = 9;
	245	if(fMCEvent){
	246	// cout << "generating MC stack"<< endl;
	247	AliStack *fMCStack = fMCEvent->Stack();
	248	if (!fMCStack) continue;
	249	TParticle *posDaughter = gamma->GetPositiveMCDaughter(fMCStack);
	250	TParticle *negDaughter = gamma->GetNegativeMCDaughter(fMCStack);
	251	// cout << "generate Daughters: "<<posDaughter << "\t" << negDaughter << endl;
	252
	253	if(posDaughter == NULL \|\| negDaughter == NULL){
	254	kind = 9; // garbage
	255	// cout << "one of the daughters not available" << endl;
	256	} else if(posDaughter->GetMother(0) != negDaughter->GetMother(0) \|\| (posDaughter->GetMother(0) == negDaughter->GetMother(0) && posDaughter->GetMother(0) ==-1)){
	257	// Not Same Mother == Combinatorial Bck
	258	kind = 1;
	259	// cout << "not the same mother" << endl;
	260	Int_t pdgCodePos;
	261	if (posDaughter->GetPdgCode()) pdgCodePos = posDaughter->GetPdgCode(); else continue;
	262	Int_t pdgCodeNeg;
	263	if (negDaughter->GetPdgCode()) pdgCodeNeg = negDaughter->GetPdgCode(); else continue;
	264	// cout << "PDG codes daughters: " << pdgCodePos << "\t" << pdgCodeNeg << endl;
	265	if(TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==11)
	266	kind = 10; //Electron Combinatorial
	267	if(TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==11 && (posDaughter->GetMother(0) == negDaughter->GetMother(0) && posDaughter->GetMother(0) ==-1))
	268	kind = 15; //direct Electron Combinatorial
	269
	270	if(TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==211)
	271	kind = 11; //Pion Combinatorial
	272	if((TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==2212) \|\|
	273	(TMath::Abs(pdgCodePos)==2212 && TMath::Abs(pdgCodeNeg)==211))
	274	kind = 12; //Pion, Proton Combinatorics
	275	if((TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==11) \|\|
	276	(TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==211))
	277	kind = 13; //Pion, Electron Combinatorics
	278	if (TMath::Abs(pdgCodePos)==321 \|\| TMath::Abs(pdgCodeNeg)==321)
	279	kind = 14; //Kaon combinatorics
	280
	281	} else {
	282	// cout << "same mother" << endl;
	283	Int_t pdgCodePos;
	284	if (posDaughter->GetPdgCode()) pdgCodePos = posDaughter->GetPdgCode(); else continue;
	285	Int_t pdgCodeNeg;
	286	if (negDaughter->GetPdgCode()) pdgCodeNeg = negDaughter->GetPdgCode(); else continue;
	287	// cout << "PDG codes daughters: " << pdgCodePos << "\t" << pdgCodeNeg << endl;
	288	Int_t pdgCode;
	289	if (gamma->GetMCParticle(fMCStack)->GetPdgCode()) pdgCode = gamma->GetMCParticle(fMCStack)->GetPdgCode(); else continue;
	290	// cout << "PDG code: " << pdgCode << endl;
	291	if(TMath::Abs(pdgCodePos)!=11 \|\| TMath::Abs(pdgCodeNeg)!=11)
	292	kind = 2; // combinatorics from hadronic decays
	293	else if ( !(pdgCodeNeg==pdgCodePos)){
	294	TParticle *truePhotonCanditate = gamma->GetMCParticle(fMCStack);
e5b6e8a6	295	Int_t motherLabelPhoton = truePhotonCanditate->GetMother(0);
ca91a3e1	296	if(pdgCode == 111)
	297	kind = 3; // pi0 Dalitz
	298	else if (pdgCode == 221)
	299	kind = 4; // eta Dalitz
	300	else if (!(negDaughter->GetUniqueID() != 5 \|\| posDaughter->GetUniqueID() !=5)){
e5b6e8a6	301	if(pdgCode == 22 && motherLabelPhoton < fMCStack->GetNprimary()){
ca91a3e1	302	kind = 0; // primary photons
	303	} else if (pdgCode == 22){
	304	kind = 5; //secondary photons
	305	}
	306	if(pdgCode == 22){
	307	Float_t mcPt = truePhotonCanditate->Pt();
	308	Float_t mcR = gamma->GetNegativeMCDaughter(fMCStack)->R();
	309	Float_t mcZ = gamma->GetNegativeMCDaughter(fMCStack)->Vz();
	310	Float_t mcPhi = gamma->GetNegativeMCDaughter(fMCStack)->Phi();
	311	Float_t mcEta = gamma->GetNegativeMCDaughter(fMCStack)->Eta();
	312	if (fStreamResolution){
	313	(*fStreamResolution)<<"Resolution"
	314	<< "ESDpt=" << gammaPt
	315	<< "ESDphi=" << gammaPhi
	316	<< "ESDeta=" << gammaEta
	317	<< "ESDr="<< gammaR
	318	<< "ESDz="<< gammaZ
	319	<< "MCpt=" << mcPt
	320	<< "MCphi=" << mcPhi
	321	<< "MCeta=" << mcEta
	322	<< "MCr="<< mcR
	323	<< "MCz="<< mcZ
	324	<< "chi2ndf=" << gammaChi2NDF
	325	<< "\n";
	326	}
	327	}
	328	} else kind = 9; //garbage
	329	} else kind = 9; //garbage
	330	}
	331	// cout << gammaPt << "\t" << gammaP<< "\t" << gammaEta<< "\t" <<gammaChi2NDF << "\t" << gammaX<< "\t" <<gammaY << "\t" << gammaZ<< "\t" << gammaR<< "\t" << gammaPhi<< "\t" <<kind << endl;
	332
	333	if (fStreamMaterial){
	334	(*fStreamMaterial)<<"ConvPointRec"
	335	<< "pt=" << gammaPt
e5b6e8a6	336	<< "theta=" << gammaTheta
ca91a3e1	337	<< "chi2ndf=" << gammaChi2NDF
ca91a3e1	338	<< "kind=" << kind
e5b6e8a6	339	<< "coord.=" << &coord
e5b6e8a6	340	<< "daughterProp.=" << &daughterProp
ca91a3e1	341	<< "\n";
	342	}
	343	} else {
	344	if (fStreamMaterial){
	345	(*fStreamMaterial)<<"ConvPointRec"
	346	<< "pt=" << gammaPt
e5b6e8a6	347	<< "theta=" << gammaTheta
ca91a3e1	348	<< "chi2ndf=" << gammaChi2NDF
e5b6e8a6	349	<< "coord.=" << &coord
e5b6e8a6	350	<< "daughterProp.=" << &daughterProp
ca91a3e1	351	<< "\n";
	352	}
	353	}
	354	}
	355	}
	356
	357	//________________________________________________________________________
ae947965	358	Int_t AliAnalysisTaskMaterial::CountTracks09(){
e5b6e8a6	359	Int_t fNumberOfESDTracks = 0;
ae947965	360	if(fInputEvent->IsA()==AliESDEvent::Class()){
	361	// Using standard function for setting Cuts
	362	Bool_t selectPrimaries=kTRUE;
	363	AliESDtrackCuts *EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
	364	EsdTrackCuts->SetMaxDCAToVertexZ(2);
	365	EsdTrackCuts->SetEtaRange(-0.9, 0.9);
	366	EsdTrackCuts->SetPtRange(0.15);
	367
	368	for(Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++){
	369	AliESDtrack* curTrack = (AliESDtrack*) fInputEvent->GetTrack(iTracks);
	370	if(!curTrack) continue;
	371	// if(fMCEvent && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){
	372	// if(!((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(abs(curTrack->GetLabel()), fMCStack)) continue;
	373	// }
	374	if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
	375	}
	376	delete EsdTrackCuts;
	377	EsdTrackCuts=0x0;
	378	}
	379	else if(fInputEvent->IsA()==AliAODEvent::Class()){
	380	for(Int_t iTracks = 0; iTracks<fInputEvent->GetNumberOfTracks(); iTracks++){
	381	AliAODTrack* curTrack = (AliAODTrack*) fInputEvent->GetTrack(iTracks);
	382	if(!curTrack->IsPrimaryCandidate()) continue;
	383	if(abs(curTrack->Eta())>0.9) continue;
	384	if(curTrack->Pt()<0.15) continue;
	385	if(abs(curTrack->ZAtDCA())>2) continue;
	386	fNumberOfESDTracks++;
	387	}
e5b6e8a6	388	}
e5b6e8a6	389
	390	return fNumberOfESDTracks;
	391	}
	392
ae947965	393	Int_t AliAnalysisTaskMaterial::CountTracks0914(){
e5b6e8a6	394
e5b6e8a6	395	// Using standard function for setting Cuts ; We use TPCOnlyTracks for outer eta region
a280ac15	396	//Bool_t selectPrimaries=kTRUE;
e5b6e8a6	397	// EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
ae947965	398	Int_t fNumberOfESDTracks = 0;
	399	if(fInputEvent->IsA()==AliESDEvent::Class()){
	400	// Using standard function for setting Cuts
	401	AliESDtrackCuts *EsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
	402	EsdTrackCuts->SetMaxDCAToVertexZ(5);
	403	EsdTrackCuts->SetEtaRange(0.9, 1.4);
	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;
	409	// if(fMCEvent && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){
	410	// if(!((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(abs(curTrack->GetLabel()), fMCStack)) continue;
	411	// }
	412	if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
	413	}
	414	EsdTrackCuts->SetEtaRange(-1.4, -0.9);
	415	for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
	416	AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
	417	if(!curTrack) continue;
	418	if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
	419	}
	420	delete EsdTrackCuts;
	421	EsdTrackCuts=0x0;
e5b6e8a6	422	}
ae947965	423	else if(fInputEvent->IsA()==AliAODEvent::Class()){
	424	for(Int_t iTracks = 0; iTracks<fInputEvent->GetNumberOfTracks(); iTracks++){
	425	AliAODTrack* curTrack = (AliAODTrack*) fInputEvent->GetTrack(iTracks);
	426	// if(!curTrack->IsPrimaryCandidate()) continue;
	427	if(abs(curTrack->Eta())<0.9 \|\| abs(curTrack->Eta())>1.4 ) continue;
	428	if(curTrack->Pt()<0.15) continue;
	429	if(abs(curTrack->ZAtDCA())>5) continue;
	430	fNumberOfESDTracks++;
	431	}
e5b6e8a6	432	}
ae947965	433
ca91a3e1	434	return fNumberOfESDTracks;
	435	}
	436
ca91a3e1	437	//________________________________________________________________________
	438	void AliAnalysisTaskMaterial::Terminate(Option_t *)
	439	{
0a2b2b4b	440	// if (fStreamMaterial){
	441	// fStreamMaterial->GetFile()->Write();
	442	// }
	443	// if (fStreamResolution){
	444	// fStreamResolution->GetFile()->Write();
	445	// }
ca91a3e1	446	}