[u/mrichter/AliRoot.git] / ANALYSIS / ANALYSISalice / AliAnalysisTaskPIDCombined.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.                  *
**************************************************************************/

///////////////////////////////////////////////////////////////////////////
//                                                                       //
//                        Basic Analysis Task                            //
//                      for PID        Analysis                          //
//                                                                       //
///////////////////////////////////////////////////////////////////////////

#include <TChain.h>
#include <TH1D.h>
#include <TH2D.h>

#include <AliVEvent.h>
#include <AliInputEventHandler.h>
#include <AliAnalysisManager.h>

#include <AliLog.h>
#include <AliPID.h>
#include <AliPIDCombined.h>
#include <AliPIDResponse.h>

#include "AliAnalysisTaskPIDCombined.h"

const char *AliAnalysisTaskPIDCombined::fgkBinMomDesc[AliAnalysisTaskPIDCombined::kPtBins] = {
  " 0 <= p < 0.5 GeV/c",
  " 0.5 <= p < 0.7 GeV/c",
  " 0.7 <= p < 1.0 GeV/c",
  " 1.0 <= p < 1.5 GeV/c",
  " 1.5 <= p < 2.0 GeV/c",
  " p >= 2.0 GeV/c"
};

ClassImp(AliAnalysisTaskPIDCombined)

//_________________________________________________________________________________
AliAnalysisTaskPIDCombined::AliAnalysisTaskPIDCombined() :
  AliAnalysisTaskSE(),
  fHistList(),
  fProbTPCnSigma(),
  fProbTOFnSigma(),
  fProbTPCTOFnSigmaTPC(),
  fProbTPC(),
  fProbTOF(),
  fProbTPCTOF(),
  fPriors(),
  fProbTPCTOFnSigTPCMom(),
  fProbTPCnSigTPCMom(),
  fProbTOFnSigTOFMom(),
  fPriorsUsed(),
  fPIDResponse(0x0),
  fPIDCombined(0x0),
  fTrackCuts(0x0),
  fTrackFilter(0x0),
  fDeDx(NULL),
  fDeDxTuned(NULL)
{
  //
  // Constructor
  //
}

//_________________________________________________________________________________
AliAnalysisTaskPIDCombined::AliAnalysisTaskPIDCombined(const char *name) :
  AliAnalysisTaskSE(name),
  fHistList(),
  fProbTPCnSigma(),
  fProbTOFnSigma(),
  fProbTPCTOFnSigmaTPC(),
  fProbTPC(),
  fProbTOF(),
  fProbTPCTOF(),
  fPriors(),
  fProbTPCTOFnSigTPCMom(),
  fProbTPCnSigTPCMom(),
  fProbTOFnSigTOFMom(),
  fPriorsUsed(),
  fPIDResponse(0x0),
  fPIDCombined(0x0),
  fTrackCuts(0x0),
  fTrackFilter(0x0),
  fDeDx(NULL),
  fDeDxTuned(NULL)
{
  //
  // Constructor
  //
  DefineInput(0,TChain::Class());
  DefineOutput(1, TList::Class());
}

