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

#include "AliSpectraBothPID.h"
#include "AliAODEvent.h"      
#include "TH1F.h"             
#include "TH2F.h"             
#include "TList.h"            
#include "AliAODTrack.h"      
#include "AliAODMCParticle.h" 
#include "AliPIDResponse.h"   
#include "AliAnalysisManager.h"
#include "AliInputEventHandler.h"
#include "AliSpectraBothTrackCuts.h"

ClassImp(AliSpectraBothPID)

AliSpectraBothPID::AliSpectraBothPID() : TNamed("PID", "PID object"), fPIDType(kNSigmacircleTPCTOF), fNSigmaPID(3), fPIDResponse(0),fshiftTPC(0),fshiftTOF(0) 
{

}

AliSpectraBothPID::AliSpectraBothPID(BothPIDType_t pidType) : TNamed("PID", "PID object"), fPIDType(pidType), fNSigmaPID(3), fPIDResponse(0), fshiftTPC(0),fshiftTOF(0)
{


}


void AliSpectraBothPID::FillQAHistos(AliSpectraBothHistoManager * hman, AliVTrack * track, AliSpectraBothTrackCuts * trackCuts) 
{

  // fill a bunch of QA histos


  // Get PID response object, if needed
  	if(!fPIDResponse) 
  	{
    		AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
   		 AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
    		fPIDResponse = inputHandler->GetPIDResponse();
  	}
	Bool_t ycut[3]={false,false,false};
   	ycut[0]=trackCuts->CheckYCut ((BothParticleSpecies_t) kSpPion);
  	ycut[1]=trackCuts->CheckYCut ((BothParticleSpecies_t) kSpKaon);
	ycut[2]=trackCuts->CheckYCut ((BothParticleSpecies_t) kSpProton);
 
 
  //Response
  	AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(track);
  
  	hman->GetPIDHistogram(kHistPIDTPC)->Fill(track->GetTPCmomentum(), track->GetTPCsignal()*track->Charge()); // PID histo
  	hman->GetPIDHistogram(kHistPIDTPCPion)->Fill(track->GetTPCmomentum(),fPIDResponse->GetTPCResponse().GetExpectedSignal(track->GetTPCmomentum(),AliPID::kPion)*track->Charge()); // Expected PIDPion histo
  	hman->GetPIDHistogram(kHistPIDTPCKaon)->Fill(track->GetTPCmomentum(),fPIDResponse->GetTPCResponse().GetExpectedSignal(track->GetTPCmomentum(),AliPID::kKaon)*track->Charge()); // Expected PIDKaon histo
  	hman->GetPIDHistogram(kHistPIDTPCProton)->Fill(track->GetTPCmomentum(),fPIDResponse->GetTPCResponse().GetExpectedSignal(track->GetTPCmomentum(),AliPID::kProton)*track->Charge()); // Expected PIDProton histo
  
  	Double_t nsigmaTPCkProton = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton)+fshiftTPC;
  	Double_t nsigmaTPCkKaon = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon)+fshiftTPC; 
  	Double_t nsigmaTPCkPion = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion)+fshiftTPC; 
  	Double_t nsigmaTOFkProton=100.0,nsigmaTOFkKaon=100.0,nsigmaTOFkPion=100.0;

  	Double_t nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
  	Double_t nsigmaTPCTOFkKaon   = TMath::Abs(nsigmaTPCkKaon);
  	Double_t nsigmaTPCTOFkPion   = TMath::Abs(nsigmaTPCkPion);
	if(fPIDType==kNSigmaTOF)
	{
		nsigmaTPCTOFkProton = 100.0;
  		nsigmaTPCTOFkKaon   = 100.0;
  		nsigmaTPCTOFkPion   =100.0;

	}

 	if(track->Pt()>trackCuts->GetPtTOFMatching())
  	{
    
   		 hman->GetPIDHistogram(kHistPIDTOF)->Fill(track->P(),(track->GetTOFsignal()/100)*track->Charge()); // PID histo
    
   		 nsigmaTOFkProton = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton)+fshiftTOF;
  		 nsigmaTOFkKaon = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon)+fshiftTOF; 
   		 nsigmaTOFkPion = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion)+fshiftTOF; 
    
    //TOF	 	 
		if(ycut[0])
    		{
    			hman->GetPtHistogram(kHistNSigPionTOF)->Fill(track->P(),nsigmaTOFkPion );
    			hman->GetPtHistogram(kHistNSigPionPtTOF)->Fill(track->Pt(),nsigmaTOFkPion );
    		}
    		if(ycut[1])
  	 	{	
    			hman->GetPtHistogram(kHistNSigKaonTOF)->Fill(track->P(),nsigmaTOFkKaon);
   	 		hman->GetPtHistogram(kHistNSigKaonPtTOF)->Fill(track->Pt(),nsigmaTOFkKaon );	 
   		}
		if(ycut[2])
		{					
   			hman->GetPtHistogram(kHistNSigProtonTOF)->Fill(track->P(),nsigmaTOFkProton ); 
    			hman->GetPtHistogram(kHistNSigProtonPtTOF)->Fill(track->Pt(),nsigmaTOFkProton );
		}
	
    
   		 nsigmaTPCTOFkProton = TMath::Sqrt((nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton)/2.0);
    		nsigmaTPCTOFkKaon = TMath::Sqrt((nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon)/2.0);
    		nsigmaTPCTOFkPion = TMath::Sqrt((nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion)/2.0);
  
  	}

   	if(ycut[0])
    	{
    		hman->GetPtHistogram(kHistNSigPionTPC)->Fill(track->P(),nsigmaTPCkPion );
    		hman->GetPtHistogram(kHistNSigPionPtTPC)->Fill(track->Pt(),nsigmaTPCkPion );
  		hman->GetPtHistogram(kHistNSigPionTPCTOF)->Fill(track->P(),nsigmaTPCTOFkPion);
  		hman->GetPtHistogram(kHistNSigPionPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkPion);
	}
    	if(ycut[1])
  	 {	
    		hman->GetPtHistogram(kHistNSigKaonTPC)->Fill(track->P(),nsigmaTPCkKaon );
   	 	hman->GetPtHistogram(kHistNSigKaonPtTPC)->Fill(track->Pt(),nsigmaTPCkKaon );	 
  		hman->GetPtHistogram(kHistNSigKaonTPCTOF)->Fill(track->P(),nsigmaTPCTOFkKaon);
  		hman->GetPtHistogram(kHistNSigKaonPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkKaon);
		
	}
	if(ycut[2])
	{					
   		hman->GetPtHistogram(kHistNSigProtonTPC)->Fill(track->P(),nsigmaTPCkProton ); 
    		hman->GetPtHistogram(kHistNSigProtonPtTPC)->Fill(track->Pt(),nsigmaTPCkProton );
  		hman->GetPtHistogram(kHistNSigProtonTPCTOF)->Fill(track->P(),nsigmaTPCTOFkProton);
  		hman->GetPtHistogram(kHistNSigProtonPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkProton);
	}

}

