[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / AliAnalysisTaskSpectraBoth.cxx

/**************************************************************************\r
 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
 *                                                                        *\r
 * Author: The ALICE Off-line Project.                                    *\r
 * Contributors are mentioned in the code where appropriate.              *\r
 *                                                                        *\r
 * Permission to use, copy, modify and distribute this software and its   *\r
 * documentation strictly for non-commercial purposes is hereby granted   *\r
 * without fee, provided that the above copyright notice appears in all   *\r
 * copies and that both the copyright notice and this permission notice   *\r
 * appear in the supporting documentation. The authors make no claims     *\r
 * about the suitability of this software for any purpose. It is          *\r
 * provided "as is" without express or implied warranty.                  *\r
 **************************************************************************/\r
\r
//-----------------------------------------------------------------\r
//         AliAnalysisTaskSpectraBoth class\r
//-----------------------------------------------------------------\r
\r
#include "TChain.h"\r
#include "TTree.h"\r
#include "TLegend.h"\r
#include "TH1F.h"\r
#include "TH2F.h"\r
#include "TCanvas.h"\r
#include "AliAnalysisTask.h"\r
#include "AliAnalysisManager.h"\r
#include "AliAODTrack.h"\r
#include "AliAODMCParticle.h"\r
#include "AliVParticle.h"\r
#include "AliAODEvent.h"\r
#include "AliAODInputHandler.h"\r
#include "AliAnalysisTaskSpectraBoth.h"\r
#include "AliAnalysisTaskESDfilter.h"\r
#include "AliAnalysisDataContainer.h"\r
#include "AliSpectraBothHistoManager.h"\r
#include "AliSpectraBothTrackCuts.h"\r
#include "AliSpectraBothEventCuts.h"\r
#include "AliCentrality.h"\r
#include "TProof.h"\r
#include "AliPID.h"\r
#include "AliVEvent.h"\r
#include "AliESDEvent.h"\r
#include "AliPIDResponse.h"\r
#include "AliStack.h"\r
#include "AliSpectraBothPID.h"\r
#include "AliGenEventHeader.h"	\r
#include <TMCProcess.h>\r
\r
#include <iostream>\r
\r
\r
\r
\r
using namespace AliSpectraNameSpaceBoth;\r
using namespace std;\r
\r
ClassImp(AliAnalysisTaskSpectraBoth)\r
\r
//________________________________________________________________________\r
AliAnalysisTaskSpectraBoth::AliAnalysisTaskSpectraBoth(const char *name) : AliAnalysisTaskSE(name), fAOD(0), fHistMan(0), fTrackCuts(0), fEventCuts(0),  fPID(0), fIsMC(0), fNRebin(0),fUseMinSigma(0),fCuts(0),fdotheMCLoopAfterEventCuts(0)\r
{\r
  // Default constructor\r
  \r
  DefineInput(0, TChain::Class());\r
  DefineOutput(1, AliSpectraBothHistoManager::Class());\r
  DefineOutput(2, AliSpectraBothEventCuts::Class());\r
  DefineOutput(3, AliSpectraBothTrackCuts::Class());\r
  DefineOutput(4, AliSpectraBothPID::Class());\r
  fNRebin=0;\r
  \r
}\r
//________________________________________________________________________\r
//________________________________________________________________________\r
void AliAnalysisTaskSpectraBoth::UserCreateOutputObjects()\r
{\r
  // create output objects\r
  fHistMan = new AliSpectraBothHistoManager("SpectraHistos",fNRebin);\r
\r
  if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");\r
  if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");\r
  if (!fPID)       AliFatal("PID object should be set in the steering macro");\r
  fTrackCuts->SetAliESDtrackCuts(fCuts);\r
  PostData(1, fHistMan  );\r
  PostData(2, fEventCuts);\r
  PostData(3, fTrackCuts);\r
  PostData(4, fPID      );\r
\r
}\r
//________________________________________________________________________\r
void AliAnalysisTaskSpectraBoth::UserExec(Option_t *)\r
{\r
  // main event loop\r
	Int_t ifAODEvent=AliSpectraBothTrackCuts::kotherobject;\r
	fAOD = dynamic_cast<AliVEvent*>(InputEvent());\r
   \r
  //	AliESDEvent* esdevent=0x0;\r
 // 	AliAODEvent* aodevent=0x0;\r
   \r
   	TString nameoftrack(fAOD->ClassName());  \r
    	if(!nameoftrack.CompareTo("AliESDEvent"))\r
    	{\r
		ifAODEvent=AliSpectraBothTrackCuts::kESDobject;\r
		//esdevent=dynamic_cast<AliESDEvent*>(fAOD);\r
	}\r
	else if(!nameoftrack.CompareTo("AliAODEvent"))\r
	{\r
		ifAODEvent=AliSpectraBothTrackCuts::kAODobject;\r
		//aodevent=dynamic_cast<AliAODEvent*>(fAOD);\r
	}\r
	else\r
		AliFatal("Not processing AODs or ESDS") ;\r
	if(fdotheMCLoopAfterEventCuts)\r
  		if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection\r
 \r
  \r
  	TClonesArray *arrayMC = 0;\r
  	Int_t npar=0;\r
  	AliStack* stack=0x0;\r
	 Double_t mcZ=-100;\r
\r
 	if (fIsMC)\r
  	{\r
		TArrayF mcVertex(3);\r
  		mcVertex[0]=9999.; mcVertex[1]=9999.; mcVertex[2]=9999.;\r
		AliMCEvent* mcEvent=(AliMCEvent*)MCEvent();\r
		AliHeader* header = mcEvent->Header();\r
    		if (!header) \r
		{\r
      			AliDebug(AliLog::kError, "Header not available");\r
      			return;\r
    		}\r
	\r
		AliGenEventHeader* genHeader = header->GenEventHeader();\r
    		if(genHeader)\r
		{\r
			genHeader->PrimaryVertex(mcVertex);\r
  			mcZ=mcVertex[2];\r
		}\r
		  if(ifAODEvent==AliSpectraBothTrackCuts::kAODobject)\r
		  {\r
			  arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());\r
			  if (!arrayMC) \r
			  {\r
					AliFatal("Error: MC particles branch not found!\n");\r
			  }\r
			  Int_t nMC = arrayMC->GetEntries();\r
			  for (Int_t iMC = 0; iMC < nMC; iMC++)\r
			  {\r
				  AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);\r
				  if(!partMC->Charge()) continue;//Skip neutrals\r
				  //if(partMC->Eta() > fTrackCuts->GetEtaMin() && partMC->Eta() < fTrackCuts->GetEtaMax()){//charged hadron are filled inside the eta acceptance\r
				  //Printf("%f     %f-%f",partMC->Eta(),fTrackCuts->GetEtaMin(),fTrackCuts->GetEtaMax());\r
				  if(partMC->Eta() > fTrackCuts->GetEtaMin() && partMC->Eta() < fTrackCuts->GetEtaMax())\r
						fHistMan->GetPtHistogram(kHistPtGen)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary());					 				 \r
				  //rapidity cut\r
				  if(partMC->Y() > fTrackCuts->GetYMax()|| partMC->Y() < fTrackCuts->GetYMin() ) \r
					continue;	\r
				  if(partMC->IsPhysicalPrimary())\r
				 	 npar++;    \r
				  // check for true PID + and fill P_t histos \r
				  Int_t charge = partMC->Charge() > 0 ? kChPos : kChNeg ;\r
				  Int_t id = fPID->GetParticleSpecie(partMC);\r
				  if(id != kSpUndefined) \r
				  {\r
					fHistMan->GetHistogram2D(kHistPtGenTruePrimary,id,charge)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary());\r
				  }\r
			  }\r
		  }\r
		  if(ifAODEvent==AliSpectraBothTrackCuts::kESDobject)\r
		  {\r
		  	//AliMCEvent* mcEvent  = (AliMCEvent*) MCEvent();\r
			//Printf("MC particles: %d", mcEvent->GetNumberOfTracks());		\r
			if (!mcEvent) \r
			{\r
				AliFatal("Error: MC particles branch not found!\n");\r
			}\r
			stack = mcEvent->Stack();\r
			Int_t nMC = stack->GetNtrack();\r
			for (Int_t iMC = 0; iMC < nMC; iMC++)\r
			{\r
				 \r
				TParticle *partMC = stack->Particle(iMC);\r
				\r
				if(!partMC)	\r
					continue;	\r
	\r
				if(!partMC->GetPDG(0))\r
					continue;\r
				if(TMath::Abs(partMC->GetPDG(0)->Charge()/3.0)<0.01) \r
					continue;//Skip neutrals\r
			 	if(partMC->Eta() > fTrackCuts->GetEtaMin() && partMC->Eta() < fTrackCuts->GetEtaMax())\r
					fHistMan->GetPtHistogram(kHistPtGen)->Fill(partMC->Pt(),stack->IsPhysicalPrimary(iMC));\r
				if(partMC->Y()   > fTrackCuts->GetYMax() ||partMC->Y()   < fTrackCuts->GetYMin()  ) \r
					continue;\r
				if(stack->IsPhysicalPrimary(iMC))\r
					 npar++;    \r
				  // check for true PID + and fill P_t histos \r
				Int_t charge = partMC->GetPDG(0)->Charge()/3.0 > 0 ? kChPos : kChNeg ;\r
				Int_t id = fPID->GetParticleSpecie(partMC);\r
				if(id != kSpUndefined) \r
				{\r
					fHistMan->GetHistogram2D(kHistPtGenTruePrimary,id,charge)->Fill(partMC->Pt(),stack->IsPhysicalPrimary(iMC));\r
				}\r
			  }\r
		  }\r
	}\r
	if(!fdotheMCLoopAfterEventCuts)\r
  		if(!fEventCuts->IsSelected(fAOD,fTrackCuts,fIsMC,mcZ))return;//event selection\r