//_________________________________________________________________________________
void AliAnalysisTaskPIDCombined::UserCreateOutputObjects()
{
  //
  // Initialise the framework objects
  //


  // ------- track cuts
  fTrackCuts = new AliESDtrackCuts("fTrackCuts", "Standard");
  fTrackCuts->SetAcceptKinkDaughters(kFALSE);
  fTrackCuts->SetMinNClustersTPC(80);
  fTrackCuts->SetMaxChi2PerClusterTPC(4);
  fTrackCuts->SetMaxDCAToVertexXY(3);
  fTrackCuts->SetMaxDCAToVertexZ(3);
  fTrackCuts->SetRequireTPCRefit(kTRUE);
  fTrackCuts->SetRequireITSRefit(kTRUE);
  fTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
  fTrackFilter = new AliAnalysisFilter("trackFilter");
  fTrackFilter->AddCuts(fTrackCuts);


  // ------- setup PIDCombined
  fPIDCombined=new AliPIDCombined;
  fPIDCombined->SetDefaultTPCPriors();
  fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);

  // no light nuclei - no need to call it, this is default
  //  fPIDCombined->SetSelectedSpecies(AliPID::kSPECIES);


  fHistList.Add(new TH1D("nEvents","Number of Evnets;Selection",2,0,2));

  for (Int_t ispec=0; ispec<5; ++ispec){

 
    fProbTPC[ispec]=new TH2D(Form("prob%s_mom_TPC",AliPID::ParticleName(ispec)),
                                   Form("%s probability vs. momentum;momentum;probability",AliPID::ParticleName(ispec)),
                                   100,0.,20.,50,0.,1.);
    fHistList.Add(fProbTPC[ispec]);
    fProbTPCnSigma[ispec]=new TH2D(Form("prob%s_nSigma_TPC",AliPID::ParticleName(ispec)),
                                   Form("%s probability vs. n#sigma;n#sigma;probability",AliPID::ParticleName(ispec)),
                                   20,-5.,5.,50,0.,1.);
    fHistList.Add(fProbTPCnSigma[ispec]);

    for (Int_t ibin=0;ibin<kPtBins;ibin++) {
      fProbTPCnSigTPCMom[ibin][ispec]=new TH2D(Form("prob%s_nSigma_TPC (%s)",AliPID::ParticleName(ispec),fgkBinMomDesc[ibin]),
					       Form("%s probability vs. n#sigma;n#sigma;probability",AliPID::ParticleName(ispec)),
					       20,-5.,5.,50,0.,1.);
      fHistList.Add(fProbTPCnSigTPCMom[ibin][ispec]);
    }


    fProbTOF[ispec]=new TH2D(Form("prob%s_mom_TOF",AliPID::ParticleName(ispec)),
                                   Form("%s probability vs. momentum;momentum;probability",AliPID::ParticleName(ispec)),
                                   100,0.,20.,50,0.,1.);
    fHistList.Add(fProbTOF[ispec]);
    fProbTOFnSigma[ispec]=new TH2D(Form("prob%s_nSigma_TOF",AliPID::ParticleName(ispec)),
                                   Form("%s probability vs. n#sigma;n#sigma;probability",AliPID::ParticleName(ispec)),
                                   20,-5.,5.,50,0.,1.);
    fHistList.Add(fProbTOFnSigma[ispec]);
    for (Int_t ibin=0;ibin<kPtBins;ibin++) {
      fProbTOFnSigTOFMom[ibin][ispec]=new TH2D(Form("prob%s_nSigma_TOF (%s)",AliPID::ParticleName(ispec),fgkBinMomDesc[ibin]),
					       Form("%s probability vs. n#sigma;n#sigma;probability",AliPID::ParticleName(ispec)),
					       20,-5.,5.,50,0.,1.);
      fHistList.Add(fProbTOFnSigTOFMom[ibin][ispec]);
    }


    fProbTPCTOF[ispec]=new TH2D(Form("prob%s_mom_TPCTOF",AliPID::ParticleName(ispec)),
                                   Form("%s probability vs. momentum;momentum;probability",AliPID::ParticleName(ispec)),
                                   100,0.,20.,50,0.,1.);
    fHistList.Add(fProbTPCTOF[ispec]);
    fProbTPCTOFnSigmaTPC[ispec]=new TH2D(Form("probTPCTOF%s_nSigma_TPC",AliPID::ParticleName(ispec)),
                                   Form("%s TPCTOF probability vs. n#sigmaTPC;n#sigma;probability",AliPID::ParticleName(ispec)),
                                   20,-5.,5.,50,0.,1.);
    fHistList.Add(fProbTPCTOFnSigmaTPC[ispec]);
    for (Int_t ibin=0;ibin<kPtBins;ibin++) {
      fProbTPCTOFnSigTPCMom[ibin][ispec]=new TH2D(Form("probTPCTOF%s_nSigma_TPC (%s)",AliPID::ParticleName(ispec),fgkBinMomDesc[ibin]),
					       Form("%s probability vs. n#sigma;n#sigma;probability",AliPID::ParticleName(ispec)),
					       20,-5.,5.,50,0.,1.);
      fHistList.Add(fProbTPCTOFnSigTPCMom[ibin][ispec]);
    }


    // basic priors
    fPriors[ispec]=new TH1F(Form("%s_priors",AliPID::ParticleName(ispec)),
			    Form("%s priors vs momentum",AliPID::ParticleName(ispec)),
			    100,0.,20.);
    fHistList.Add(fPriors[ispec]);
    switch (ispec) {
    case AliPID::kElectron:
      for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.02);
      break;
    case AliPID::kMuon:
      for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.02);
      break;
    case AliPID::kPion:
      for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.56);
      break;
    case AliPID::kKaon:
      for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.20);
      break;
    case AliPID::kProton:
      for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.20);
      break;
    default:
      break;
    }
    fPIDCombined->SetPriorDistribution((AliPID::EParticleType)ispec,fPriors[ispec]);

    // priors used
    fPriorsUsed[ispec] = new TH2D(Form("%s_priorsUsed",AliPID::ParticleName(ispec)),
			    Form("%s priors vs transverse momentum;p_{t} (GeV/c);priors",AliPID::ParticleName(ispec)),
				  100,0.,20.,101,0,1.01);      
    fHistList.Add(fPriorsUsed[ispec]);
  }


  fDeDx = new TH2D("hDeDx",";p_{TPC};dE/dx (a.u.)",500,0,5,500,0,500);
  fHistList.Add(fDeDx);
  fDeDxTuned = new TH2D("hDeDxTuned",";p_{TPC};dE/dx (a.u.)",500,0,5,500,0,500);
  fHistList.Add(fDeDxTuned);

  fHistList.SetOwner();
  PostData(1,&fHistList);


}