Int_t AliSpectraBothPID::GetParticleSpecie(AliAODMCParticle * part) 
{
  // return PID according to MC truth
  	switch(TMath::Abs(part->PdgCode()))
	{
  		case 2212:
   		return kSpProton;
   		 break;
  		case 321:
    		return kSpKaon;
    		break;
  		case 211:
    		return kSpPion;
    		break;
  		default:
    		return kSpUndefined;
  	} 
}

Int_t AliSpectraBothPID::GetParticleSpecie(TParticle* part) 
{
  // return PID according to MC truth
  switch(TMath::Abs(part->GetPdgCode()))
  {
  	case 2212:
    	return kSpProton;
   	 break;
  	case 321:
    	return kSpKaon;
    	break;
  	case 211:
    	return kSpPion;
    	break;
  	default:
    	return kSpUndefined;
  } 
}

Int_t AliSpectraBothPID::GetParticleSpecie(AliSpectraBothHistoManager * hman,AliVTrack      * trk, AliSpectraBothTrackCuts * trackCuts, Bool_t* rec) 
{
  // return PID according to detectors
  // Get PID response object, if needed

	  // guess the particle based on the smaller nsigma
  	rec[kSpPion]=false;
 	 rec[kSpKaon]=false;
  	rec[kSpProton]=false;

 	if(!fPIDResponse) 
 	 {
    		AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
    		AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
   		 fPIDResponse = inputHandler->GetPIDResponse();
  	}

  	if(!fPIDResponse) 
	{
    		AliFatal("Cannot get pid response");
    	return 0;
  	}


  // Compute nsigma for each hypthesis
  	AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(trk);

  // --- TPC
  	Double_t nsigmaTPCkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton)+fshiftTPC);
  	Double_t nsigmaTPCkKaon   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon)+fshiftTPC); 
 	 Double_t nsigmaTPCkPion   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion)+fshiftTPC);
  
 
	Double_t nsigmaTOFkProton=0.0,nsigmaTOFkKaon=0.0,nsigmaTOFkPion=0.0;


  	Double_t nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
  	Double_t nsigmaTPCTOFkKaon   = TMath::Abs(nsigmaTPCkKaon);
  	Double_t nsigmaTPCTOFkPion   = TMath::Abs(nsigmaTPCkPion);
	if(fPIDType==kNSigmaTOF)
	{
		nsigmaTPCTOFkProton = 100.0;
  		nsigmaTPCTOFkKaon   = 100.0;
  		nsigmaTPCTOFkPion   =100.0;
		nsigmaTOFkProton = 100.0;
  		nsigmaTOFkKaon   = 100.0;
  		nsigmaTOFkPion   =100.0;

	}
	

  	if(trk->Pt()>trackCuts->GetPtTOFMatching())
  	{
    		nsigmaTOFkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton)+fshiftTOF);
   	 	nsigmaTOFkKaon   = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon)+fshiftTOF); 
    		nsigmaTOFkPion   = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion)+fshiftTOF); 
    		// the TOF info is used in combined
    		nsigmaTPCTOFkProton = TMath::Sqrt((nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton)/2.);
		nsigmaTPCTOFkKaon   = TMath::Sqrt((nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon)/2.);
		nsigmaTPCTOFkPion   = TMath::Sqrt((nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion)/2.);
  	}
	 else
	{ 
		if (fPIDType==kNSigmaTOF)
			return kSpUndefined;

	}
	
  	// select the nsigma to be used for the actual PID
  	Double_t nsigmaPion=100; 
	Double_t nsigmaKaon=100; 
	Double_t nsigmaProton=100;

  	switch (fPIDType) 
	{
  		case kNSigmaTPC:
   		nsigmaProton  =  nsigmaTPCkProton;
    		nsigmaKaon	  =  nsigmaTPCkKaon  ;
    		nsigmaPion    =  nsigmaTPCkPion  ;
    		break;
  		case kNSigmaTOF:
    		nsigmaProton  =  nsigmaTOFkProton;
    		nsigmaKaon	  =  nsigmaTOFkKaon  ;
    		nsigmaPion    =  nsigmaTOFkPion  ;
    		break;
  		case kNSigmacircleTPCTOF:
    		nsigmaProton  =  nsigmaTPCTOFkProton;
    		nsigmaKaon	  =  nsigmaTPCTOFkKaon  ;
    		nsigmaPion    =  nsigmaTPCTOFkPion  ;
    		break;
		case kNSigmasquareTPCTOF:
		nsigmaProton  =  TMath::Max(nsigmaTPCkProton,nsigmaTOFkProton);
    		nsigmaKaon	  =  TMath::Max(nsigmaTPCkKaon,nsigmaTOFkKaon);
    		nsigmaPion    =  TMath::Max(nsigmaTPCkPion,nsigmaTOFkPion)  ;
    		break;

  	}	
  
  	if(nsigmaPion   < fNSigmaPID)
		rec[kSpPion]=true;
  	if(nsigmaKaon   < fNSigmaPID)
		rec[kSpKaon]=true;
 	 if(nsigmaProton   < fNSigmaPID)
		rec[kSpProton]=true;
	
  	if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) 
		return kSpUndefined;//if is the default value for the three
  	if( ( nsigmaKaon   < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon   < fNSigmaPID))
	{
		if(hman->GetPIDHistogram(kHistPIDTPCKaonRec))
    			hman->GetPIDHistogram(kHistPIDTPCKaonRec)->Fill(trk->GetTPCmomentum(), trk->GetTPCsignal()*trk->Charge()); // Reconstructed PIDKaon histo
    		return kSpKaon;
  	}
  	if( ( nsigmaPion   < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion   < fNSigmaPID))
	{
		if(hman->GetPIDHistogram(kHistPIDTPCPionRec))
    			hman->GetPIDHistogram(kHistPIDTPCPionRec)->Fill(trk->GetTPCmomentum(), trk->GetTPCsignal()*trk->Charge()); // Reconstructed PIDPion histo
    		return kSpPion;
  	}
  	if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID))
	{
		if(hman->GetPIDHistogram(kHistPIDTPCProtonRec))
    			hman->GetPIDHistogram(kHistPIDTPCProtonRec)->Fill(trk->GetTPCmomentum(), trk->GetTPCsignal()*trk->Charge()); // Reconstructed PIDProton histo
    		return kSpProton;
  	}
  	// else, return undefined
  	return kSpUndefined;

}