\r
  	//main loop on tracks\r
	Int_t ntracks=0;\r
  	//cout<<fAOD->GetNumberOfTracks()<<endl;\r
  	for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) \r
	{\r
		AliVTrack* track = dynamic_cast<AliVTrack*>(fAOD->GetTrack(iTracks));\r
      		AliAODTrack* aodtrack=0;\r
  		AliESDtrack* esdtrack=0;\r
  		Float_t dca=-999.;\r
  		if(ifAODEvent==AliSpectraBothTrackCuts::kESDobject)\r
  		{\r
			esdtrack=dynamic_cast<AliESDtrack*>(track);\r
			if(!esdtrack)\r
				continue;\r
			Float_t dcaz=0.0;\r
			esdtrack->GetImpactParameters(dca,dcaz);\r
  		}\r
  		else if (ifAODEvent==AliSpectraBothTrackCuts::kAODobject)\r
  		{\r
			aodtrack=dynamic_cast<AliAODTrack*>(track);\r
			if(!aodtrack)\r
				continue;		\r
			dca=aodtrack->DCA();\r
  		}\r
  		else\r
			continue;\r
    		if (!fTrackCuts->IsSelected(track,kTRUE)) \r
			continue;\r
    		ntracks++;\r
    		fPID->FillQAHistos(fHistMan, track, fTrackCuts);\r
    \r
    		//calculate DCA for AOD track\r
    		if(dca==-999.)\r
		{// track->DCA() does not work in old AOD production\r
     	 		Double_t d[2], covd[3];\r
      			AliVTrack* track_clone=(AliVTrack*)track->Clone("track_clone"); // need to clone because PropagateToDCA updates the track parameters\r
      			Bool_t isDCA = track_clone->PropagateToDCA(fAOD->GetPrimaryVertex(),fAOD->GetMagneticField(),9999.,d,covd);\r
      			delete track_clone;\r
      			if(!isDCA)\r
				d[0]=-999.;\r
      			dca=d[0];\r
    		}\r
     		fHistMan->GetPtHistogram(kHistPtRec)->Fill(track->Pt(),dca);  // PT histo\r
    		// get identity and charge\r
    		Bool_t rec[3]={false,false,false};\r
    		Int_t idRec  = fPID->GetParticleSpecie(fHistMan,track, fTrackCuts,rec);\r
		for(int irec=kSpPion;irec<kNSpecies;irec++)\r
    		{\r
   \r
			if(fUseMinSigma)\r
			{\r
				if(irec>kSpPion)\r
					break;\r
			}\r
			else\r
			{	\r
				if(!rec[irec]) \r
					idRec = kSpUndefined;\r
				else	\r
					idRec=irec;\r
			}		\r
   \r
			Int_t charge = track->Charge() > 0 ? kChPos : kChNeg;\r
		\r
			// Fill histograms, only if inside y and nsigma acceptance\r
			if(idRec != kSpUndefined && fTrackCuts->CheckYCut ((BothParticleSpecies_t)idRec))\r
				fHistMan->GetHistogram2D(kHistPtRecSigma,idRec,charge)->Fill(track->Pt(),dca);\r
			//can't put a continue because we still have to fill allcharged primaries, done later\r
		\r
			/* MC Part */\r
			if (arrayMC||stack) \r
			{\r
				Bool_t isPrimary           = kFALSE;\r
			 	Bool_t isSecondaryMaterial = kFALSE; \r
			 	Bool_t isSecondaryWeak     = kFALSE; \r
			 	Int_t idGen     =kSpUndefined;\r
			 	Int_t pdgcode=0;\r
				if (ifAODEvent==AliSpectraBothTrackCuts::kAODobject)\r
				{\r
					AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));\r
				  	if (!partMC) \r
					{ \r
						AliError("Cannot get MC particle");\r
						continue; \r
				  	}\r
				  	// Check if it is primary, secondary from material or secondary from weak decay\r
				  	isPrimary           = partMC->IsPhysicalPrimary();\r
					isSecondaryWeak     = partMC->IsSecondaryFromWeakDecay();\r
					isSecondaryMaterial      = partMC->IsSecondaryFromMaterial();\r
					//cout<<"AOD tagging "<<isPrimary<<" "<<isSecondaryWeak<<isSecondaryMaterial<<" "<<partMC->GetMCProcessCode()<<endl;\r
\r
				  	if(!isPrimary&&!isSecondaryWeak&&!isSecondaryMaterial) \r
				  	{\r
						AliError("old tagging");\r
						Int_t mfl=-999,codemoth=-999;\r
						Int_t indexMoth=partMC->GetMother(); // FIXME ignore fakes? TO BE CHECKED, on ESD is GetFirstMother()\r
						if(indexMoth>=0)\r
						{//is not fake\r
					  		AliAODMCParticle* moth = (AliAODMCParticle*) arrayMC->At(indexMoth);\r
					  		codemoth = TMath::Abs(moth->GetPdgCode());\r
					  		mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));\r
						}\r
						//Int_t uniqueID = partMC->GetUniqueID();\r
						//cout<<"uniqueID: "<<partMC->GetUniqueID()<<"       "<<kPDecay<<endl;\r
						//cout<<"status: "<<partMC->GetStatus()<<"       "<<kPDecay<<endl;\r
						// if(uniqueID == kPDecay)Printf("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");\r
						if(mfl==3) \r
							isSecondaryWeak     = kTRUE; // add if(partMC->GetStatus() & kPDecay)? FIXME\r
						else       \r
							isSecondaryMaterial = kTRUE;\r
				  	}\r
					//cout<<"AOD 2 tagging "<<isPrimary<<" "<<isSecondaryWeak<<isSecondaryMaterial<<" "<<partMC->GetMCProcessCode()<<endl;\r