//_________________________________________________________________________________
void AliAnalysisTaskPIDCombined::UserExec(Option_t *)
{
  //
  // Main loop. Called for every event
  //
  AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
  AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
  fPIDResponse=inputHandler->GetPIDResponse();
  if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");

  //  Printf(" ---------------------- UserExec PID task ---------------------");
  
  FillHistogram("nEvents",0.);

  AliVEvent *event=InputEvent();
  AliVTrack *track=0x0;
  Int_t ntracks=event->GetNumberOfTracks();

  Double_t probTPC[AliPID::kSPECIES]={0.};
  Double_t probTOF[AliPID::kSPECIES]={0.};
  Double_t probTPCTOF[AliPID::kSPECIES]={0.};
  
  //loop over all tracks
  for (Int_t itrack=0; itrack<ntracks; ++itrack){

    track=(AliVTrack*)event->GetTrack(itrack);

    if ( fTrackFilter->IsSelected(track) ) {

      Double_t mom=track->GetTPCmomentum();
      Int_t ibin=GetMomBin(mom);

      fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
      UInt_t detUsed = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);
     
      if (detUsed != 0) {  // TPC is available
	for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec) {
	  Double_t nSigmaTPC = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)ispec);
	  fProbTPC[ispec]->Fill(mom,probTPC[ispec]);
	  fProbTPCnSigma[ispec]->Fill(nSigmaTPC,probTPC[ispec]);
	  fProbTPCnSigTPCMom[ibin][ispec]->Fill(nSigmaTPC,probTPC[ispec]);
	}

	// compute priors for TPC+TOF, even if we ask just TOF for PID
	fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
	detUsed = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);
	Double_t priors[5]; 	// check priors used for TOF
	fPIDCombined->GetPriors(track,priors,fPIDResponse,detUsed);
	for(Int_t ispec=0;ispec<5;ispec++) fPriorsUsed[ispec]->Fill(TMath::Abs(track->Pt()),priors[ispec]);

	if (detUsed != 0) {  // TOF is available
	  for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec) {
	    Double_t nSigmaTOF = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)ispec);
	    fProbTOF[ispec]->Fill(mom,probTOF[ispec]);
	    fProbTOFnSigma[ispec]->Fill(nSigmaTOF,probTOF[ispec]);
	    fProbTOFnSigTOFMom[ibin][ispec]->Fill(nSigmaTOF,probTOF[ispec]);
	  }
	}

	fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
	detUsed = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);
	if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask() ) {
	  for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec) {
	    Double_t nSigmaTPC = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)ispec);
	    fProbTPCTOF[ispec]->Fill(mom,probTPCTOF[ispec]);
	    fProbTPCTOFnSigmaTPC[ispec]->Fill(nSigmaTPC,probTPCTOF[ispec]);
	    fProbTPCTOFnSigTPCMom[ibin][ispec]->Fill(nSigmaTPC,probTPCTOF[ispec]);
	  }
	}

      }

      fPIDResponse->GetTPCsignalTunedOnData(track);

      fDeDx->Fill(mom,track->GetTPCsignal());
      fDeDxTuned->Fill(mom,track->GetTPCsignalTunedOnData());

    }
  }

  PostData(1, &fHistList);
}