Long64_t AliSpectraBothPID::Merge(TCollection* list)
{
  // Merging interface.
  // Returns the number of merged objects (including this).

  Printf("Merging");

  if (!list)
    return 0;

  if (list->IsEmpty())
    return 1;

  TIterator* iter = list->MakeIterator();
  TObject* obj;

  // Actually, we don't do anything here...
  // collections of all histograms
  //  TList collections;

  Int_t count = 0;

  while ((obj = iter->Next())) {
    AliSpectraBothPID* entry = dynamic_cast<AliSpectraBothPID*> (obj);
    if (entry == 0) 
      continue;

    // TH1I * histo = entry->GetHistoCuts();      
    // collections.Add(histo);
    count++;
  }
  
  //  fHistoCuts->Merge(&collections);
  
  delete iter;
  Printf("OK");
  return count+1;
}
Commit	Line	Data
239a080a	1	#include "AliSpectraBothPID.h"
	2	#include "AliAODEvent.h"
	3	#include "TH1F.h"
	4	#include "TH2F.h"
	5	#include "TList.h"
	6	#include "AliAODTrack.h"
	7	#include "AliAODMCParticle.h"
	8	#include "AliPIDResponse.h"
	9	#include "AliAnalysisManager.h"
	10	#include "AliInputEventHandler.h"
	11	#include "AliSpectraBothTrackCuts.h"
	12
	13	ClassImp(AliSpectraBothPID)
	14