\r
\r
				  	idGen     = fPID->GetParticleSpecie(partMC);\r
				  	pdgcode=partMC->GetPdgCode(); \r
				}\r
				else if (ifAODEvent==AliSpectraBothTrackCuts::kESDobject)\r
				{\r
					TParticle *partMC =stack->Particle(TMath::Abs(track->GetLabel()));\r
					if (!partMC) \r
					{ \r
						AliError("Cannot get MC particle");\r
						continue; \r
				  	}\r
				  	isPrimary           = stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()));\r
					isSecondaryWeak     = stack->IsSecondaryFromWeakDecay(TMath::Abs(track->GetLabel()));\r
					isSecondaryMaterial      = stack->IsSecondaryFromMaterial(TMath::Abs(track->GetLabel()));\r
					//cout<<"ESD tagging "<<isPrimary<<" "<<isSecondaryWeak<<isSecondaryMaterial<<endl;\r
					\r
				   	idGen     = fPID->GetParticleSpecie(partMC);\r
				   	pdgcode=partMC->GetPdgCode(); \r
				}\r
				else\r
					return;\r
			\r
			//	  cout<<isPrimary<<" "<<isSecondaryWeak<<" "<<isSecondaryMaterial<<endl;\r
			//	  cout<<" functions "<<partMC->IsPhysicalPrimary()<<" "<<partMC->IsSecondaryFromWeakDecay()<<" "<<partMC->IsSecondaryFromMaterial()<<endl;\r
		  \r
		  		if (isPrimary&&irec==kSpPion)\r
					fHistMan->GetPtHistogram(kHistPtRecPrimaryAll)->Fill(track->Pt(),dca);  // PT histo of primaries\r
		  \r
		  	//nsigma cut (reconstructed nsigma)\r
		  		if(idRec == kSpUndefined) \r
					continue;\r
		  \r
		  	// rapidity cut (reconstructed pt and identity)\r
		 		 if(!fTrackCuts->CheckYCut ((BothParticleSpecies_t)idRec)) continue;\r
		  \r
		  // Get true ID\r
		  \r
		  //if(TMath::Abs(partMC->Y())   > fTrackCuts->GetY()  ) continue;	    // FIXME: do we need a rapidity cut on the generated?\r
		  // Fill histograms for primaries\r
		  \r
		 		 if (idRec == idGen) fHistMan->GetHistogram2D(kHistPtRecTrue,  idGen, charge)->Fill(track->Pt(),dca); \r
		  \r
		  		if (isPrimary) \r
				{\r
					fHistMan->GetHistogram2D(kHistPtRecSigmaPrimary, idRec, charge)->Fill(track->Pt(),dca); \r
					if(idGen != kSpUndefined) \r
					{\r
		  				fHistMan->GetHistogram2D(kHistPtRecPrimary,      idGen, charge)->Fill(track->Pt(),dca);\r
		  				if (idRec == idGen) \r
							fHistMan->GetHistogram2D(kHistPtRecTruePrimary,  idGen, charge)->Fill(track->Pt(),dca); \r
					}\r
		  		}\r
		 		 //25th Apr - Muons are added to Pions -- FIXME\r
		 		if ( pdgcode == 13 && idRec == kSpPion) \r
				{ \r
					fHistMan->GetPtHistogram(kHistPtRecTrueMuonPlus)->Fill(track->Pt(),dca); \r
					if(isPrimary)\r
		  				fHistMan->GetPtHistogram(kHistPtRecTruePrimaryMuonPlus)->Fill(track->Pt(),dca); \r
		  		}\r
		  		if ( pdgcode == -13 && idRec == kSpPion) \r
				{ \r
					fHistMan->GetPtHistogram(kHistPtRecTrueMuonMinus)->Fill(track->Pt(),dca); \r
					if (isPrimary) \r
					{\r
		  				fHistMan->GetPtHistogram(kHistPtRecTruePrimaryMuonMinus)->Fill(track->Pt(),dca); \r
					}\r
		  		}\r
		  \r
		 		 ///..... END FIXME\r
		  \r
		  		// Fill secondaries\r
		  		if(isSecondaryWeak    )  \r
					fHistMan->GetHistogram2D(kHistPtRecSigmaSecondaryWeakDecay, idRec, charge)->Fill(track->Pt(),dca);\r
		  		if(isSecondaryMaterial)  \r
					fHistMan->GetHistogram2D(kHistPtRecSigmaSecondaryMaterial , idRec, charge)->Fill(track->Pt(),dca);\r
		  \r
			}//end if(arrayMC)\r
		}\r
	\r
	\r
  	} // end loop on tracks\r
 // cout<< ntracks<<endl;\r
  fHistMan->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->Fill(npar,ntracks);\r
  PostData(1, fHistMan  );\r
  PostData(2, fEventCuts);\r
  PostData(3, fTrackCuts);\r
  PostData(4, fPID      );\r
}\r
\r
//_________________________________________________________________\r
void   AliAnalysisTaskSpectraBoth::Terminate(Option_t *)\r
{\r
  // Terminate\r
}\r
Commit	Line	Data
239a080a	1	/**************************************************************************\r
	2	* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
	3	* *\r
	4	* Author: The ALICE Off-line Project. *\r
	5	* Contributors are mentioned in the code where appropriate. *\r
	6	* *\r
	7	* Permission to use, copy, modify and distribute this software and its *\r
	8	* documentation strictly for non-commercial purposes is hereby granted *\r
	9	* without fee, provided that the above copyright notice appears in all *\r
	10	* copies and that both the copyright notice and this permission notice *\r
	11	* appear in the supporting documentation. The authors make no claims *\r
	12	* about the suitability of this software for any purpose. It is *\r
	13	* provided "as is" without express or implied warranty. *\r
	14	**************************************************************************/\r
	15	\r
	16	//-----------------------------------------------------------------\r
	17	// AliAnalysisTaskSpectraBoth class\r
	18	//-----------------------------------------------------------------\r
	19	\r
	20	#include "TChain.h"\r
	21	#include "TTree.h"\r
	22	#include "TLegend.h"\r
	23	#include "TH1F.h"\r
	24	#include "TH2F.h"\r
	25	#include "TCanvas.h"\r
	26	#include "AliAnalysisTask.h"\r
	27	#include "AliAnalysisManager.h"\r
	28	#include "AliAODTrack.h"\r
	29	#include "AliAODMCParticle.h"\r
	30	#include "AliVParticle.h"\r
	31	#include "AliAODEvent.h"\r
	32	#include "AliAODInputHandler.h"\r
	33	#include "AliAnalysisTaskSpectraBoth.h"\r
	34	#include "AliAnalysisTaskESDfilter.h"\r
	35	#include "AliAnalysisDataContainer.h"\r
	36	#include "AliSpectraBothHistoManager.h"\r
	37	#include "AliSpectraBothTrackCuts.h"\r
	38	#include "AliSpectraBothEventCuts.h"\r
	39	#include "AliCentrality.h"\r
	40	#include "TProof.h"\r
	41	#include "AliPID.h"\r
	42	#include "AliVEvent.h"\r
	43	#include "AliESDEvent.h"\r
	44	#include "AliPIDResponse.h"\r
	45	#include "AliStack.h"\r
	46	#include "AliSpectraBothPID.h"\r