//_________________________________________________________________________________
void AliAnalysisTaskPIDCombined::FillHistogram(const char* name, Double_t x, Double_t weight)
{
  //
  // Fill 1D histogram by name
  //
  ((TH1*)fHistList.FindObject(name))->Fill(x,weight);
}

//_________________________________________________________________________________
void AliAnalysisTaskPIDCombined::FillHistogram(const char* name, Double_t x, Double_t y, Double_t weight)
{
  //
  // Fill 2D histogram by name
  //
  ((TH2*)fHistList.FindObject(name))->Fill(x,y,weight);
}


//_________________________________________________________________________________
Int_t AliAnalysisTaskPIDCombined::GetMomBin(Float_t mom)
{
  //
  // Given momentum return histogram to be filled
  //
  if (mom>0. && mom < 0.5) return 0;
  if (mom>=0.5 && mom < 0.7) return 1;
  if (mom>=0.7 && mom < 1.0) return 2;
  if (mom>=1.0 && mom < 1.5) return 3;
  if (mom>=1.5 && mom < 2.0) return 4;
  return kPtBins-1;
}
Commit	Line	Data
a65a7e70	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	// //
	18	// Basic Analysis Task //
	19	// for PID Analysis //
	20	// //
	21	///////////////////////////////////////////////////////////////////////////
	22
	23	#include <TChain.h>
	24	#include <TH1D.h>
	25	#include <TH2D.h>
	26
	27	#include <AliVEvent.h>
	28	#include <AliInputEventHandler.h>
	29	#include <AliAnalysisManager.h>
	30
	31	#include <AliLog.h>
	32	#include <AliPID.h>
	33	#include <AliPIDCombined.h>
	34	#include <AliPIDResponse.h>
	35
	36	#include "AliAnalysisTaskPIDCombined.h"
	37
	38	const char *AliAnalysisTaskPIDCombined::fgkBinMomDesc[AliAnalysisTaskPIDCombined::kPtBins] = {
	39	" 0 <= p < 0.5 GeV/c",
	40	" 0.5 <= p < 0.7 GeV/c",
	41	" 0.7 <= p < 1.0 GeV/c",
	42	" 1.0 <= p < 1.5 GeV/c",
	43	" 1.5 <= p < 2.0 GeV/c",
	44	" p >= 2.0 GeV/c"
	45	};
	46
	47	ClassImp(AliAnalysisTaskPIDCombined)
	48
	49	//_________________________________________________________________________________
	50	AliAnalysisTaskPIDCombined::AliAnalysisTaskPIDCombined() :
	51	AliAnalysisTaskSE(),
	52	fHistList(),
	53	fProbTPCnSigma(),
	54	fProbTOFnSigma(),
	55	fProbTPCTOFnSigmaTPC(),
	56	fProbTPC(),
	57	fProbTOF(),
	58	fProbTPCTOF(),
	59	fPriors(),
	60	fProbTPCTOFnSigTPCMom(),
	61	fProbTPCnSigTPCMom(),
	62	fProbTOFnSigTOFMom(),
	63	fPriorsUsed(),
	64	fPIDResponse(0x0),
65	fPIDCombined(0x0),
66	fTrackCuts(0x0),
67	fTrackFilter(0x0),
539a5a59	68	fDeDx(NULL),
539a5a59	69	fDeDxTuned(NULL)
a65a7e70	70	{
	71	//
	72	// Constructor
	73	//
	74	}
	75
	76	//_________________________________________________________________________________
	77	AliAnalysisTaskPIDCombined::AliAnalysisTaskPIDCombined(const char *name) :
	78	AliAnalysisTaskSE(name),
	79	fHistList(),
	80	fProbTPCnSigma(),
	81	fProbTOFnSigma(),
	82	fProbTPCTOFnSigmaTPC(),
	83	fProbTPC(),
	84	fProbTOF(),
	85	fProbTPCTOF(),
	86	fPriors(),
	87	fProbTPCTOFnSigTPCMom(),
	88	fProbTPCnSigTPCMom(),
	89	fProbTOFnSigTOFMom(),
	90	fPriorsUsed(),
	91	fPIDResponse(0x0),
	92	fPIDCombined(0x0),
	93	fTrackCuts(0x0),
	94	fTrackFilter(0x0),