2d98dd91	15	AliSpectraBothPID::AliSpectraBothPID() : TNamed("PID", "PID object"), fPIDType(kNSigmacircleTPCTOF), fNSigmaPID(3), fPIDResponse(0),fshiftTPC(0),fshiftTOF(0)
239a080a	16	{
	17
	18	}
	19
	20	AliSpectraBothPID::AliSpectraBothPID(BothPIDType_t pidType) : TNamed("PID", "PID object"), fPIDType(pidType), fNSigmaPID(3), fPIDResponse(0), fshiftTPC(0),fshiftTOF(0)
	21	{
	22
	23
	24
	25	}
	26
	27
	28
	29	void AliSpectraBothPID::FillQAHistos(AliSpectraBothHistoManager * hman, AliVTrack * track, AliSpectraBothTrackCuts * trackCuts)
	30	{
	31
	32	// fill a bunch of QA histos
	33
	34
	35	// Get PID response object, if needed
	36	if(!fPIDResponse)
	37	{
	38	AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
	39	AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
	40	fPIDResponse = inputHandler->GetPIDResponse();
	41	}
b3ea73e1	42	Bool_t ycut[3]={false,false,false};
	43	ycut[0]=trackCuts->CheckYCut ((BothParticleSpecies_t) kSpPion);
	44	ycut[1]=trackCuts->CheckYCut ((BothParticleSpecies_t) kSpKaon);
	45	ycut[2]=trackCuts->CheckYCut ((BothParticleSpecies_t) kSpProton);
	46
	47
239a080a	48	//Response
	49	AliVParticle inEvHMain = dynamic_cast<AliVParticle >(track);
	50
	51	hman->GetPIDHistogram(kHistPIDTPC)->Fill(track->GetTPCmomentum(), track->GetTPCsignal()*track->Charge()); // PID histo
	52	hman->GetPIDHistogram(kHistPIDTPCPion)->Fill(track->GetTPCmomentum(),fPIDResponse->GetTPCResponse().GetExpectedSignal(track->GetTPCmomentum(),AliPID::kPion)*track->Charge()); // Expected PIDPion histo
	53	hman->GetPIDHistogram(kHistPIDTPCKaon)->Fill(track->GetTPCmomentum(),fPIDResponse->GetTPCResponse().GetExpectedSignal(track->GetTPCmomentum(),AliPID::kKaon)*track->Charge()); // Expected PIDKaon histo
	54	hman->GetPIDHistogram(kHistPIDTPCProton)->Fill(track->GetTPCmomentum(),fPIDResponse->GetTPCResponse().GetExpectedSignal(track->GetTPCmomentum(),AliPID::kProton)*track->Charge()); // Expected PIDProton histo
	55
b3ea73e1	56	Double_t nsigmaTPCkProton = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton)+fshiftTPC;
	57	Double_t nsigmaTPCkKaon = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon)+fshiftTPC;
	58	Double_t nsigmaTPCkPion = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion)+fshiftTPC;