8bb435b0	47	#include "AliGenEventHeader.h" \r
239a080a	48	#include <TMCProcess.h>\r
	49	\r
	50	#include <iostream>\r
	51	\r
	52	\r
	53	\r
	54	\r
	55	using namespace AliSpectraNameSpaceBoth;\r
	56	using namespace std;\r
	57	\r
	58	ClassImp(AliAnalysisTaskSpectraBoth)\r
	59	\r
	60	//________________________________________________________________________\r
	61	AliAnalysisTaskSpectraBoth::AliAnalysisTaskSpectraBoth(const char *name) : AliAnalysisTaskSE(name), fAOD(0), fHistMan(0), fTrackCuts(0), fEventCuts(0), fPID(0), fIsMC(0), fNRebin(0),fUseMinSigma(0),fCuts(0),fdotheMCLoopAfterEventCuts(0)\r
	62	{\r
	63	// Default constructor\r
	64	\r
	65	DefineInput(0, TChain::Class());\r
	66	DefineOutput(1, AliSpectraBothHistoManager::Class());\r
	67	DefineOutput(2, AliSpectraBothEventCuts::Class());\r
	68	DefineOutput(3, AliSpectraBothTrackCuts::Class());\r
	69	DefineOutput(4, AliSpectraBothPID::Class());\r
	70	fNRebin=0;\r
	71	\r
	72	}\r
	73	//________________________________________________________________________\r
	74	//________________________________________________________________________\r
	75	void AliAnalysisTaskSpectraBoth::UserCreateOutputObjects()\r
	76	{\r
	77	// create output objects\r
	78	fHistMan = new AliSpectraBothHistoManager("SpectraHistos",fNRebin);\r
	79	\r
	80	if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");\r
	81	if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");\r
	82	if (!fPID) AliFatal("PID object should be set in the steering macro");\r
	83	fTrackCuts->SetAliESDtrackCuts(fCuts);\r
	84	PostData(1, fHistMan );\r
	85	PostData(2, fEventCuts);\r
	86	PostData(3, fTrackCuts);\r
	87	PostData(4, fPID );\r
	88	\r
	89	}\r
	90	//________________________________________________________________________\r
	91	void AliAnalysisTaskSpectraBoth::UserExec(Option_t *)\r
	92	{\r
	93	// main event loop\r
	94	Int_t ifAODEvent=AliSpectraBothTrackCuts::kotherobject;\r
	95	fAOD = dynamic_cast<AliVEvent*>(InputEvent());\r
	96	\r
	97	// AliESDEvent* esdevent=0x0;\r
	98	// AliAODEvent* aodevent=0x0;\r
	99	\r
	100	TString nameoftrack(fAOD->ClassName()); \r
	101	if(!nameoftrack.CompareTo("AliESDEvent"))\r
	102	{\r
	103	ifAODEvent=AliSpectraBothTrackCuts::kESDobject;\r
	104	//esdevent=dynamic_cast<AliESDEvent*>(fAOD);\r
	105	}\r
	106	else if(!nameoftrack.CompareTo("AliAODEvent"))\r
	107	{\r
	108	ifAODEvent=AliSpectraBothTrackCuts::kAODobject;\r
	109	//aodevent=dynamic_cast<AliAODEvent*>(fAOD);\r
	110	}\r
	111	else\r
112	AliFatal("Not processing AODs or ESDS") ;\r
113	if(fdotheMCLoopAfterEventCuts)\r
114	if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection\r
115	\r
116	\r
117	TClonesArray *arrayMC = 0;\r
118	Int_t npar=0;\r
119	AliStack* stack=0x0;\r
8bb435b0	120	Double_t mcZ=-100;\r
8bb435b0	121	\r
239a080a	122	if (fIsMC)\r
239a080a	123	{\r
8bb435b0	124	TArrayF mcVertex(3);\r
	125	mcVertex[0]=9999.; mcVertex[1]=9999.; mcVertex[2]=9999.;\r
	126	AliMCEvent* mcEvent=(AliMCEvent*)MCEvent();\r
	127	AliHeader* header = mcEvent->Header();\r
	128	if (!header) \r
	129	{\r
	130	AliDebug(AliLog::kError, "Header not available");\r
	131	return;\r
	132	}\r
	133	\r
	134	AliGenEventHeader* genHeader = header->GenEventHeader();\r
	135	if(genHeader)\r
	136	{\r
	137	genHeader->PrimaryVertex(mcVertex);\r
	138	mcZ=mcVertex[2];\r
	139	}\r
239a080a	140	if(ifAODEvent==AliSpectraBothTrackCuts::kAODobject)\r
8bb435b0	141	{\r
239a080a	142	arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());\r
	143	if (!arrayMC) \r
	144	{\r
	145	AliFatal("Error: MC particles branch not found!\n");\r
	146	}\r
	147	Int_t nMC = arrayMC->GetEntries();\r
	148	for (Int_t iMC = 0; iMC < nMC; iMC++)\r
	149	{\r
	150	AliAODMCParticle partMC = (AliAODMCParticle) arrayMC->At(iMC);\r
	151	if(!partMC->Charge()) continue;//Skip neutrals\r
	152	//if(partMC->Eta() > fTrackCuts->GetEtaMin() && partMC->Eta() < fTrackCuts->GetEtaMax()){//charged hadron are filled inside the eta acceptance\r
	153	//Printf("%f %f-%f",partMC->Eta(),fTrackCuts->GetEtaMin(),fTrackCuts->GetEtaMax());\r
	154	if(partMC->Eta() > fTrackCuts->GetEtaMin() && partMC->Eta() < fTrackCuts->GetEtaMax())\r
8bb435b0	155	fHistMan->GetPtHistogram(kHistPtGen)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary()); \r
239a080a	156	//rapidity cut\r
b3ea73e1	157	if(partMC->Y() > fTrackCuts->GetYMax()\|\| partMC->Y() < fTrackCuts->GetYMin() ) \r
b3ea73e1	158	continue; \r
239a080a	159	if(partMC->IsPhysicalPrimary())\r
	160	npar++; \r
	161	// check for true PID + and fill P_t histos \r
	162	Int_t charge = partMC->Charge() > 0 ? kChPos : kChNeg ;\r
	163	Int_t id = fPID->GetParticleSpecie(partMC);\r
	164	if(id != kSpUndefined) \r
	165	{\r
	166	fHistMan->GetHistogram2D(kHistPtGenTruePrimary,id,charge)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary());\r
	167	}\r
	168	}\r
	169	}\r
	170	if(ifAODEvent==AliSpectraBothTrackCuts::kESDobject)\r