539a5a59	95	fDeDx(NULL),
539a5a59	96	fDeDxTuned(NULL)
a65a7e70	97	{
	98	//
	99	// Constructor
	100	//
	101	DefineInput(0,TChain::Class());
	102	DefineOutput(1, TList::Class());
	103	}
	104
	105	//_________________________________________________________________________________
	106	void AliAnalysisTaskPIDCombined::UserCreateOutputObjects()
	107	{
	108	//
	109	// Initialise the framework objects
	110	//
	111
	112
	113	// ------- track cuts
	114	fTrackCuts = new AliESDtrackCuts("fTrackCuts", "Standard");
	115	fTrackCuts->SetAcceptKinkDaughters(kFALSE);
	116	fTrackCuts->SetMinNClustersTPC(80);
	117	fTrackCuts->SetMaxChi2PerClusterTPC(4);
	118	fTrackCuts->SetMaxDCAToVertexXY(3);
	119	fTrackCuts->SetMaxDCAToVertexZ(3);
	120	fTrackCuts->SetRequireTPCRefit(kTRUE);
	121	fTrackCuts->SetRequireITSRefit(kTRUE);
	122	fTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
	123	fTrackFilter = new AliAnalysisFilter("trackFilter");
	124	fTrackFilter->AddCuts(fTrackCuts);
	125
	126
	127
	128	// ------- setup PIDCombined
	129	fPIDCombined=new AliPIDCombined;
	130	fPIDCombined->SetDefaultTPCPriors();
	131	fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
	132
	133	// no light nuclei - no need to call it, this is default
	134	// fPIDCombined->SetSelectedSpecies(AliPID::kSPECIES);
	135
	136
	137	fHistList.Add(new TH1D("nEvents","Number of Evnets;Selection",2,0,2));
	138
	139	for (Int_t ispec=0; ispec<5; ++ispec){
	140
	141
	142	fProbTPC[ispec]=new TH2D(Form("prob%s_mom_TPC",AliPID::ParticleName(ispec)),
	143	Form("%s probability vs. momentum;momentum;probability",AliPID::ParticleName(ispec)),
	144	100,0.,20.,50,0.,1.);
	145	fHistList.Add(fProbTPC[ispec]);
	146	fProbTPCnSigma[ispec]=new TH2D(Form("prob%s_nSigma_TPC",AliPID::ParticleName(ispec)),
	147	Form("%s probability vs. n#sigma;n#sigma;probability",AliPID::ParticleName(ispec)),
	148	20,-5.,5.,50,0.,1.);
	149	fHistList.Add(fProbTPCnSigma[ispec]);
	150
	151	for (Int_t ibin=0;ibin<kPtBins;ibin++) {
	152	fProbTPCnSigTPCMom[ibin][ispec]=new TH2D(Form("prob%s_nSigma_TPC (%s)",AliPID::ParticleName(ispec),fgkBinMomDesc[ibin]),
	153	Form("%s probability vs. n#sigma;n#sigma;probability",AliPID::ParticleName(ispec)),
	154	20,-5.,5.,50,0.,1.);
	155	fHistList.Add(fProbTPCnSigTPCMom[ibin][ispec]);
	156	}
	157
	158
	159
	160	fProbTOF[ispec]=new TH2D(Form("prob%s_mom_TOF",AliPID::ParticleName(ispec)),
161	Form("%s probability vs. momentum;momentum;probability",AliPID::ParticleName(ispec)),
162	100,0.,20.,50,0.,1.);
163	fHistList.Add(fProbTOF[ispec]);
164	fProbTOFnSigma[ispec]=new TH2D(Form("prob%s_nSigma_TOF",AliPID::ParticleName(ispec)),
165	Form("%s probability vs. n#sigma;n#sigma;probability",AliPID::ParticleName(ispec)),
166	20,-5.,5.,50,0.,1.);
167	fHistList.Add(fProbTOFnSigma[ispec]);
168	for (Int_t ibin=0;ibin<kPtBins;ibin++) {
169	fProbTOFnSigTOFMom[ibin][ispec]=new TH2D(Form("prob%s_nSigma_TOF (%s)",AliPID::ParticleName(ispec),fgkBinMomDesc[ibin]),
170	Form("%s probability vs. n#sigma;n#sigma;probability",AliPID::ParticleName(ispec)),
171	20,-5.,5.,50,0.,1.);
172	fHistList.Add(fProbTOFnSigTOFMom[ibin][ispec]);
173	}
174
175
176
177	fProbTPCTOF[ispec]=new TH2D(Form("prob%s_mom_TPCTOF",AliPID::ParticleName(ispec)),
178	Form("%s probability vs. momentum;momentum;probability",AliPID::ParticleName(ispec)),
179	100,0.,20.,50,0.,1.);
180	fHistList.Add(fProbTPCTOF[ispec]);
181	fProbTPCTOFnSigmaTPC[ispec]=new TH2D(Form("probTPCTOF%s_nSigma_TPC",AliPID::ParticleName(ispec)),
182	Form("%s TPCTOF probability vs. n#sigmaTPC;n#sigma;probability",AliPID::ParticleName(ispec)),
183	20,-5.,5.,50,0.,1.);
184	fHistList.Add(fProbTPCTOFnSigmaTPC[ispec]);
185	for (Int_t ibin=0;ibin<kPtBins;ibin++) {
186	fProbTPCTOFnSigTPCMom[ibin][ispec]=new TH2D(Form("probTPCTOF%s_nSigma_TPC (%s)",AliPID::ParticleName(ispec),fgkBinMomDesc[ibin]),
187	Form("%s probability vs. n#sigma;n#sigma;probability",AliPID::ParticleName(ispec)),
188	20,-5.,5.,50,0.,1.);
189	fHistList.Add(fProbTPCTOFnSigTPCMom[ibin][ispec]);
190	}
191
192
193
194	// basic priors
195	fPriors[ispec]=new TH1F(Form("%s_priors",AliPID::ParticleName(ispec)),
196	Form("%s priors vs momentum",AliPID::ParticleName(ispec)),
197	100,0.,20.);
198	fHistList.Add(fPriors[ispec]);
199	switch (ispec) {
200	case AliPID::kElectron:
201	for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.02);
202	break;
203	case AliPID::kMuon:
204	for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.02);
205	break;
206	case AliPID::kPion:
207	for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.56);
208	break;
209	case AliPID::kKaon:
210	for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.20);
211	break;
212	case AliPID::kProton:
213	for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.20);
214	break;
215	default:
216	break;
217	}
218	fPIDCombined->SetPriorDistribution((AliPID::EParticleType)ispec,fPriors[ispec]);
219
220	// priors used
221	fPriorsUsed[ispec] = new TH2D(Form("%s_priorsUsed",AliPID::ParticleName(ispec)),
222	Form("%s priors vs transverse momentum;p_{t} (GeV/c);priors",AliPID::ParticleName(ispec)),
223	100,0.,20.,101,0,1.01);
224	fHistList.Add(fPriorsUsed[ispec]);
225	}
226
227
539a5a59	228	fDeDx = new TH2D("hDeDx",";p_{TPC};dE/dx (a.u.)",500,0,5,500,0,500);
	229	fHistList.Add(fDeDx);
	230	fDeDxTuned = new TH2D("hDeDxTuned",";p_{TPC};dE/dx (a.u.)",500,0,5,500,0,500);
	231	fHistList.Add(fDeDxTuned);
	232