239a080a	59	Double_t nsigmaTOFkProton=100.0,nsigmaTOFkKaon=100.0,nsigmaTOFkPion=100.0;
	60
	61	Double_t nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
	62	Double_t nsigmaTPCTOFkKaon = TMath::Abs(nsigmaTPCkKaon);
	63	Double_t nsigmaTPCTOFkPion = TMath::Abs(nsigmaTPCkPion);
	64	if(fPIDType==kNSigmaTOF)
	65	{
	66	nsigmaTPCTOFkProton = 100.0;
	67	nsigmaTPCTOFkKaon = 100.0;
	68	nsigmaTPCTOFkPion =100.0;
	69
	70	}
	71
	72	if(track->Pt()>trackCuts->GetPtTOFMatching())
	73	{
	74
	75	hman->GetPIDHistogram(kHistPIDTOF)->Fill(track->P(),(track->GetTOFsignal()/100)*track->Charge()); // PID histo
	76
b3ea73e1	77	nsigmaTOFkProton = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton)+fshiftTOF;
	78	nsigmaTOFkKaon = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon)+fshiftTOF;
	79	nsigmaTOFkPion = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion)+fshiftTOF;
239a080a	80
b3ea73e1	81	//TOF
	82	if(ycut[0])
	83	{
	84	hman->GetPtHistogram(kHistNSigPionTOF)->Fill(track->P(),nsigmaTOFkPion );
	85	hman->GetPtHistogram(kHistNSigPionPtTOF)->Fill(track->Pt(),nsigmaTOFkPion );
	86	}
	87	if(ycut[1])
	88	{
	89	hman->GetPtHistogram(kHistNSigKaonTOF)->Fill(track->P(),nsigmaTOFkKaon);
	90	hman->GetPtHistogram(kHistNSigKaonPtTOF)->Fill(track->Pt(),nsigmaTOFkKaon );
	91	}
	92	if(ycut[2])
	93	{
	94	hman->GetPtHistogram(kHistNSigProtonTOF)->Fill(track->P(),nsigmaTOFkProton );
	95	hman->GetPtHistogram(kHistNSigProtonPtTOF)->Fill(track->Pt(),nsigmaTOFkProton );
	96	}
	97