8bb435b0	171	{\r
8bb435b0	172	//AliMCEvent* mcEvent = (AliMCEvent*) MCEvent();\r
239a080a	173	//Printf("MC particles: %d", mcEvent->GetNumberOfTracks()); \r
	174	if (!mcEvent) \r
	175	{\r
	176	AliFatal("Error: MC particles branch not found!\n");\r
	177	}\r
	178	stack = mcEvent->Stack();\r
	179	Int_t nMC = stack->GetNtrack();\r
	180	for (Int_t iMC = 0; iMC < nMC; iMC++)\r
	181	{\r
	182	\r
	183	TParticle *partMC = stack->Particle(iMC);\r
	184	\r
	185	if(!partMC) \r
	186	continue; \r
	187	\r
	188	if(!partMC->GetPDG(0))\r
	189	continue;\r
	190	if(TMath::Abs(partMC->GetPDG(0)->Charge()/3.0)<0.01) \r
	191	continue;//Skip neutrals\r
	192	if(partMC->Eta() > fTrackCuts->GetEtaMin() && partMC->Eta() < fTrackCuts->GetEtaMax())\r
	193	fHistMan->GetPtHistogram(kHistPtGen)->Fill(partMC->Pt(),stack->IsPhysicalPrimary(iMC));\r
b3ea73e1	194	if(partMC->Y() > fTrackCuts->GetYMax() \|\|partMC->Y() < fTrackCuts->GetYMin() ) \r
239a080a	195	continue;\r
	196	if(stack->IsPhysicalPrimary(iMC))\r
	197	npar++; \r
	198	// check for true PID + and fill P_t histos \r
	199	Int_t charge = partMC->GetPDG(0)->Charge()/3.0 > 0 ? kChPos : kChNeg ;\r
	200	Int_t id = fPID->GetParticleSpecie(partMC);\r
	201	if(id != kSpUndefined) \r
	202	{\r
	203	fHistMan->GetHistogram2D(kHistPtGenTruePrimary,id,charge)->Fill(partMC->Pt(),stack->IsPhysicalPrimary(iMC));\r
	204	}\r
	205	}\r
	206	}\r
	207	}\r