a65a7e70	233	fHistList.SetOwner();
	234	PostData(1,&fHistList);
	235
	236
	237	}
	238
	239	//_________________________________________________________________________________
	240	void AliAnalysisTaskPIDCombined::UserExec(Option_t *)
	241	{
	242	//
	243	// Main loop. Called for every event
	244	//
	245	AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
	246	AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
	247	fPIDResponse=inputHandler->GetPIDResponse();
	248	if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");
	249
	250	// Printf(" ---------------------- UserExec PID task ---------------------");
	251
	252	FillHistogram("nEvents",0.);
	253
	254	AliVEvent *event=InputEvent();
	255	AliVTrack *track=0x0;
	256	Int_t ntracks=event->GetNumberOfTracks();
	257
	258	Double_t probTPC[AliPID::kSPECIES]={0.};
	259	Double_t probTOF[AliPID::kSPECIES]={0.};
	260	Double_t probTPCTOF[AliPID::kSPECIES]={0.};
	261
	262	//loop over all tracks
	263	for (Int_t itrack=0; itrack<ntracks; ++itrack){
	264
	265	track=(AliVTrack*)event->GetTrack(itrack);
	266
	267	if ( fTrackFilter->IsSelected(track) ) {
	268
	269	Double_t mom=track->GetTPCmomentum();
	270	Int_t ibin=GetMomBin(mom);
	271
	272	fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
	273	UInt_t detUsed = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);
	274
	275	if (detUsed != 0) { // TPC is available
	276	for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec) {
	277	Double_t nSigmaTPC = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)ispec);
	278	fProbTPC[ispec]->Fill(mom,probTPC[ispec]);
	279	fProbTPCnSigma[ispec]->Fill(nSigmaTPC,probTPC[ispec]);
	280	fProbTPCnSigTPCMom[ibin][ispec]->Fill(nSigmaTPC,probTPC[ispec]);
	281	}
	282
	283	// compute priors for TPC+TOF, even if we ask just TOF for PID
	284	fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
	285	detUsed = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);
	286	Double_t priors[5]; // check priors used for TOF
	287	fPIDCombined->GetPriors(track,priors,fPIDResponse,detUsed);
	288	for(Int_t ispec=0;ispec<5;ispec++) fPriorsUsed[ispec]->Fill(TMath::Abs(track->Pt()),priors[ispec]);
	289
	290	if (detUsed != 0) { // TOF is available
	291	for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec) {
	292	Double_t nSigmaTOF = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)ispec);
	293	fProbTOF[ispec]->Fill(mom,probTOF[ispec]);
	294	fProbTOFnSigma[ispec]->Fill(nSigmaTOF,probTOF[ispec]);
	295	fProbTOFnSigTOFMom[ibin][ispec]->Fill(nSigmaTOF,probTOF[ispec]);
	296	}