239a080a	98
2d98dd91	99	nsigmaTPCTOFkProton = TMath::Sqrt((nsigmaTPCkProtonnsigmaTPCkProton+nsigmaTOFkProtonnsigmaTOFkProton)/2.0);
	100	nsigmaTPCTOFkKaon = TMath::Sqrt((nsigmaTPCkKaonnsigmaTPCkKaon+nsigmaTOFkKaonnsigmaTOFkKaon)/2.0);
	101	nsigmaTPCTOFkPion = TMath::Sqrt((nsigmaTPCkPionnsigmaTPCkPion+nsigmaTOFkPionnsigmaTOFkPion)/2.0);
239a080a	102
239a080a	103	}
239a080a	104
b3ea73e1	105	if(ycut[0])
	106	{
	107	hman->GetPtHistogram(kHistNSigPionTPC)->Fill(track->P(),nsigmaTPCkPion );
	108	hman->GetPtHistogram(kHistNSigPionPtTPC)->Fill(track->Pt(),nsigmaTPCkPion );
	109	hman->GetPtHistogram(kHistNSigPionTPCTOF)->Fill(track->P(),nsigmaTPCTOFkPion);
	110	hman->GetPtHistogram(kHistNSigPionPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkPion);
	111	}
	112	if(ycut[1])
	113	{
	114	hman->GetPtHistogram(kHistNSigKaonTPC)->Fill(track->P(),nsigmaTPCkKaon );
	115	hman->GetPtHistogram(kHistNSigKaonPtTPC)->Fill(track->Pt(),nsigmaTPCkKaon );
	116	hman->GetPtHistogram(kHistNSigKaonTPCTOF)->Fill(track->P(),nsigmaTPCTOFkKaon);
	117	hman->GetPtHistogram(kHistNSigKaonPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkKaon);
	118
	119	}
	120	if(ycut[2])
	121	{
	122	hman->GetPtHistogram(kHistNSigProtonTPC)->Fill(track->P(),nsigmaTPCkProton );
	123	hman->GetPtHistogram(kHistNSigProtonPtTPC)->Fill(track->Pt(),nsigmaTPCkProton );
	124	hman->GetPtHistogram(kHistNSigProtonTPCTOF)->Fill(track->P(),nsigmaTPCTOFkProton);
	125	hman->GetPtHistogram(kHistNSigProtonPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkProton);
	126	}