239a080a	208	if(!fdotheMCLoopAfterEventCuts)\r
8bb435b0	209	if(!fEventCuts->IsSelected(fAOD,fTrackCuts,fIsMC,mcZ))return;//event selection\r
239a080a	210	\r
	211	//main loop on tracks\r
	212	Int_t ntracks=0;\r
	213	//cout<<fAOD->GetNumberOfTracks()<<endl;\r
	214	for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) \r
	215	{\r
	216	AliVTrack* track = dynamic_cast<AliVTrack*>(fAOD->GetTrack(iTracks));\r
	217	AliAODTrack* aodtrack=0;\r
	218	AliESDtrack* esdtrack=0;\r
	219	Float_t dca=-999.;\r
	220	if(ifAODEvent==AliSpectraBothTrackCuts::kESDobject)\r
	221	{\r
	222	esdtrack=dynamic_cast<AliESDtrack*>(track);\r
737ad8ac	223	if(!esdtrack)\r
737ad8ac	224	continue;\r
239a080a	225	Float_t dcaz=0.0;\r
	226	esdtrack->GetImpactParameters(dca,dcaz);\r
	227	}\r
	228	else if (ifAODEvent==AliSpectraBothTrackCuts::kAODobject)\r
	229	{\r
737ad8ac	230	aodtrack=dynamic_cast<AliAODTrack*>(track);\r
	231	if(!aodtrack)\r
	232	continue; \r
239a080a	233	dca=aodtrack->DCA();\r
	234	}\r
	235	else\r
	236	continue;\r
	237	if (!fTrackCuts->IsSelected(track,kTRUE)) \r
	238	continue;\r
	239	ntracks++;\r
	240	fPID->FillQAHistos(fHistMan, track, fTrackCuts);\r
	241	\r
	242	//calculate DCA for AOD track\r
	243	if(dca==-999.)\r
	244	{// track->DCA() does not work in old AOD production\r
	245	Double_t d[2], covd[3];\r
	246	AliVTrack* track_clone=(AliVTrack*)track->Clone("track_clone"); // need to clone because PropagateToDCA updates the track parameters\r
	247	Bool_t isDCA = track_clone->PropagateToDCA(fAOD->GetPrimaryVertex(),fAOD->GetMagneticField(),9999.,d,covd);\r
	248	delete track_clone;\r
	249	if(!isDCA)\r
	250	d[0]=-999.;\r
	251	dca=d[0];\r
	252	}\r
	253	fHistMan->GetPtHistogram(kHistPtRec)->Fill(track->Pt(),dca); // PT histo\r
	254	// get identity and charge\r
	255	Bool_t rec[3]={false,false,false};\r
	256	Int_t idRec = fPID->GetParticleSpecie(fHistMan,track, fTrackCuts,rec);\r
	257	for(int irec=kSpPion;irec<kNSpecies;irec++)\r
	258	{\r
	259	\r
	260	if(fUseMinSigma)\r
	261	{\r
	262	if(irec>kSpPion)\r
	263	break;\r
	264	}\r
	265	else\r
	266	{ \r
	267	if(!rec[irec]) \r
	268	idRec = kSpUndefined;\r
	269	else \r
	270	idRec=irec;\r
	271	} \r
	272	\r
	273	Int_t charge = track->Charge() > 0 ? kChPos : kChNeg;\r
	274	\r
	275	// Fill histograms, only if inside y and nsigma acceptance\r
	276	if(idRec != kSpUndefined && fTrackCuts->CheckYCut ((BothParticleSpecies_t)idRec))\r
	277	fHistMan->GetHistogram2D(kHistPtRecSigma,idRec,charge)->Fill(track->Pt(),dca);\r
	278	//can't put a continue because we still have to fill allcharged primaries, done later\r
	279	\r
	280	/* MC Part */\r
	281	if (arrayMC\|\|stack) \r
	282	{\r
	283	Bool_t isPrimary = kFALSE;\r
	284	Bool_t isSecondaryMaterial = kFALSE; \r
	285	Bool_t isSecondaryWeak = kFALSE; \r
	286	Int_t idGen =kSpUndefined;\r
	287	Int_t pdgcode=0;\r
	288	if (ifAODEvent==AliSpectraBothTrackCuts::kAODobject)\r
	289	{\r
	290	AliAODMCParticle partMC = (AliAODMCParticle) arrayMC->At(TMath::Abs(track->GetLabel()));\r
	291	if (!partMC) \r
	292	{ \r
	293	AliError("Cannot get MC particle");\r
	294	continue; \r
	295	}\r
	296	// Check if it is primary, secondary from material or secondary from weak decay\r
297	isPrimary = partMC->IsPhysicalPrimary();\r
298	isSecondaryWeak = partMC->IsSecondaryFromWeakDecay();\r
299	isSecondaryMaterial = partMC->IsSecondaryFromMaterial();\r
300	//cout<<"AOD tagging "<<isPrimary<<" "<<isSecondaryWeak<<isSecondaryMaterial<<" "<<partMC->GetMCProcessCode()<<endl;\r
301	\r
302	if(!isPrimary&&!isSecondaryWeak&&!isSecondaryMaterial) \r
303	{\r
304	AliError("old tagging");\r
305	Int_t mfl=-999,codemoth=-999;\r
306	Int_t indexMoth=partMC->GetMother(); // FIXME ignore fakes? TO BE CHECKED, on ESD is GetFirstMother()\r
307	if(indexMoth>=0)\r
308	{//is not fake\r
309	AliAODMCParticle* moth = (AliAODMCParticle*) arrayMC->At(indexMoth);\r
310	codemoth = TMath::Abs(moth->GetPdgCode());\r
311	mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));\r
312	}\r
313	//Int_t uniqueID = partMC->GetUniqueID();\r
314	//cout<<"uniqueID: "<<partMC->GetUniqueID()<<" "<<kPDecay<<endl;\r
315	//cout<<"status: "<<partMC->GetStatus()<<" "<<kPDecay<<endl;\r
316	// if(uniqueID == kPDecay)Printf("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");\r
317	if(mfl==3) \r
318	isSecondaryWeak = kTRUE; // add if(partMC->GetStatus() & kPDecay)? FIXME\r
319	else \r
320	isSecondaryMaterial = kTRUE;\r
321	}\r
322	//cout<<"AOD 2 tagging "<<isPrimary<<" "<<isSecondaryWeak<<isSecondaryMaterial<<" "<<partMC->GetMCProcessCode()<<endl;\r
323	\r
324	\r
325	idGen = fPID->GetParticleSpecie(partMC);\r
326	pdgcode=partMC->GetPdgCode(); \r
327	}\r
328	else if (ifAODEvent==AliSpectraBothTrackCuts::kESDobject)\r
329	{\r
330	TParticle *partMC =stack->Particle(TMath::Abs(track->GetLabel()));\r
331	if (!partMC) \r
332	{ \r
333	AliError("Cannot get MC particle");\r
334	continue; \r
335	}\r
336	isPrimary = stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()));\r
337	isSecondaryWeak = stack->IsSecondaryFromWeakDecay(TMath::Abs(track->GetLabel()));\r
338	isSecondaryMaterial = stack->IsSecondaryFromMaterial(TMath::Abs(track->GetLabel()));\r
339	//cout<<"ESD tagging "<<isPrimary<<" "<<isSecondaryWeak<<isSecondaryMaterial<<endl;\r
340	\r
341	idGen = fPID->GetParticleSpecie(partMC);\r
342	pdgcode=partMC->GetPdgCode(); \r
343	}\r
344	else\r
345	return;\r
346	\r
347	// cout<<isPrimary<<" "<<isSecondaryWeak<<" "<<isSecondaryMaterial<<endl;\r
348	// cout<<" functions "<<partMC->IsPhysicalPrimary()<<" "<<partMC->IsSecondaryFromWeakDecay()<<" "<<partMC->IsSecondaryFromMaterial()<<endl;\r
349	\r
350	if (isPrimary&&irec==kSpPion)\r
94a8ed71	351	fHistMan->GetPtHistogram(kHistPtRecPrimaryAll)->Fill(track->Pt(),dca); // PT histo of primaries\r
239a080a	352	\r
	353	//nsigma cut (reconstructed nsigma)\r
	354	if(idRec == kSpUndefined) \r
	355	continue;\r
	356	\r
	357	// rapidity cut (reconstructed pt and identity)\r
	358	if(!fTrackCuts->CheckYCut ((BothParticleSpecies_t)idRec)) continue;\r
	359	\r
	360	// Get true ID\r
	361	\r
	362	//if(TMath::Abs(partMC->Y()) > fTrackCuts->GetY() ) continue; // FIXME: do we need a rapidity cut on the generated?\r
	363	// Fill histograms for primaries\r
	364	\r
	365	if (idRec == idGen) fHistMan->GetHistogram2D(kHistPtRecTrue, idGen, charge)->Fill(track->Pt(),dca); \r
	366	\r
	367	if (isPrimary) \r
	368	{\r
	369	fHistMan->GetHistogram2D(kHistPtRecSigmaPrimary, idRec, charge)->Fill(track->Pt(),dca); \r
	370	if(idGen != kSpUndefined) \r
	371	{\r
	372	fHistMan->GetHistogram2D(kHistPtRecPrimary, idGen, charge)->Fill(track->Pt(),dca);\r
	373	if (idRec == idGen) \r
	374	fHistMan->GetHistogram2D(kHistPtRecTruePrimary, idGen, charge)->Fill(track->Pt(),dca); \r
	375	}\r
	376	}\r
	377	//25th Apr - Muons are added to Pions -- FIXME\r
	378	if ( pdgcode == 13 && idRec == kSpPion) \r
	379	{ \r
	380	fHistMan->GetPtHistogram(kHistPtRecTrueMuonPlus)->Fill(track->Pt(),dca); \r
	381	if(isPrimary)\r
	382	fHistMan->GetPtHistogram(kHistPtRecTruePrimaryMuonPlus)->Fill(track->Pt(),dca); \r
	383	}\r
	384	if ( pdgcode == -13 && idRec == kSpPion) \r
	385	{ \r
	386	fHistMan->GetPtHistogram(kHistPtRecTrueMuonMinus)->Fill(track->Pt(),dca); \r
	387	if (isPrimary) \r
	388	{\r
	389	fHistMan->GetPtHistogram(kHistPtRecTruePrimaryMuonMinus)->Fill(track->Pt(),dca); \r
	390	}\r
	391	}\r
	392	\r
	393	///..... END FIXME\r
	394	\r
	395	// Fill secondaries\r
	396	if(isSecondaryWeak ) \r
	397	fHistMan->GetHistogram2D(kHistPtRecSigmaSecondaryWeakDecay, idRec, charge)->Fill(track->Pt(),dca);\r
	398	if(isSecondaryMaterial) \r
	399	fHistMan->GetHistogram2D(kHistPtRecSigmaSecondaryMaterial , idRec, charge)->Fill(track->Pt(),dca);\r
	400	\r
	401	}//end if(arrayMC)\r
	402	}\r
	403	\r
	404	\r
	405	} // end loop on tracks\r
	406	// cout<< ntracks<<endl;\r
	407	fHistMan->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->Fill(npar,ntracks);\r
	408	PostData(1, fHistMan );\r
	409	PostData(2, fEventCuts);\r
	410	PostData(3, fTrackCuts);\r
	411	PostData(4, fPID );\r
	412	}\r
	413	\r
	414	//_________________________________________________________________\r
	415	void AliAnalysisTaskSpectraBoth::Terminate(Option_t *)\r
416	{\r
417	// Terminate\r
418	}\r