fPIDResponse(0x0),\r
fPIDCombined(0x0),\r
fTrackCuts(0x0),\r
- fTrackFilter(0x0)\r
+ fTrackFilter(0x0),\r
+ fDeDx(NULL),
+ fDeDxTuned(NULL)
{\r
//\r
// Constructor\r
fPIDResponse(0x0),\r
fPIDCombined(0x0),\r
fTrackCuts(0x0),\r
- fTrackFilter(0x0)\r
+ fTrackFilter(0x0),\r
+ fDeDx(NULL),
+ fDeDxTuned(NULL)
{\r
//\r
// Constructor\r
}\r
\r
\r
-\r
+ 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();\r
PostData(1,&fHistList);\r
\r
\r
}\r
\r
+ fPIDResponse->GetTPCsignalTunedOnData(track);
+
+ fDeDx->Fill(mom,track->GetTPCsignal());
+ fDeDxTuned->Fill(mom,track->GetTPCsignalTunedOnData());
\r
}\r
}\r
AliPIDCombined *fPIDCombined; //! combined PID object\r
AliESDtrackCuts *fTrackCuts; //! track selection\r
AliAnalysisFilter *fTrackFilter; //! track filter\r
+
+ TH2D *fDeDx; //! histo with the dedx
+ TH2D *fDeDxTuned; //! histo to check the dedx tuning in MC
+
\r
AliAnalysisTaskPIDCombined(const AliAnalysisTaskPIDCombined &c);\r
AliAnalysisTaskPIDCombined& operator= (const AliAnalysisTaskPIDCombined &c);\r
fPIDResponse(0x0),
fRun(0),
fOldRun(0),
-fRecoPass(0)
+fRecoPass(0),
+fIsTunedOnData(kFALSE),
+fRecoPassTuned(0)
{
//
// Dummy constructor
fPIDResponse(0x0),
fRun(0),
fOldRun(0),
-fRecoPass(0)
+fRecoPass(0),
+fIsTunedOnData(kFALSE),
+fRecoPassTuned(0)
{
//
// Default constructor
fPIDResponse->SetOADBPath(AliAnalysisManager::GetOADBPath());
if (!fOADBPath.IsNull()) fPIDResponse->SetOADBPath(fOADBPath.Data());
+
+ if(fIsTunedOnData) fPIDResponse->SetTunedOnData(kTRUE,fRecoPassTuned);
}
//______________________________________________________________________________
void SetOADBPath(const char* path) {fOADBPath=path;}
const char* GetOADBPath() const { return fOADBPath.Data(); }
+ void SetTuneOnData(Bool_t flag,Int_t recopass){fIsTunedOnData=flag;fRecoPassTuned=recopass;};
private:
Bool_t fIsMC; // If we run on MC data
Int_t fRun; //! current run number
Int_t fOldRun; //! current run number
Int_t fRecoPass; //! reconstruction pass
+
+ Bool_t fIsTunedOnData; // flag to tune MC on data (dE/dx)
+ Int_t fRecoPassTuned; // Reco pass tuned on data for MC
//
void SetRecoInfo();
fProdVertex(NULL),
fTrackPhiOnEMCal(-999),
fTrackEtaOnEMCal(-999),
+ fTPCsignalTuned(0),
fAODEvent(NULL)
{
// default constructor
fProdVertex(prodVertex),
fTrackPhiOnEMCal(-999),
fTrackEtaOnEMCal(-999),
+ fTPCsignalTuned(0),
fAODEvent(NULL)
{
// constructor
fProdVertex(prodVertex),
fTrackPhiOnEMCal(-999),
fTrackEtaOnEMCal(-999),
+ fTPCsignalTuned(0),
fAODEvent(NULL)
{
// constructor
fProdVertex(trk.fProdVertex),
fTrackPhiOnEMCal(trk.fTrackPhiOnEMCal),
fTrackEtaOnEMCal(trk.fTrackEtaOnEMCal),
+ fTPCsignalTuned(trk.fTPCsignalTuned),
fAODEvent(trk.fAODEvent)
{
// Copy constructor
fCaloIndex = trk.fCaloIndex;
fTrackPhiOnEMCal = trk.fTrackPhiOnEMCal;
fTrackEtaOnEMCal = trk.fTrackEtaOnEMCal;
+ fTPCsignalTuned = trk.fTPCsignalTuned;
delete fCovMatrix;
if(trk.fCovMatrix) fCovMatrix=new AliAODRedCov<6>(*trk.fCovMatrix);
//pid signal interface
Double_t GetITSsignal() const { return fDetPid?fDetPid->GetITSsignal():0.; }
Double_t GetTPCsignal() const { return fDetPid?fDetPid->GetTPCsignal():0.; }
+ Double_t GetTPCsignalTunedOnData() const { return fTPCsignalTuned;}
+ void SetTPCsignalTunedOnData(Double_t signal) {fTPCsignalTuned = signal;}
UShort_t GetTPCsignalN() const { return fDetPid?fDetPid->GetTPCsignalN():0; }
virtual AliTPCdEdxInfo* GetTPCdEdxInfo() const {return fDetPid?fDetPid->GetTPCdEdxInfo():0;}
Double_t GetTPCmomentum() const { return fDetPid?fDetPid->GetTPCmomentum():0.; }
UChar_t GetTRDntrackletsPID() const;
void GetHMPIDpid(Double_t */*p*/) const { return; } // TODO: To be implemented properly with the new HMPID object
- const AliAODEvent* GetAODEvent(){return fAODEvent;}
+ const AliAODEvent* GetAODEvent() const {return fAODEvent;}
void SetAODEvent(const AliAODEvent* ptr){fAODEvent = ptr;}
AliAODPid *GetDetPid() const { return fDetPid; }
Double_t fTrackPhiOnEMCal; // phi of track after being propagated to 430cm
Double_t fTrackEtaOnEMCal; // eta of track after being propagated to 430cm
+ Double_t fTPCsignalTuned; //! TPC signal tuned on data when using MC
+
const AliAODEvent* fAODEvent; //!
- ClassDef(AliAODTrack, 17);
+ ClassDef(AliAODTrack, 18);
};
inline Bool_t AliAODTrack::IsPrimaryCandidate() const
// Origin: Rosa Romita, GSI, r.romita@gsi.de
//-----------------------------------------------------------------
+#include "TRandom.h"
#include "AliLog.h"
#include "AliPID.h"
#include "AliAODpidUtil.h"
#include "AliAODPid.h"
#include "AliTRDPIDResponse.h"
#include "AliESDtrack.h"
+#include "AliAODMCHeader.h"
+#include "AliAODMCParticle.h"
ClassImp(AliAODpidUtil)
return 0;
}
//_________________________________________________________________________
+Float_t AliAODpidUtil::GetTPCsignalTunedOnData(const AliVTrack *t) const {
+ AliAODTrack *track = (AliAODTrack *) t;
+ Float_t dedx = track->GetTPCsignalTunedOnData();
+ if(dedx > 0) return dedx;
+
+ Double_t mom = t->GetTPCmomentum();
+
+ dedx = t->GetTPCsignal();
+ track->SetTPCsignalTunedOnData(dedx);
+
+ if(dedx < 20) return dedx;
+
+
+ AliPID::EParticleType type = AliPID::kPion;
+
+ AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(track->GetAODEvent()->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
+ if (mcHeader) {
+ TClonesArray *mcArray = (TClonesArray*)track->GetAODEvent()->GetList()->FindObject(AliAODMCParticle::StdBranchName());
+
+ Bool_t kGood = kTRUE;
+
+ Int_t iS = TMath::Abs(((AliAODMCParticle*)mcArray->At(TMath::Abs(t->GetLabel())))->GetPdgCode());
+ if(iS==AliPID::ParticleCode(AliPID::kElectron)){
+ type = AliPID::kElectron;
+ }
+ else if(iS==AliPID::ParticleCode(AliPID::kMuon)){
+ type = AliPID::kMuon;
+ }
+ else if(iS==AliPID::ParticleCode(AliPID::kPion)){
+ type = AliPID::kPion;
+ }
+ else if(iS==AliPID::ParticleCode(AliPID::kKaon)){
+ type = AliPID::kKaon;
+ }
+ else if(iS==AliPID::ParticleCode(AliPID::kProton)){
+ type = AliPID::kProton;
+ }
+ else if(iS==AliPID::ParticleCode(AliPID::kDeuteron)){ // d
+ type = AliPID::kDeuteron;
+ }
+ else if(iS==AliPID::ParticleCode(AliPID::kTriton)){ // t
+ type = AliPID::kTriton;
+ }
+ else if(iS==AliPID::ParticleCode(AliPID::kHe3)){ // 3He
+ type = AliPID::kHe3;
+ }
+ else if(iS==AliPID::ParticleCode(AliPID::kAlpha)){ // 4He
+ type = AliPID::kAlpha;
+ }
+ else
+ kGood = kFALSE;
+
+ if(kGood){
+ Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
+ Double_t sigma=fTPCResponse.GetExpectedSigma(mom,t->GetTPCsignalN(),type);
+ dedx = gRandom->Gaus(bethe,sigma);
+
+ if(iS == AliPID::ParticleCode(AliPID::kHe3) || iS == AliPID::ParticleCode(AliPID::kAlpha)) dedx *= 5;
+ }
+
+ }
+
+ track->SetTPCsignalTunedOnData(dedx);
+ return dedx;
+}
+//_________________________________________________________________________
void AliAODpidUtil::MakeTPCPID(const AliAODTrack *track,Double_t *p) const
{
//
virtual Float_t NumberOfSigmasTOF(const AliVParticle *vtrack, AliPID::EParticleType type) const;
+ Float_t GetTPCsignalTunedOnData(const AliVTrack *t) const;
+
private:
ClassDef(AliAODpidUtil,3) // PID calculation class
#include "TArrayI.h"
#include "TArrayF.h"
+#include "TRandom.h"
#include "AliLog.h"
#include "AliPID.h"
#include "AliTOFHeader.h"
#include "AliESDpid.h"
#include "AliESDEvent.h"
#include "AliESDtrack.h"
+#include "AliMCEvent.h"
+#include "AliMCEventHandler.h"
+#include "AliAnalysisManager.h"
+
ClassImp(AliESDpid)
return 0;
}
//_________________________________________________________________________
+Float_t AliESDpid::GetTPCsignalTunedOnData(const AliVTrack *t) const {
+ AliESDtrack *track = (AliESDtrack *) t;
+ Float_t dedx = track->GetTPCsignalTunedOnData();
+
+ if(dedx > 0) return dedx;
+
+ Double_t mom = t->GetTPCmomentum();
+ dedx = t->GetTPCsignal();
+ track->SetTPCsignalTunedOnData(dedx);
+ if(dedx < 20) return dedx;
+
+ AliPID::EParticleType type = AliPID::kPion;
+
+ AliMCEventHandler* eventHandler=NULL;
+ eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
+ if (eventHandler) {
+ AliMCEvent* mcEvent = eventHandler->MCEvent();
+ if(mcEvent){
+ Bool_t kGood = kTRUE;
+ AliMCParticle *MCpart = (AliMCParticle *) mcEvent->GetTrack(TMath::Abs(t->GetLabel()));
+ TParticle *part = MCpart->Particle();
+
+ Int_t iS = TMath::Abs(part->GetPdgCode());
+
+ if(iS==AliPID::ParticleCode(AliPID::kElectron)){
+ type = AliPID::kElectron;
+ }
+ else if(iS==AliPID::ParticleCode(AliPID::kMuon)){
+ type = AliPID::kMuon;
+ }
+ else if(iS==AliPID::ParticleCode(AliPID::kPion)){
+ type = AliPID::kPion;
+ }
+ else if(iS==AliPID::ParticleCode(AliPID::kKaon)){
+ type = AliPID::kKaon;
+ }
+ else if(iS==AliPID::ParticleCode(AliPID::kProton)){
+ type = AliPID::kProton;
+ }
+ else if(iS==AliPID::ParticleCode(AliPID::kDeuteron)){ // d
+ type = AliPID::kDeuteron;
+ }
+ else if(iS==AliPID::ParticleCode(AliPID::kTriton)){ // t
+ type = AliPID::kTriton;
+ }
+ else if(iS==AliPID::ParticleCode(AliPID::kHe3)){ // 3He
+ type = AliPID::kHe3;
+ }
+ else if(iS==AliPID::ParticleCode(AliPID::kAlpha)){ // 4He
+ type = AliPID::kAlpha;
+ }
+ else
+ kGood = kFALSE;
+
+ if(kGood){
+ Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
+ Double_t sigma=fTPCResponse.GetExpectedSigma(mom,t->GetTPCsignalN(),type);
+ dedx = gRandom->Gaus(bethe,sigma);
+ if(iS == AliPID::ParticleCode(AliPID::kHe3) || iS == AliPID::ParticleCode(AliPID::kAlpha)) dedx *= 5;
+ }
+ }
+ }
+
+ track->SetTPCsignalTunedOnData(dedx);
+ return dedx;
+}
+//_________________________________________________________________________
void AliESDpid::MakeTPCPID(AliESDtrack *track) const
{
//
void SetNMaxSigmaTOFTPCMismatch(Float_t range) {fRangeTOFMismatch=range;}
Float_t GetNMaxSigmaTOFTPCMismatch() const {return fRangeTOFMismatch;}
+ Float_t GetTPCsignalTunedOnData(const AliVTrack *t) const;
+
private:
Float_t fRangeTOFMismatch; // nSigma max for TOF matching with TPC
fGlobalChi2(0),
fITSsignal(0),
fTPCsignal(0),
+ fTPCsignalTuned(0),
fTPCsignalS(0),
fTPCdEdxInfo(0),
fTRDsignal(0),
fGlobalChi2(track.fGlobalChi2),
fITSsignal(track.fITSsignal),
fTPCsignal(track.fTPCsignal),
+ fTPCsignalTuned(track.fTPCsignalTuned),
fTPCsignalS(track.fTPCsignalS),
fTPCdEdxInfo(0),
fTRDsignal(track.fTRDsignal),
fGlobalChi2(0),
fITSsignal(0),
fTPCsignal(0),
+ fTPCsignalTuned(0),
fTPCsignalS(0),
fTPCdEdxInfo(0),
fTRDsignal(0),
fGlobalChi2(0),
fITSsignal(0),
fTPCsignal(0),
+ fTPCsignalTuned(0),
fTPCsignalS(0),
fTPCdEdxInfo(0),
fTRDsignal(0),
fITSsignal = source.fITSsignal;
for (Int_t i=0;i<4;i++) {fITSdEdxSamples[i]=source.fITSdEdxSamples[i];}
fTPCsignal = source.fTPCsignal;
+ fTPCsignalTuned = source.fTPCsignalTuned;
fTPCsignalS = source.fTPCsignalS;
for(int i = 0; i< 4;++i){
fTPCPoints[i] = source.fTPCPoints[i];
track.fTPCchi2 = fTPCchi2;
track.fTPCchi2Iter1 = fTPCchi2Iter1;
track.fTPCsignal = fTPCsignal;
+ track.fTPCsignalTuned = fTPCsignalTuned;
track.fTPCsignalS = fTPCsignalS;
for(int i = 0;i<4;++i)track.fTPCPoints[i] = fTPCPoints[i];
fTPCClusterMap = 0;
fTPCSharedMap = 0;
fTPCsignal= 0;
+ fTPCsignalTuned= 0;
fTPCsignalS= 0;
fTPCsignalN= 0;
for (Int_t i=0;i<AliPID::kSPECIES;i++) fTPCr[i]=0;
void SetTPCsignal(Float_t signal, Float_t sigma, UChar_t npoints){
fTPCsignal = signal; fTPCsignalS = sigma; fTPCsignalN = npoints;
}
+ void SetTPCsignalTunedOnData(Float_t signal){
+ fTPCsignalTuned = signal;
+ }
void SetTPCdEdxInfo(AliTPCdEdxInfo * dEdxInfo);
AliTPCdEdxInfo * GetTPCdEdxInfo() const {return fTPCdEdxInfo;}
Double_t GetTPCsignal() const {return fTPCsignal;}
+ Double_t GetTPCsignalTunedOnData() const {return fTPCsignalTuned;}
Double_t GetTPCsignalSigma() const {return fTPCsignalS;}
UShort_t GetTPCsignalN() const {return fTPCsignalN;}
Double_t GetTPCmomentum() const {return fIp?fIp->GetP():GetP();}
Double32_t fITSdEdxSamples[4]; // [0.,0.,10] ITS dE/dx samples
Double32_t fTPCsignal; // [0.,0.,10] detector's PID signal
+ Double32_t fTPCsignalTuned; //! [0.,0.,10] detector's PID signal tuned on data when using MC
Double32_t fTPCsignalS; // [0.,0.,10] RMS of dEdx measurement
AliTPCdEdxInfo * fTPCdEdxInfo; // object containing dE/dx information for different pad regions
Double32_t fTPCPoints[4]; // [0.,0.,10] TPC points -first, max. dens, last and max density
static bool fgkOnlineMode; //! indicate the online mode to skip some of the functionality
AliESDtrack & operator=(const AliESDtrack & );
- ClassDef(AliESDtrack,65) //ESDtrack
+ ClassDef(AliESDtrack,66) //ESDtrack
};
fTOFPIDParams(0x0),
fEMCALPIDParams(0x0),
fCurrentEvent(0x0),
-fCurrCentrality(0.0)
+fCurrCentrality(0.0),
+fTuneMConData(kFALSE)
{
//
// default ctor
fTOFPIDParams(0x0),
fEMCALPIDParams(0x0),
fCurrentEvent(0x0),
-fCurrCentrality(0.0)
+fCurrCentrality(0.0),
+fTuneMConData(kFALSE)
{
//
// copy ctor
Double_t dedx=track->GetTPCsignal();
Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
+ if(fTuneMConData) dedx = this->GetTPCsignalTunedOnData(track);
+
for (Int_t j=0; j<AliPID::kSPECIES; j++) {
AliPID::EParticleType type=AliPID::EParticleType(j);
Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
TString datatype="DATA";
//in case of mc fRecoPass is per default 1
if (fIsMC) {
- datatype="MC";
- fRecoPass=1;
+ datatype="MC";
+ if(!fTuneMConData) datatype="MC";
+ fRecoPass=1;
}
//
// period
TString period=fLHCperiod;
- if (fIsMC) period=fMCperiodTPC;
+ if (fIsMC && !fTuneMConData) period=fMCperiodTPC;
AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
Bool_t found=kFALSE;
//set the new PID splines
//
if (fArrPidResponseMaster){
+ Int_t recopass = fRecoPass;
+ if(fTuneMConData) recopass = fRecoPassUser;
TObject *grAll=0x0;
//for MC don't use period information
// if (fIsMC) period="[A-Z0-9]*";
//for MC use MC period information
//pattern for the default entry (valid for all particles)
- TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
-
+ TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
+
//loop over entries and filter them
for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp){
TObject *responseFunction=fArrPidResponseMaster->At(iresp);
void SetTOFtail(Float_t tail=1.1){if(tail > 0) fTOFtail=tail; else printf("TOF tail should be greater than 0 (nothing done)\n");};
void SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option);
+ virtual Float_t GetTPCsignalTunedOnData(const AliVTrack *t) const {return t->GetTPCsignal();};
+ Bool_t IsTunedOnData() const {return fTuneMConData;};
+ void SetTunedOnData(Bool_t flag=kTRUE,Int_t recoPass=0){fTuneMConData = flag; if(recoPass>0) fRecoPassUser = recoPass;};
+
AliPIDResponse(const AliPIDResponse &other);
AliPIDResponse& operator=(const AliPIDResponse &other);
Float_t fCurrCentrality; //! current centrality
+ Bool_t fTuneMConData; // switch to force the MC to be similar to data (dE/dx)
+
void ExecNewRun();
//
//
void SetRecoInfo();
- ClassDef(AliPIDResponse, 7); //PID response handling
+ ClassDef(AliPIDResponse, 8); //PID response handling
};
inline Float_t AliPIDResponse::NumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const {
Double_t sig = track->GetTPCsignal();
UInt_t sigN = track->GetTPCsignalN();
+ if(fTuneMConData) sig = GetTPCsignalTunedOnData(track);
+
Double_t nSigma = -999.;
if (sigN>0) nSigma=fTPCResponse.GetNumberOfSigmas(mom,sig,sigN,type);
// assigned clusters and/or the track dip angle, for example.
//
- Double_t mass=AliPID::ParticleMass(n);
+ Double_t mass=AliPID::ParticleMassZ(n);
if (!fUseDatabase||fResponseFunctions.GetEntriesFast()>AliPID::kUnknown) return Bethe(mom/mass);
//
TSpline3 * responseFunction = (TSpline3 *) fResponseFunctions.UncheckedAt(n);
\r
virtual Double_t GetITSsignal() const {return 0.;}\r
virtual Double_t GetTPCsignal() const {return 0.;}\r
+ virtual Double_t GetTPCsignalTunedOnData() const {return 0.;}
virtual UShort_t GetTPCsignalN() const {return 0 ;}\r
virtual Double_t GetTPCmomentum() const {return 0.;}\r
virtual Double_t GetTOFsignal() const {return 0.;}\r