239a080a	127
	128	}
	129
	130	Int_t AliSpectraBothPID::GetParticleSpecie(AliAODMCParticle * part)
	131	{
	132	// return PID according to MC truth
	133	switch(TMath::Abs(part->PdgCode()))
	134	{
	135	case 2212:
	136	return kSpProton;
	137	break;
	138	case 321:
	139	return kSpKaon;
	140	break;
	141	case 211:
	142	return kSpPion;
	143	break;
	144	default:
	145	return kSpUndefined;
	146	}
	147	}
	148
	149	Int_t AliSpectraBothPID::GetParticleSpecie(TParticle* part)
	150	{
	151	// return PID according to MC truth
	152	switch(TMath::Abs(part->GetPdgCode()))
	153	{
	154	case 2212:
	155	return kSpProton;
	156	break;
	157	case 321:
	158	return kSpKaon;
	159	break;
	160	case 211:
	161	return kSpPion;
	162	break;
	163	default:
	164	return kSpUndefined;
	165	}
	166	}
	167
	168	Int_t AliSpectraBothPID::GetParticleSpecie(AliSpectraBothHistoManager * hman,AliVTrack * trk, AliSpectraBothTrackCuts * trackCuts, Bool_t* rec)
	169	{
	170	// return PID according to detectors
	171	// Get PID response object, if needed
	172
	173	// guess the particle based on the smaller nsigma
239a080a	174	rec[kSpPion]=false;
	175	rec[kSpKaon]=false;
	176	rec[kSpProton]=false;
	177
	178	if(!fPIDResponse)
	179	{
	180	AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
	181	AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
	182	fPIDResponse = inputHandler->GetPIDResponse();
	183	}
	184
	185	if(!fPIDResponse)
	186	{
	187	AliFatal("Cannot get pid response");
	188	return 0;
	189	}
	190
	191
	192	// Compute nsigma for each hypthesis
	193	AliVParticle inEvHMain = dynamic_cast<AliVParticle >(trk);
	194
	195	// --- TPC
	196	Double_t nsigmaTPCkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton)+fshiftTPC);
	197	Double_t nsigmaTPCkKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon)+fshiftTPC);
	198	Double_t nsigmaTPCkPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion)+fshiftTPC);
	199
	200
	201
2d98dd91	202	Double_t nsigmaTOFkProton=0.0,nsigmaTOFkKaon=0.0,nsigmaTOFkPion=0.0;
	203
	204
