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

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

//-----------------------------------------------------------------
//         AliAnalysisTaskSpectraBoth class
//-----------------------------------------------------------------

#include "TChain.h"
#include "TTree.h"
#include "TLegend.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h"
#include "TCanvas.h"
#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"
#include "AliAODTrack.h"
#include "AliAODMCParticle.h"
#include "AliVParticle.h"
#include "AliAODEvent.h"
#include "AliAODInputHandler.h"
#include "AliAnalysisTaskSpectraBoth.h"
#include "AliAnalysisTaskESDfilter.h"
#include "AliAnalysisDataContainer.h"
#include "AliSpectraBothHistoManager.h"
#include "AliSpectraBothTrackCuts.h"
#include "AliSpectraBothEventCuts.h"
#include "AliCentrality.h"
#include "TProof.h"
#include "AliPID.h"
#include "AliVEvent.h"
#include "AliESDEvent.h"
#include "AliPIDResponse.h"
#include "AliStack.h"
#include "AliSpectraBothPID.h"
#include "AliGenEventHeader.h"	
#include <TMCProcess.h>

#include <iostream>


using namespace AliSpectraNameSpaceBoth;
using namespace std;

ClassImp(AliAnalysisTaskSpectraBoth)

//________________________________________________________________________
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),
fmakePIDQAhisto(1),fMotherWDPDGcode(-1)