297	}
298
299	fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF\|AliPIDResponse::kDetTPC);
300	detUsed = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);
301	if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask() ) {
302	for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec) {
303	Double_t nSigmaTPC = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)ispec);
304	fProbTPCTOF[ispec]->Fill(mom,probTPCTOF[ispec]);
305	fProbTPCTOFnSigmaTPC[ispec]->Fill(nSigmaTPC,probTPCTOF[ispec]);
306	fProbTPCTOFnSigTPCMom[ibin][ispec]->Fill(nSigmaTPC,probTPCTOF[ispec]);
307	}
308	}
309
310	}
311
539a5a59	312	fPIDResponse->GetTPCsignalTunedOnData(track);
	313
	314	fDeDx->Fill(mom,track->GetTPCsignal());
	315	fDeDxTuned->Fill(mom,track->GetTPCsignalTunedOnData());
a65a7e70	316
	317	}
	318	}
	319
	320	PostData(1, &fHistList);
	321	}
	322
	323	//_________________________________________________________________________________
	324	void AliAnalysisTaskPIDCombined::FillHistogram(const char* name, Double_t x, Double_t weight)
	325	{
	326	//
	327	// Fill 1D histogram by name
	328	//
	329	((TH1*)fHistList.FindObject(name))->Fill(x,weight);
	330	}
	331
	332	//_________________________________________________________________________________
	333	void AliAnalysisTaskPIDCombined::FillHistogram(const char* name, Double_t x, Double_t y, Double_t weight)
	334	{
	335	//
	336	// Fill 2D histogram by name
	337	//
	338	((TH2*)fHistList.FindObject(name))->Fill(x,y,weight);
	339	}
	340
	341
	342	//_________________________________________________________________________________
	343	Int_t AliAnalysisTaskPIDCombined::GetMomBin(Float_t mom)
	344	{
	345	//
	346	// Given momentum return histogram to be filled
	347	//
	348	if (mom>0. && mom < 0.5) return 0;
	349	if (mom>=0.5 && mom < 0.7) return 1;
	350	if (mom>=0.7 && mom < 1.0) return 2;
	351	if (mom>=1.0 && mom < 1.5) return 3;
	352	if (mom>=1.5 && mom < 2.0) return 4;
	353	return kPtBins-1;
	354	}
	355