239a080a	205
	206	Double_t nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
	207	Double_t nsigmaTPCTOFkKaon = TMath::Abs(nsigmaTPCkKaon);
	208	Double_t nsigmaTPCTOFkPion = TMath::Abs(nsigmaTPCkPion);
	209	if(fPIDType==kNSigmaTOF)
	210	{
	211	nsigmaTPCTOFkProton = 100.0;
	212	nsigmaTPCTOFkKaon = 100.0;
	213	nsigmaTPCTOFkPion =100.0;
2d98dd91	214	nsigmaTOFkProton = 100.0;
	215	nsigmaTOFkKaon = 100.0;
	216	nsigmaTOFkPion =100.0;
239a080a	217
	218	}
	219
	220
	221
	222	if(trk->Pt()>trackCuts->GetPtTOFMatching())
	223	{
	224	nsigmaTOFkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton)+fshiftTOF);
	225	nsigmaTOFkKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon)+fshiftTOF);
	226	nsigmaTOFkPion = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion)+fshiftTOF);
	227	// the TOF info is used in combined
	228	nsigmaTPCTOFkProton = TMath::Sqrt((nsigmaTPCkProtonnsigmaTPCkProton+nsigmaTOFkProtonnsigmaTOFkProton)/2.);
	229	nsigmaTPCTOFkKaon = TMath::Sqrt((nsigmaTPCkKaonnsigmaTPCkKaon+nsigmaTOFkKaonnsigmaTOFkKaon)/2.);
	230	nsigmaTPCTOFkPion = TMath::Sqrt((nsigmaTPCkPionnsigmaTPCkPion+nsigmaTOFkPionnsigmaTOFkPion)/2.);
	231	}
	232	else
	233	{
	234	if (fPIDType==kNSigmaTOF)
	235	return kSpUndefined;
	236
	237	}
	238
	239	// select the nsigma to be used for the actual PID
	240	Double_t nsigmaPion=100;
	241	Double_t nsigmaKaon=100;
	242	Double_t nsigmaProton=100;
	243
	244	switch (fPIDType)
	245	{
	246	case kNSigmaTPC:
	247	nsigmaProton = nsigmaTPCkProton;
	248	nsigmaKaon = nsigmaTPCkKaon ;
	249	nsigmaPion = nsigmaTPCkPion ;
	250	break;
	251	case kNSigmaTOF:
	252	nsigmaProton = nsigmaTOFkProton;
	253	nsigmaKaon = nsigmaTOFkKaon ;
	254	nsigmaPion = nsigmaTOFkPion ;
	255	break;
2d98dd91	256	case kNSigmacircleTPCTOF:
239a080a	257	nsigmaProton = nsigmaTPCTOFkProton;
	258	nsigmaKaon = nsigmaTPCTOFkKaon ;
	259	nsigmaPion = nsigmaTPCTOFkPion ;
	260	break;
2d98dd91	261	case kNSigmasquareTPCTOF:
	262	nsigmaProton = TMath::Max(nsigmaTPCkProton,nsigmaTOFkProton);
	263	nsigmaKaon = TMath::Max(nsigmaTPCkKaon,nsigmaTOFkKaon);
	264	nsigmaPion = TMath::Max(nsigmaTPCkPion,nsigmaTOFkPion) ;
	265	break;
	266
239a080a	267	}
	268
	269	if(nsigmaPion < fNSigmaPID)
	270	rec[kSpPion]=true;
	271	if(nsigmaKaon < fNSigmaPID)
	272	rec[kSpKaon]=true;
	273	if(nsigmaProton < fNSigmaPID)
	274	rec[kSpProton]=true;
	275
	276	if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton ))
	277	return kSpUndefined;//if is the default value for the three
	278	if( ( nsigmaKaon < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon < fNSigmaPID))
	279	{
be25efef	280	if(hman->GetPIDHistogram(kHistPIDTPCKaonRec))
be25efef	281	hman->GetPIDHistogram(kHistPIDTPCKaonRec)->Fill(trk->GetTPCmomentum(), trk->GetTPCsignal()*trk->Charge()); // Reconstructed PIDKaon histo
239a080a	282	return kSpKaon;
	283	}
	284	if( ( nsigmaPion < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion < fNSigmaPID))
	285	{
be25efef	286	if(hman->GetPIDHistogram(kHistPIDTPCPionRec))
be25efef	287	hman->GetPIDHistogram(kHistPIDTPCPionRec)->Fill(trk->GetTPCmomentum(), trk->GetTPCsignal()*trk->Charge()); // Reconstructed PIDPion histo
239a080a	288	return kSpPion;
	289	}
	290	if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID))
	291	{
be25efef	292	if(hman->GetPIDHistogram(kHistPIDTPCProtonRec))
be25efef	293	hman->GetPIDHistogram(kHistPIDTPCProtonRec)->Fill(trk->GetTPCmomentum(), trk->GetTPCsignal()*trk->Charge()); // Reconstructed PIDProton histo
239a080a	294	return kSpProton;
	295	}
	296	// else, return undefined
	297	return kSpUndefined;
	298
	299	}
	300
	301	Long64_t AliSpectraBothPID::Merge(TCollection* list)
	302	{
	303	// Merging interface.
	304	// Returns the number of merged objects (including this).
	305
	306	Printf("Merging");
	307
	308	if (!list)
	309	return 0;
	310
	311	if (list->IsEmpty())
	312	return 1;
	313
	314	TIterator* iter = list->MakeIterator();
	315	TObject* obj;
	316
	317	// Actually, we don't do anything here...
	318	// collections of all histograms
	319	// TList collections;
	320
	321	Int_t count = 0;
	322
	323	while ((obj = iter->Next())) {
	324	AliSpectraBothPID* entry = dynamic_cast<AliSpectraBothPID*> (obj);
	325	if (entry == 0)
	326	continue;
	327
	328	// TH1I * histo = entry->GetHistoCuts();
	329	// collections.Add(histo);
	330	count++;
	331	}
	332
	333	// fHistoCuts->Merge(&collections);
	334
	335	delete iter;
	336	Printf("OK");
	337	return count+1;
	338	}
	339