{
  // Default constructor
  
  DefineInput(0, TChain::Class());
  DefineOutput(1, AliSpectraBothHistoManager::Class());
  DefineOutput(2, AliSpectraBothEventCuts::Class());
  DefineOutput(3, AliSpectraBothTrackCuts::Class());
  DefineOutput(4, AliSpectraBothPID::Class());
  fNRebin=0;
  
}
//________________________________________________________________________
//________________________________________________________________________
void AliAnalysisTaskSpectraBoth::UserCreateOutputObjects()
{
  // create output objects
  fHistMan = new AliSpectraBothHistoManager("SpectraHistos",fNRebin,fmakePIDQAhisto);

  if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");
  if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");
  if (!fPID)       AliFatal("PID object should be set in the steering macro");
  fTrackCuts->SetAliESDtrackCuts(fCuts);
  fEventCuts->InitHisto();
  fTrackCuts->InitHisto();

  PostData(1, fHistMan  );
  PostData(2, fEventCuts);
  PostData(3, fTrackCuts);
  PostData(4, fPID      );
}
//________________________________________________________________________
void AliAnalysisTaskSpectraBoth::UserExec(Option_t *)
{
  // main event loop
	Int_t ifAODEvent=AliSpectraBothTrackCuts::kotherobject;
	fAOD = dynamic_cast<AliVEvent*>(InputEvent());
  //	AliESDEvent* esdevent=0x0;
 // 	AliAODEvent* aodevent=0x0;
   
   	TString nameoftrack(fAOD->ClassName());  
    	if(!nameoftrack.CompareTo("AliESDEvent"))
    	{
		ifAODEvent=AliSpectraBothTrackCuts::kESDobject;
		//esdevent=dynamic_cast<AliESDEvent*>(fAOD);
	}
	else if(!nameoftrack.CompareTo("AliAODEvent"))
	{
		ifAODEvent=AliSpectraBothTrackCuts::kAODobject;
		//aodevent=dynamic_cast<AliAODEvent*>(fAOD);
	}
	else
		AliFatal("Not processing AODs or ESDS") ;
	if(fdotheMCLoopAfterEventCuts)
  		if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection
  	TClonesArray *arrayMC = 0;
  	Int_t npar=0;
  	AliStack* stack=0x0;
	 Double_t mcZ=-100;

 	if (fIsMC)
  	{
		TArrayF mcVertex(3);
  		mcVertex[0]=9999.; mcVertex[1]=9999.; mcVertex[2]=9999.;
		AliMCEvent* mcEvent=(AliMCEvent*)MCEvent();
		if (!mcEvent) 
		{
			AliFatal("Error: MC particles branch not found!\n");
		}
		AliHeader* header = mcEvent->Header();
    		if (!header) 
		{
      			AliDebug(AliLog::kError, "Header not available");
      			return;
    		}
	
		AliGenEventHeader* genHeader = header->GenEventHeader();
    		if(genHeader)
		{
			genHeader->PrimaryVertex(mcVertex);
  			mcZ=mcVertex[2];
		}
		  if(ifAODEvent==AliSpectraBothTrackCuts::kAODobject)
		  {
			  arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
			  if (!arrayMC) 
			  {
					AliFatal("Error: MC particles branch not found!\n");
			  }
			  Int_t nMC = arrayMC->GetEntries();
			  for (Int_t iMC = 0; iMC < nMC; iMC++)
			  {
				  AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
				  if(!partMC->Charge()) continue;//Skip neutrals
				  //if(partMC->Eta() > fTrackCuts->GetEtaMin() && partMC->Eta() < fTrackCuts->GetEtaMax()){//charged hadron are filled inside the eta acceptance
				  //Printf("%f     %f-%f",partMC->Eta(),fTrackCuts->GetEtaMin(),fTrackCuts->GetEtaMax());
				  if(partMC->Eta() > fTrackCuts->GetEtaMin() && partMC->Eta() < fTrackCuts->GetEtaMax())
						fHistMan->GetPtHistogram(kHistPtGen)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary());					 				 
				  //rapidity cut
				  if(partMC->Y() > fTrackCuts->GetYMax()|| partMC->Y() < fTrackCuts->GetYMin() ) 
					continue;	
				  if(partMC->IsPhysicalPrimary())
				 	 npar++;    
				  // check for true PID + and fill P_t histos 
				  Int_t charge = partMC->Charge() > 0 ? kChPos : kChNeg ;
				  Int_t id = fPID->GetParticleSpecie(partMC);
				  if(id != kSpUndefined) 
				  {
					fHistMan->GetHistogram2D(kHistPtGenTruePrimary,id,charge)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary());
				  }
			  }
		  }
		  if(ifAODEvent==AliSpectraBothTrackCuts::kESDobject)
		  {
			stack = mcEvent->Stack();
			Int_t nMC = stack->GetNtrack();
			for (Int_t iMC = 0; iMC < nMC; iMC++)
			{
				 
				TParticle *partMC = stack->Particle(iMC);
				
				if(!partMC)	
					continue;	
	
				if(!partMC->GetPDG(0))
					continue;
				if(TMath::Abs(partMC->GetPDG(0)->Charge()/3.0)<0.01) 
					continue;//Skip neutrals
			 	if(partMC->Eta() > fTrackCuts->GetEtaMin() && partMC->Eta() < fTrackCuts->GetEtaMax())
					fHistMan->GetPtHistogram(kHistPtGen)->Fill(partMC->Pt(),stack->IsPhysicalPrimary(iMC));
				if(partMC->Y()   > fTrackCuts->GetYMax() ||partMC->Y()   < fTrackCuts->GetYMin()  ) 
					continue;
				if(stack->IsPhysicalPrimary(iMC))
					 npar++;    
				  // check for true PID + and fill P_t histos 
				Int_t charge = partMC->GetPDG(0)->Charge()/3.0 > 0 ? kChPos : kChNeg ;
				Int_t id = fPID->GetParticleSpecie(partMC);
				if(id != kSpUndefined) 
				{
					fHistMan->GetHistogram2D(kHistPtGenTruePrimary,id,charge)->Fill(partMC->Pt(),stack->IsPhysicalPrimary(iMC));
				}
			  }
		  }
	}
	if(!fdotheMCLoopAfterEventCuts)
  		if(!fEventCuts->IsSelected(fAOD,fTrackCuts,fIsMC,mcZ))return;//event selection
  	//main loop on tracks
	Int_t ntracks=0;
  	//cout<<fAOD->GetNumberOfTracks()<<endl;
  	for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) 
	{
		AliVTrack* track = dynamic_cast<AliVTrack*>(fAOD->GetTrack(iTracks));
      		AliAODTrack* aodtrack=0;
  		AliESDtrack* esdtrack=0;
  		Float_t dca=-999.;
		Float_t dcaz=-999.;
		Short_t ncls=-1; 
		Float_t chi2perndf=-1.0;
  		if(ifAODEvent==AliSpectraBothTrackCuts::kESDobject)
  		{
			esdtrack=dynamic_cast<AliESDtrack*>(track);
			if(!esdtrack)
				continue;
			esdtrack->GetImpactParameters(dca,dcaz);
			ncls=esdtrack->GetTPCNcls();
			if ( ncls > 5) 
       				chi2perndf=(esdtrack->GetTPCchi2()/Float_t(ncls - 5));
  
			else 
       				chi2perndf=-1;
    			
  		}
  		else if (ifAODEvent==AliSpectraBothTrackCuts::kAODobject)
  		{
			aodtrack=dynamic_cast<AliAODTrack*>(track);
			if(!aodtrack)
				continue;		
			dca=aodtrack->DCA();
			dcaz=aodtrack->ZAtDCA();
			ncls=aodtrack->GetTPCNcls();
			chi2perndf=aodtrack->Chi2perNDF();
  		}
  		else
			continue;
    		if (!fTrackCuts->IsSelected(track,kTRUE)) 
			continue;
    		ntracks++;
		if(fmakePIDQAhisto)
    			fPID->FillQAHistos(fHistMan, track, fTrackCuts);
    		
		//calculate DCA for AOD track
    		if(dca==-999.)
		{// track->DCA() does not work in old AOD production
     	 		Double_t d[2], covd[3];
      			AliVTrack* track_clone=(AliVTrack*)track->Clone("track_clone"); // need to clone because PropagateToDCA updates the track parameters
      			Bool_t isDCA = track_clone->PropagateToDCA(fAOD->GetPrimaryVertex(),fAOD->GetMagneticField(),9999.,d,covd);
      			delete track_clone;
      			if(!isDCA)
				d[0]=-999.;
      			dca=d[0];
			dcaz=d[1];	
    		}
	
     		fHistMan->GetPtHistogram(kHistPtRec)->Fill(track->Pt(),dca);  // PT histo
    		// get identity and charge
    		Bool_t rec[3]={false,false,false};
    		Int_t idRec  = fPID->GetParticleSpecie(fHistMan,track, fTrackCuts,rec);
		for(int irec=kSpPion;irec<kNSpecies;irec++)
    		{
   
			if(fUseMinSigma)
			{
				if(irec>kSpPion)
					break;
			}
			else
			{	
				if(!rec[irec]) 
					idRec = kSpUndefined;
				else	
					idRec=irec;
			}		
   
			Int_t charge = track->Charge() > 0 ? kChPos : kChNeg;
		
			// Fill histograms, only if inside y and nsigma acceptance
			if(idRec != kSpUndefined && fTrackCuts->CheckYCut ((BothParticleSpecies_t)idRec))
			{
				fHistMan->GetHistogram2D(kHistPtRecSigma,idRec,charge)->Fill(track->Pt(),dca);
				fTrackCuts->GetHistoDCAzQA()->Fill(idRec,track->Pt(),dcaz);
				fTrackCuts->GetHistoNclustersQA()->Fill(idRec,track->Pt(),ncls);
				fTrackCuts->GetHistochi2perNDFQA()->Fill(idRec,track->Pt(),chi2perndf);
			}
			//can't put a continue because we still have to fill allcharged primaries, done later
		
			/* MC Part */
			if (arrayMC||stack) 
			{
				Bool_t isPrimary           = kFALSE;
			 	Bool_t isSecondaryMaterial = kFALSE; 
			 	Bool_t isSecondaryWeak     = kFALSE; 
			 	Int_t idGen     =kSpUndefined;
			 	Int_t pdgcode=0;
				Int_t motherpdg=-1;
				if (ifAODEvent==AliSpectraBothTrackCuts::kAODobject)
				{
					AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));
				  	if (!partMC) 
					{ 
						AliError("Cannot get MC particle");
						continue; 
				  	}
				  	// Check if it is primary, secondary from material or secondary from weak decay
				  	isPrimary           = partMC->IsPhysicalPrimary();
					isSecondaryWeak     = partMC->IsSecondaryFromWeakDecay();
					isSecondaryMaterial      = partMC->IsSecondaryFromMaterial();
					//cout<<"AOD tagging "<<isPrimary<<" "<<isSecondaryWeak<<isSecondaryMaterial<<" "<<partMC->GetMCProcessCode()<<endl;

				  	if(!isPrimary&&!isSecondaryWeak&&!isSecondaryMaterial)//old tagging for old AODs 
				  	{
						AliError("old tagging");
						Int_t mfl=-999,codemoth=-999;
						Int_t indexMoth=partMC->GetMother(); // FIXME ignore fakes? TO BE CHECKED, on ESD is GetFirstMother()
						if(indexMoth>=0)
						{//is not fake
					  		AliAODMCParticle* moth = (AliAODMCParticle*) arrayMC->At(indexMoth);
					  		codemoth = TMath::Abs(moth->GetPdgCode());
					  		mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));
						}
						//Int_t uniqueID = partMC->GetUniqueID();
						//cout<<"uniqueID: "<<partMC->GetUniqueID()<<"       "<<kPDecay<<endl;
						//cout<<"status: "<<partMC->GetStatus()<<"       "<<kPDecay<<endl;
						// if(uniqueID == kPDecay)Printf("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
						if(mfl==3) 
							isSecondaryWeak     = kTRUE; // add if(partMC->GetStatus() & kPDecay)? FIXME
						else       
							isSecondaryMaterial = kTRUE;
				  	}
					//cout<<"AOD 2 tagging "<<isPrimary<<" "<<isSecondaryWeak<<isSecondaryMaterial<<" "<<partMC->GetMCProcessCode()<<endl;
					if(isSecondaryWeak)
					{	
						Int_t indexMoth=partMC->GetMother(); // FIXME ignore fakes? TO BE CHECKED, on ESD is GetFirstMother()
						if(indexMoth>=0)
						{
					  		AliAODMCParticle* moth = (AliAODMCParticle*) arrayMC->At(indexMoth);
							if(!moth)
					  			motherpdg=TMath::Abs(moth->GetPdgCode());
						}
					}

				  	idGen     = fPID->GetParticleSpecie(partMC);
				  	pdgcode=partMC->GetPdgCode(); 
				}
				else if (ifAODEvent==AliSpectraBothTrackCuts::kESDobject)
				{
					TParticle *partMC =stack->Particle(TMath::Abs(track->GetLabel()));
					if (!partMC) 
					{ 
						AliError("Cannot get MC particle");
						continue; 
				  	}
				  	isPrimary           = stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()));
					isSecondaryWeak     = stack->IsSecondaryFromWeakDecay(TMath::Abs(track->GetLabel()));
					isSecondaryMaterial      = stack->IsSecondaryFromMaterial(TMath::Abs(track->GetLabel()));
					//cout<<"ESD tagging "<<isPrimary<<" "<<isSecondaryWeak<<isSecondaryMaterial<<endl;
					
					if(isSecondaryWeak)	
					{

						TParticle* moth=stack->Particle(TMath::Abs(partMC->GetFirstMother()));
						if(moth)
							 motherpdg = TMath::Abs(moth->GetPdgCode());

					}
					

				   	idGen     = fPID->GetParticleSpecie(partMC);
				   	pdgcode=partMC->GetPdgCode(); 
				}
				else
					return;
			
			//	  cout<<isPrimary<<" "<<isSecondaryWeak<<" "<<isSecondaryMaterial<<endl;
			//	  cout<<" functions "<<partMC->IsPhysicalPrimary()<<" "<<partMC->IsSecondaryFromWeakDecay()<<" "<<partMC->IsSecondaryFromMaterial()<<endl;
		  
		  		if (isPrimary&&irec==kSpPion)
					fHistMan->GetPtHistogram(kHistPtRecPrimaryAll)->Fill(track->Pt(),dca);  // PT histo of reconstrutsed primaries in defined eta
		  
		  	//nsigma cut (reconstructed nsigma)
		  		if(idRec == kSpUndefined) 
					continue;
		  
		  	// rapidity cut (reconstructed pt and identity)
		 		 if(!fTrackCuts->CheckYCut ((BothParticleSpecies_t)idRec)) continue;
		  
		  // Get true ID
		  
		  
		 		 if (idRec == idGen) fHistMan->GetHistogram2D(kHistPtRecTrue,  idGen, charge)->Fill(track->Pt(),dca); 
		  
		  		if (isPrimary) 
				{
					fHistMan->GetHistogram2D(kHistPtRecSigmaPrimary, idRec, charge)->Fill(track->Pt(),dca); 
					if(idGen != kSpUndefined) 
					{
		  				fHistMan->GetHistogram2D(kHistPtRecPrimary,      idGen, charge)->Fill(track->Pt(),dca);
		  				if (idRec == idGen) 
							fHistMan->GetHistogram2D(kHistPtRecTruePrimary,  idGen, charge)->Fill(track->Pt(),dca); 
					}
		  		}
		 		 //25th Apr - Muons are added to Pions -- FIXME
		 		if ( pdgcode == 13 && idRec == kSpPion) 
				{ 
					fHistMan->GetPtHistogram(kHistPtRecTrueMuonPlus)->Fill(track->Pt(),dca); 
					if(isPrimary)
		  				fHistMan->GetPtHistogram(kHistPtRecTruePrimaryMuonPlus)->Fill(track->Pt(),dca); 
		  		}
		  		if ( pdgcode == -13 && idRec == kSpPion) 
				{ 
					fHistMan->GetPtHistogram(kHistPtRecTrueMuonMinus)->Fill(track->Pt(),dca); 
					if (isPrimary) 
					{
		  				fHistMan->GetPtHistogram(kHistPtRecTruePrimaryMuonMinus)->Fill(track->Pt(),dca); 
					}
		  		}
		  
		 		 ///..... END FIXME
		  
		  		// Fill secondaries
		  		if(isSecondaryWeak)
				{  
					if(fMotherWDPDGcode>0)
					{
						if(motherpdg==fMotherWDPDGcode)
							fHistMan->GetHistogram2D(kHistPtRecSigmaSecondaryWeakDecay, idRec, charge)->Fill(track->Pt(),dca);
					}	
					else	
						fHistMan->GetHistogram2D(kHistPtRecSigmaSecondaryWeakDecay, idRec, charge)->Fill(track->Pt(),dca);
				}
		  		if(isSecondaryMaterial)  
					fHistMan->GetHistogram2D(kHistPtRecSigmaSecondaryMaterial , idRec, charge)->Fill(track->Pt(),dca);
		  
			}//end if(arrayMC)
		}
	
	
  	} // end loop on tracks

 // cout<< ntracks<<endl;
  fHistMan->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->Fill(npar,ntracks);
  PostData(1, fHistMan  );
  PostData(2, fEventCuts);
  PostData(3, fTrackCuts);
  PostData(4, fPID      );
}

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