diff --git a/input/MPDTPCPidResponseVhelle.root b/input/MPDTPCPidResponseVhelle.root
new file mode 100644
index 0000000000000000000000000000000000000000..235ccee44f1332175ed7c795d493c12babf6cd68
Binary files /dev/null and b/input/MPDTPCPidResponseVhelle.root differ
diff --git a/mpddata/MpdFillDstTask.cxx b/mpddata/MpdFillDstTask.cxx
index c52653c55bd388c903a2220b7d0cecf98ff9b9bf..fd057ca6f4536f973cb8ad93eccdfbe915a88704 100644
--- a/mpddata/MpdFillDstTask.cxx
+++ b/mpddata/MpdFillDstTask.cxx
@@ -34,231 +34,230 @@ using std::cout;
using std::endl;
// ----- constructor with names ------------------------------------
+
MpdFillDstTask::MpdFillDstTask(const char *name, const char *title)
- : FairTask(name)
-{
+: FairTask(name) {
fEvent = NULL;
fZdcSkeletonesSaved = 0;
}
// ----- Destructor -----------------------------------------------
-MpdFillDstTask::~MpdFillDstTask()
-{
- if (fEvent) delete fEvent;
+
+MpdFillDstTask::~MpdFillDstTask() {
+ if (fEvent) delete fEvent;
}
// -------------------------------------------------------------------
-InitStatus MpdFillDstTask::Init()
-{
- fEvent = new MpdEvent();
-
- FairRootManager *manager = FairRootManager::Instance();
-
- fKFTracks = (TClonesArray *) manager->GetObject("TpcKalmanTrack");
- fKFEctTracks = (TClonesArray *) manager->GetObject("EctTrack");
- fMCTracks = (TClonesArray *) manager->GetObject("MCTrack");
- fMCEventHeader = (FairMCEventHeader*) manager->GetObject("MCEventHeader.");
- fTofMatching = (TClonesArray *) manager->GetObject("TOFMatching");
- fEtofMatching = (TClonesArray *) manager->GetObject("ETOFMatching");
- fVertex = (TClonesArray *) manager->GetObject("Vertex"); //AZ
-
- fELossZdc1Value = (TClonesArray *) manager->GetObject("ELossZdc1Value"); //EL
- fELossZdc2Value = (TClonesArray *) manager->GetObject("ELossZdc2Value"); //EL
- fELossZdc1Histo = (TClonesArray *) manager->GetObject("ELossZdc1Histo"); //EL
- fELossZdc2Histo = (TClonesArray *) manager->GetObject("ELossZdc2Histo"); //EL
- fHistZdc1En = (TH2F *) manager->GetObject("HistZdc1En"); //EL
- fHistZdc2En = (TH2F *) manager->GetObject("HistZdc2En"); //EL
- fZdcSkeletonesSaved=0;
-
- FairRootManager::Instance()->Register("MCEventHeader.","MC", fMCEventHeader, kTRUE);
- FairRootManager::Instance()->Register("MCTrack", "MC", fMCTracks, kTRUE);
- //FairRootManager::Instance()->Register("MpdEvent","MpdEvents", fEvents, kTRUE);
- FairRootManager::Instance()->Register("MPDEvent.","MpdEvent", fEvent, kTRUE);
- // FairRootManager::Instance()->Register("EZdc1","EZdc1", fELossZdc1Histo, kTRUE);
- // FairRootManager::Instance()->Register("EZdc2","EZdc2", fELossZdc2Histo, kTRUE);
-
- //fhTrackMotherId = new TH1F("TrackMotherId","",1000,-5,995);
- //fhTrackPrimaryPDG = new TH1F("TrackPrimaryPDG","",1000,-500,500);
- //fhTrackVertex = new TH2F("TrackVertex","",1000,-500,500,1000,-500,500);
- //fhTruthVertex = new TH2F("TruthVertex","",1000,-500,500,1000,-500,500);
-
- return kSUCCESS;
+
+InitStatus MpdFillDstTask::Init() {
+ fEvent = new MpdEvent();
+
+ FairRootManager *manager = FairRootManager::Instance();
+
+ fKFTracks = (TClonesArray *) manager->GetObject("TpcKalmanTrack");
+ fKFEctTracks = (TClonesArray *) manager->GetObject("EctTrack");
+ fMCTracks = (TClonesArray *) manager->GetObject("MCTrack");
+ fMCEventHeader = (FairMCEventHeader*) manager->GetObject("MCEventHeader.");
+ fTofMatching = (TClonesArray *) manager->GetObject("TOFMatching");
+ fEtofMatching = (TClonesArray *) manager->GetObject("ETOFMatching");
+ fVertex = (TClonesArray *) manager->GetObject("Vertex"); //AZ
+
+ fELossZdc1Value = (TClonesArray *) manager->GetObject("ELossZdc1Value"); //EL
+ fELossZdc2Value = (TClonesArray *) manager->GetObject("ELossZdc2Value"); //EL
+ fELossZdc1Histo = (TClonesArray *) manager->GetObject("ELossZdc1Histo"); //EL
+ fELossZdc2Histo = (TClonesArray *) manager->GetObject("ELossZdc2Histo"); //EL
+ fHistZdc1En = (TH2F *) manager->GetObject("HistZdc1En"); //EL
+ fHistZdc2En = (TH2F *) manager->GetObject("HistZdc2En"); //EL
+ fZdcSkeletonesSaved = 0;
+
+ FairRootManager::Instance()->Register("MCEventHeader.", "MC", fMCEventHeader, kTRUE);
+ FairRootManager::Instance()->Register("MCTrack", "MC", fMCTracks, kTRUE);
+ //FairRootManager::Instance()->Register("MpdEvent","MpdEvents", fEvents, kTRUE);
+ FairRootManager::Instance()->Register("MPDEvent.", "MpdEvent", fEvent, kTRUE);
+ // FairRootManager::Instance()->Register("EZdc1","EZdc1", fELossZdc1Histo, kTRUE);
+ // FairRootManager::Instance()->Register("EZdc2","EZdc2", fELossZdc2Histo, kTRUE);
+
+ //fhTrackMotherId = new TH1F("TrackMotherId","",1000,-5,995);
+ //fhTrackPrimaryPDG = new TH1F("TrackPrimaryPDG","",1000,-500,500);
+ //fhTrackVertex = new TH2F("TrackVertex","",1000,-500,500,1000,-500,500);
+ //fhTruthVertex = new TH2F("TruthVertex","",1000,-500,500,1000,-500,500);
+
+ return kSUCCESS;
}
// -------------------------------------------------------------------
-void MpdFillDstTask::Exec(Option_t * option)
-{
- Reset(); // Clear previos event information
-
- if (!fZdcSkeletonesSaved) { // empty skeletones saved only once
- FairRootManager* ioman = FairRootManager::Instance();
- if (fHistZdc1En) {
- fHistZdc1En->SetDirectory((TFile*)ioman->GetOutFile());
- fHistZdc2En->SetDirectory((TFile*)ioman->GetOutFile());
- fHistZdc1En->Write();
- fHistZdc2En->Write();
- fZdcSkeletonesSaved=1;
+
+void MpdFillDstTask::Exec(Option_t * option) {
+ Reset(); // Clear previos event information
+
+ if (!fZdcSkeletonesSaved) { // empty skeletones saved only once
+ FairRootManager* ioman = FairRootManager::Instance();
+ if (fHistZdc1En) {
+ fHistZdc1En->SetDirectory((TFile*) ioman->GetOutFile());
+ fHistZdc2En->SetDirectory((TFile*) ioman->GetOutFile());
+ fHistZdc1En->Write();
+ fHistZdc2En->Write();
+ fZdcSkeletonesSaved = 1;
+ }
+ }
+
+ Int_t nReco = fKFTracks ? fKFTracks->GetEntriesFast() : 0;
+ Int_t nEctReco = fKFEctTracks ? fKFEctTracks->GetEntriesFast() : 0;
+ cout << "\n-I- [MpdFillDstTask::Exec] " << nReco + nEctReco << " reconstruced tracks to write" << endl;
+
+ FairRunAna *fRun = FairRunAna::Instance();
+ fEvent->SetRunInfoRunId(fRun->GetRunId());
+
+ FairField *fPar = fRun->GetField();
+ fEvent->SetRunInfoMagneticFieldZ(fPar->GetType());
+
+ if (fVertex) {
+ MpdVertex *vertex = (MpdVertex*) fVertex->UncheckedAt(0);
+ fEvent->SetPrimaryVerticesX(vertex->GetX());
+ fEvent->SetPrimaryVerticesY(vertex->GetY());
+ fEvent->SetPrimaryVerticesZ(vertex->GetZ());
+ } else {
+ // Vertex info is not available
+ fEvent->SetPrimaryVerticesX(0.);
+ fEvent->SetPrimaryVerticesY(0.);
+ fEvent->SetPrimaryVerticesZ(0.);
+ }
+
+ // check clone track into ECT
+ typedef std::set<Int_t> trackSet;
+ trackSet EctTrackSet;
+ for (Int_t index = 0; index < nEctReco; index++) // cycle by ECT KF tracks
+ {
+ MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fKFEctTracks->UncheckedAt(index);
+ if (track->IsFromTpc()) EctTrackSet.insert(trackSet::value_type(track->GetTpcIndex()));
+ }
+
+ Int_t nMatching = fTofMatching ? fTofMatching->GetEntriesFast() : 0;
+ MpdTofMatchingData *pMatchingData;
+ bool matchingDataExist;
+
+ //AZ MpdParticleIdentification *identificator = new MpdParticleIdentification();
+ MpdParticleIdentification identificator;
+
+ for (Int_t i = 0; i < nReco; i++) {
+ // check clone track into ECT
+ if (nEctReco > 0) {
+ trackSet::iterator it = EctTrackSet.find(i);
+ if (it != EctTrackSet.end()) continue;
+ }
+
+ MpdKalmanTrack *kftrack = (MpdKalmanTrack*) fKFTracks->UncheckedAt(i);
+
+ MpdTrack *track = fEvent->AddGlobalTrack();
+ track->SetID(kftrack->GetTrackID());
+ track->SetNofHits(kftrack->GetNofHits());
+ track->SetdEdXTPC(kftrack->GetPartID());
+ Float_t Ppi, Pk, Pe, Pp;
+
+ if (!identificator.GetTpcProbs(kftrack->Momentum3().Mag(), kftrack->GetPartID(), kftrack->GetNofHits(), Ppi, Pk, Pp, Pe, 0)) { //0 - equal bayesian coefficients
+ track->SetTPCpidProb(Pe, Ppi, Pk, Pp, BIT(2));
+ }
+ matchingDataExist = false;
+ for (Int_t tofIndex = 0; tofIndex < nMatching; tofIndex++) {
+ pMatchingData = (MpdTofMatchingData*) fTofMatching->UncheckedAt(tofIndex);
+ if (pMatchingData->GetKFTrackIndex() == i) {
+ matchingDataExist = true;
+ break;
+ } // first matching
+ }
+
+ if (matchingDataExist) {
+ track->SetTofBeta(pMatchingData->GetBeta());
+ track->SetTofMass2(pMatchingData->GetMass2());
+ track->SetTofHitIndex(pMatchingData->GetTofHitIndex());
+
+ if (!identificator.GetTofProbs(pMatchingData->GetMomentum().Mag(), pMatchingData->GetBeta(), Ppi, Pk, Pp, Pe, 0)) {
+ track->SetTOFpidProb(Pe, Ppi, Pk, Pp, BIT(1));
+ }
+ }
+ Float_t tpcProbs[4] = {track->GetTPCPidProbPion(), track->GetTPCPidProbKaon(), track->GetTPCPidProbProton(), track->GetTPCPidProbElectron()};
+ Float_t tofProbs[4] = {track->GetTOFPidProbPion(), track->GetTOFPidProbKaon(), track->GetTOFPidProbProton(), track->GetTOFPidProbElectron()};
+ Float_t combProbs[4]; //probabilities combined from TOF & TPC
+ identificator.GetCombinedProbs(tofProbs, tpcProbs, combProbs, 4);
+ Ppi = combProbs[0];
+ Pk = combProbs[1];
+ Pp = combProbs[2];
+ Pe = combProbs[3];
+ track->SetCombPidProb(Pe, Ppi, Pk, Pp);
+
+ if (kftrack->GetParam(4) == 0.) track->SetPt(TMath::Sqrt(-1)); /*NaN*/
+ else track->SetPt(1. / kftrack->GetParam(4)); /*signed Pt*/
+
+ track->SetTheta(TMath::PiOver2() - kftrack->GetParam(3)); // Theta: angle from beam line
+ track->SetPhi(kftrack->GetParam(2)); // Phi
+
+ TMatrixD Cov = *kftrack->GetCovariance(); // Error matrix
+ track->SetPtError(TMath::Sqrt(Cov(4, 4))); // Pt error
+ track->SetThetaError(TMath::Sqrt(Cov(3, 3))); // Theta error
+ track->SetPhiError(TMath::Sqrt(Cov(2, 2))); // Phi error
+
+ track->SetChi2(kftrack->GetChi2());
+
+ Double_t phi = kftrack->GetParam(0) / kftrack->GetPosNew();
+ track->SetFirstPointX(kftrack->GetPosNew() * TMath::Cos(phi)); // closest to beam line
+ track->SetFirstPointY(kftrack->GetPosNew() * TMath::Sin(phi));
+ track->SetFirstPointZ(kftrack->GetParam(1));
+
+ track->SetLastPointX(0.); // AZ - currently not available
+ track->SetLastPointY(0.); // AZ - currently not available
+ track->SetLastPointZ(0.); // AZ - currently not available
+ }
+
+
+ Int_t nEMatching = fEtofMatching ? fEtofMatching->GetEntries() : 0;
+ MpdEtofMatchingData *pEMatchingData;
+
+ for (Int_t i = 0; i < nEctReco; i++) {
+ MpdKalmanTrack *kftrack = (MpdKalmanTrack*) fKFEctTracks->UncheckedAt(i);
+
+ MpdTrack *track = fEvent->AddGlobalTrack();
+ track->SetID(kftrack->GetTrackID());
+ track->SetNofHits(kftrack->GetNofHits());
+ matchingDataExist = false;
+ for (Int_t etofIndex = 0; etofIndex < nEMatching; etofIndex++) {
+ pEMatchingData = (MpdEtofMatchingData*) fEtofMatching->UncheckedAt(etofIndex);
+ if (pEMatchingData->GetKFTrackIndex() == i) {
+ matchingDataExist = true;
+ break;
+ } // first matching
+ }
+
+ if (matchingDataExist) {
+ track->SetTofMass2(pEMatchingData->GetMass2());
+ track->SetTofHitIndex(pEMatchingData->GetTofHitIndex());
+ }
+
+ if (kftrack->GetParam(4) == 0.) track->SetPt(TMath::Sqrt(-1)); /*NaN*/
+ else track->SetPt(1. / kftrack->GetParam(4)); /*signed Pt*/
+
+ track->SetTheta(TMath::PiOver2() - kftrack->GetParam(3)); // Theta: angle from beam line
+ track->SetPhi(kftrack->GetParam(2)); // Phi
+
+ TMatrixD Cov = *kftrack->GetCovariance(); // Error matrix
+ track->SetPtError(TMath::Sqrt(Cov(4, 4))); // Pt error
+ track->SetThetaError(TMath::Sqrt(Cov(3, 3))); // Theta error
+ track->SetPhiError(TMath::Sqrt(Cov(2, 2))); // Phi error
+
+ track->SetChi2(kftrack->GetChi2());
+
+ Double_t phi = kftrack->GetParam(0) / kftrack->GetPosNew();
+ track->SetFirstPointX(kftrack->GetPosNew() * TMath::Cos(phi)); // closest to beam line
+ track->SetFirstPointY(kftrack->GetPosNew() * TMath::Sin(phi));
+ track->SetFirstPointZ(kftrack->GetParam(1));
+
+ track->SetLastPointX(0.); // AZ - currently not available
+ track->SetLastPointY(0.); // AZ - currently not available
+ track->SetLastPointZ(0.); // AZ - currently not available
}
- }
-
- Int_t nReco = fKFTracks ? fKFTracks->GetEntriesFast() : 0;
- Int_t nEctReco = fKFEctTracks ? fKFEctTracks->GetEntriesFast() : 0;
- cout<<"\n-I- [MpdFillDstTask::Exec] "<< nReco + nEctReco<< " reconstruced tracks to write" <<endl;
-
- FairRunAna *fRun = FairRunAna::Instance();
- fEvent->SetRunInfoRunId(fRun->GetRunId());
-
- FairField *fPar = fRun->GetField();
- fEvent->SetRunInfoMagneticFieldZ(fPar->GetType());
-
- if (fVertex)
- {
- MpdVertex *vertex = (MpdVertex*) fVertex->UncheckedAt(0);
- fEvent->SetPrimaryVerticesX(vertex->GetX());
- fEvent->SetPrimaryVerticesY(vertex->GetY());
- fEvent->SetPrimaryVerticesZ(vertex->GetZ());
- } else {
- // Vertex info is not available
- fEvent->SetPrimaryVerticesX(0.);
- fEvent->SetPrimaryVerticesY(0.);
- fEvent->SetPrimaryVerticesZ(0.);
- }
-
- // check clone track into ECT
- typedef std::set<Int_t> trackSet;
- trackSet EctTrackSet;
- for(Int_t index = 0; index < nEctReco; index++) // cycle by ECT KF tracks
- {
- MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fKFEctTracks->UncheckedAt(index);
- if(track->IsFromTpc()) EctTrackSet.insert(trackSet::value_type(track->GetTpcIndex()));
- }
-
- Int_t nMatching = fTofMatching ? fTofMatching->GetEntriesFast() : 0;
- MpdTofMatchingData *pMatchingData;
- bool matchingDataExist;
-
- //AZ MpdParticleIdentification *identificator = new MpdParticleIdentification();
- MpdParticleIdentification identificator;
-
- for(Int_t i = 0; i < nReco; i++)
- {
- // check clone track into ECT
- if(nEctReco>0)
- {
- trackSet::iterator it = EctTrackSet.find(i);
- if(it != EctTrackSet.end()) continue;
- }
-
- MpdKalmanTrack *kftrack = (MpdKalmanTrack*) fKFTracks->UncheckedAt(i);
-
- MpdTrack *track = fEvent->AddGlobalTrack();
- track->SetID(kftrack->GetTrackID());
- track->SetNofHits(kftrack->GetNofHits());
- track->SetdEdXTPC(kftrack->GetPartID());
- Float_t Ppi, Pk, Pe, Pp;
-
- if(!identificator.GetTpcProbs(kftrack->Momentum3().Mag(), kftrack->GetPartID(), kftrack->GetNofHits(), Ppi, Pk, Pp, Pe, 1)) {
- track->SetTPCpidProb(Pe, Ppi, Pk, Pp, BIT(2));
- }
- matchingDataExist = false;
- for (Int_t tofIndex = 0; tofIndex < nMatching; tofIndex++)
- {
- pMatchingData = (MpdTofMatchingData*) fTofMatching->UncheckedAt(tofIndex);
- if(pMatchingData->GetKFTrackIndex() == i){ matchingDataExist = true; break;} // first matching
- }
-
- if(matchingDataExist)
- {
- track->SetTofBeta(pMatchingData->GetBeta());
- track->SetTofMass2(pMatchingData->GetMass2());
- track->SetTofHitIndex(pMatchingData->GetTofHitIndex());
-
- if(!identificator.GetTofProbs(pMatchingData->GetMomentum().Mag(), pMatchingData->GetBeta(), Ppi, Pk, Pp, Pe, 0)) {
- track->SetTOFpidProb(Pe, Ppi, Pk, Pp, BIT(1));
- }
- }
- Float_t tpcProbs[4] = {track->GetTPCPidProbPion(), track->GetTPCPidProbKaon(), track->GetTPCPidProbProton(), track->GetTPCPidProbElectron()};
- Float_t tofProbs[4] = {track->GetTOFPidProbPion(), track->GetTOFPidProbKaon(), track->GetTOFPidProbProton(), track->GetTOFPidProbElectron()};
- Float_t combProbs[4]; //probabilities combined from TOF & TPC
- identificator.GetCombinedProbs(tofProbs, tpcProbs, combProbs, 4);
- Ppi = combProbs[0];
- Pk = combProbs[1];
- Pp = combProbs[2];
- Pe = combProbs[3];
- track->SetCombPidProb(Pe, Ppi, Pk, Pp);
-
- if ( kftrack->GetParam(4) == 0.) track->SetPt(TMath::Sqrt(-1)); /*NaN*/
- else track->SetPt(1./kftrack->GetParam(4)); /*signed Pt*/
-
- track->SetTheta(TMath::PiOver2()-kftrack->GetParam(3)); // Theta: angle from beam line
- track->SetPhi(kftrack->GetParam(2)); // Phi
-
- TMatrixD Cov = *kftrack->GetCovariance(); // Error matrix
- track->SetPtError(TMath::Sqrt(Cov(4,4))); // Pt error
- track->SetThetaError(TMath::Sqrt(Cov(3,3))); // Theta error
- track->SetPhiError(TMath::Sqrt(Cov(2,2))); // Phi error
-
- track->SetChi2(kftrack->GetChi2());
-
- Double_t phi = kftrack->GetParam(0) / kftrack->GetPosNew();
- track->SetFirstPointX(kftrack->GetPosNew()*TMath::Cos(phi)); // closest to beam line
- track->SetFirstPointY(kftrack->GetPosNew()*TMath::Sin(phi));
- track->SetFirstPointZ(kftrack->GetParam(1));
-
- track->SetLastPointX(0.); // AZ - currently not available
- track->SetLastPointY(0.); // AZ - currently not available
- track->SetLastPointZ(0.); // AZ - currently not available
- }
-
-
- Int_t nEMatching = fEtofMatching ? fEtofMatching->GetEntries() : 0;
- MpdEtofMatchingData *pEMatchingData;
-
- for (Int_t i = 0; i < nEctReco; i++)
- {
- MpdKalmanTrack *kftrack = (MpdKalmanTrack*) fKFEctTracks->UncheckedAt(i);
-
- MpdTrack *track = fEvent->AddGlobalTrack();
- track->SetID(kftrack->GetTrackID());
- track->SetNofHits(kftrack->GetNofHits());
- matchingDataExist = false;
- for (Int_t etofIndex = 0; etofIndex < nEMatching; etofIndex++)
- {
- pEMatchingData = (MpdEtofMatchingData*) fEtofMatching->UncheckedAt(etofIndex);
- if(pEMatchingData->GetKFTrackIndex() == i){ matchingDataExist = true; break;} // first matching
- }
-
- if(matchingDataExist)
- {
- track->SetTofMass2(pEMatchingData->GetMass2());
- track->SetTofHitIndex(pEMatchingData->GetTofHitIndex());
- }
-
- if ( kftrack->GetParam(4) == 0.) track->SetPt(TMath::Sqrt(-1)); /*NaN*/
- else track->SetPt(1./kftrack->GetParam(4)); /*signed Pt*/
-
- track->SetTheta(TMath::PiOver2()-kftrack->GetParam(3)); // Theta: angle from beam line
- track->SetPhi(kftrack->GetParam(2)); // Phi
-
- TMatrixD Cov = *kftrack->GetCovariance(); // Error matrix
- track->SetPtError(TMath::Sqrt(Cov(4,4))); // Pt error
- track->SetThetaError(TMath::Sqrt(Cov(3,3))); // Theta error
- track->SetPhiError(TMath::Sqrt(Cov(2,2))); // Phi error
-
- track->SetChi2(kftrack->GetChi2());
-
- Double_t phi = kftrack->GetParam(0) / kftrack->GetPosNew();
- track->SetFirstPointX(kftrack->GetPosNew()*TMath::Cos(phi)); // closest to beam line
- track->SetFirstPointY(kftrack->GetPosNew()*TMath::Sin(phi));
- track->SetFirstPointZ(kftrack->GetParam(1));
-
- track->SetLastPointX(0.); // AZ - currently not available
- track->SetLastPointY(0.); // AZ - currently not available
- track->SetLastPointZ(0.); // AZ - currently not available
- }
}
// -------------------------------------------------------------------
//Delete MC tracks being outside the MPD
-void MpdFillDstTask::CleanMC(){
- cout<<"-I- [MpdFillDstTask::Exec] Cleaning from outer decays..."<<flush;
+
+void MpdFillDstTask::CleanMC() {
+ cout << "-I- [MpdFillDstTask::Exec] Cleaning from outer decays..." << flush;
Int_t nMC = fMCTracks->GetEntriesFast(), motherId;
FairMCTrack *track;
@@ -272,24 +271,21 @@ void MpdFillDstTask::CleanMC(){
vector<int> removedMother;
Int_t i, j;
- for (i = 0; i < nMC; i++)
- {
+ for (i = 0; i < nMC; i++) {
track = (FairMCTrack*) fMCTracks->UncheckedAt(i);
//filling hostograms
motherId = track->GetMotherId();
//fhTrackMotherId->Fill(motherId);
- if (motherId == -1){
+ if (motherId == -1) {
//fhTrackPrimaryPDG->Fill(track->GetPdgCode());
- }
- else{
+ } else {
//fhTrackVertex->Fill(track->GetStartX(),track->GetStartY());
//if decay was out of MPD
- if (sqrt(track->GetStartX()*track->GetStartX()+track->GetStartY()*track->GetStartY()) > rMax){
+ if (sqrt(track->GetStartX() * track->GetStartX() + track->GetStartY() * track->GetStartY()) > rMax) {
fMCTracks->Remove(track);
removedMother.push_back(i);
- }
- else{
+ } else {
//motherId < childId (i) REQUIRED
/*for (j = 0; j < cntRemovedMother; j++){
if (motherId == removedMother[j]){
@@ -320,7 +316,7 @@ void MpdFillDstTask::CleanMC(){
cntRemovedMother++;
}*/
//binary search
- if (binary_search (removedMother.begin(), removedMother.end(), motherId)){
+ if (binary_search(removedMother.begin(), removedMother.end(), motherId)) {
fMCTracks->Remove(track);
removedMother.push_back(i);
}
@@ -331,35 +327,36 @@ void MpdFillDstTask::CleanMC(){
}
fMCTracks->Compress();
- cout<<endl;
+ cout << endl;
}
// -------------------------------------------------------------------
-void MpdFillDstTask::Reset()
-{
+
+void MpdFillDstTask::Reset() {
fEvent->Reset();
}
// -------------------------------------------------------------------
-void MpdFillDstTask::Finish()
-{
- //cout<<"\n-I- [MpdFillDstTask::Finish] "<< endl;
- //fEvents->Dump();
- //cout << "\n";
-
- /* !this code gives some problems! -- !add extra event and so segfault!
- FairRootManager *fManager = FairRootManager::Instance();
- fManager->Fill();
- */
- //fhTrackMotherId->Write();
- //fhTrackPrimaryPDG->Write();
- //fhTrackVertex->Write();
- //fhTruthVertex->Write();
-
- delete fEvent;
- fEvent = NULL;
+
+void MpdFillDstTask::Finish() {
+ //cout<<"\n-I- [MpdFillDstTask::Finish] "<< endl;
+ //fEvents->Dump();
+ //cout << "\n";
+
+ /* !this code gives some problems! -- !add extra event and so segfault!
+ FairRootManager *fManager = FairRootManager::Instance();
+ fManager->Fill();
+ */
+ //fhTrackMotherId->Write();
+ //fhTrackPrimaryPDG->Write();
+ //fhTrackVertex->Write();
+ //fhTruthVertex->Write();
+
+ delete fEvent;
+ fEvent = NULL;
}
// -------------------------------------------------------------------
+
/*MpdEvent *MpdFillDstTask::AddEvent(Option_t * option)
{
TClonesArray &events = *fEvents;
@@ -368,7 +365,7 @@ void MpdFillDstTask::Finish()
return event;
}*/
-MpdTrack *MpdFillDstTask::AddPrimaryTrack(){
+MpdTrack *MpdFillDstTask::AddPrimaryTrack() {
return NULL;
}
// -------------------------------------------------------------------
diff --git a/tpc/CMakeLists.txt b/tpc/CMakeLists.txt
index bfd935bd02a93106af38f8219813fcb326c8bcce..63dbdc24ad0ee2541747dab201665ec27da4cd87 100644
--- a/tpc/CMakeLists.txt
+++ b/tpc/CMakeLists.txt
@@ -50,6 +50,7 @@ MpdTpcHitProducer.cxx
MpdTpcSectorGeo.cxx
MpdTpcDigitizerAZ.cxx
MpdTpcClusterFinderAZ.cxx
+MpdTPCpid.cxx
)
# fill list of header files from list of source files
diff --git a/tpc/MpdParticleIdentification.cxx b/tpc/MpdParticleIdentification.cxx
index 34439f1f9219e72e5f1c79814fc771a16f9cdeb1..554a64a42b957baa563caf66600edc881302e0fc 100644
--- a/tpc/MpdParticleIdentification.cxx
+++ b/tpc/MpdParticleIdentification.cxx
@@ -32,361 +32,6 @@ MpdParticleIdentification::MpdParticleIdentification() {
MpdParticleIdentification::~MpdParticleIdentification() {
}
-Int_t MpdParticleIdentification::GetTpcProbs(Float_t P, Float_t dedx, Int_t nHits, Float_t& Ppi, Float_t& Pk, Float_t& Pp, Float_t& Pe, Int_t method) {
-
- /*Parameters for fitting*/
- //AZ Float_t *ProtonPar, *PionPar, *KaonPar, *ElectronPar;
- /************************/
-
- dedx *= 1000000; // converting to keV/cm
- const Int_t Ntypes = 4; //order: pion, kaon, proton, electron
- Float_t ProtonPar[9], PionPar[9], KaonPar[9], ElectronPar[9], resultProb[Ntypes]; //AZ
- const Int_t Nintervals = 10;
- Float_t sigmas[Ntypes][Nintervals] = {//array of sigmas for different particles and different intervals of momentum. They were got from 100000 BOX events
- {0.113, 0.107, 0.100, 0.107, 0.097, 0.111, 0.122, 0.105, 0.115, 0.118},
- {0.407, 0.276, 0.260, 0.159, 0.185, 0.154, 0.122, 0.105, 0.115, 0.118},
- {0.507, 0.488, 0.394, 0.330, 0.268, 0.244, 0.260, 0.172, 0.193, 0.118},
- {0.163, 0.160, 0.149, 0.159, 0.172, 0.205, 0.260, 0.172, 0.193, 0.156}
- };
-
- Float_t sigma = 1.0 / TMath::Sqrt(nHits); //for gaussian
- Float_t sum = 0.0; // for normalizing
-
- /*A priori coefficients for Bayes approach */
- Float_t bayesAprioriCoefficients[Ntypes];
- /*A "measured" probabilities */
- Float_t gausProb[Ntypes];
-
- Float_t sigPi, sigPr, sigKa, sigEl; //sigmas for gaussians
-
- switch (method) {
-
- case 0: //bethe-bloch approximation with "standard" sigmas and equal bayesian coefficients
-
- //parameters were got from Bethe-Bloch approximation for 100000 BOX events
- /*AZ
- ProtonPar = new Float_t[4];
- PionPar = new Float_t[4];
- KaonPar = new Float_t[4];
- ElectronPar = new Float_t[4];
- */
- ProtonPar[0] = -3.957;
- ProtonPar[1] = 3.525;
- ProtonPar[2] = 16.468;
- ProtonPar[3] = 0.815;
- ProtonPar[4] = 5.247;
- PionPar[0] = -0.739;
- PionPar[1] = 7.373;
- PionPar[2] = 3904.790;
- PionPar[3] = 0.274;
- PionPar[4] = 5.497;
- KaonPar[0] = -2.590;
- KaonPar[1] = 4.918;
- KaonPar[2] = 79.722;
- KaonPar[3] = 0.357;
- KaonPar[4] = 4.511;
- ElectronPar[0] = -1.552;
- ElectronPar[1] = 1.748;
- ElectronPar[2] = 7.425;
- ElectronPar[3] = 0.980;
- ElectronPar[4] = 1.604;
-
- sigma = 0.07 * dedx;
- sigPi = sigPr = sigKa = sigEl = sigma;
- bayesAprioriCoefficients[0] = 1.0;
- bayesAprioriCoefficients[1] = 1.0;
- bayesAprioriCoefficients[2] = 1.0;
- bayesAprioriCoefficients[3] = 1.0;
-
- gausProb[0] = Gaus(dedx, BetheBlochFunction(P, PionPar), sigPi, kTRUE);
- gausProb[1] = Gaus(dedx, BetheBlochFunction(P, KaonPar), sigKa, kTRUE);
- gausProb[2] = Gaus(dedx, BetheBlochFunction(P, ProtonPar), sigPr, kTRUE);
- gausProb[3] = Gaus(dedx, BetheBlochFunction(P, ElectronPar), sigEl, kTRUE);
-
- sum = gausProb[0] + gausProb[1] + gausProb[2] + gausProb[3];
- if (sum == 0.) return 1;
- gausProb[0] /= sum;
- gausProb[1] /= sum;
- gausProb[2] /= sum;
- gausProb[3] /= sum;
- break;
-
- case 1: //bethe-bloch approximation with special different sigmas and byesian coefficients
-
- //parameters were got from Bethe-Bloch approximation for 100000 BOX events
- /*AZ
- ProtonPar = new Float_t[4];
- PionPar = new Float_t[4];
- KaonPar = new Float_t[4];
- ElectronPar = new Float_t[4];
- */
- ProtonPar[0] = -3.957;
- ProtonPar[1] = 3.525;
- ProtonPar[2] = 16.468;
- ProtonPar[3] = 0.815;
- ProtonPar[4] = 5.247;
- PionPar[0] = -0.739;
- PionPar[1] = 7.373;
- PionPar[2] = 3904.790;
- PionPar[3] = 0.274;
- PionPar[4] = 5.497;
- KaonPar[0] = -2.590;
- KaonPar[1] = 4.918;
- KaonPar[2] = 79.722;
- KaonPar[3] = 0.357;
- KaonPar[4] = 4.511;
- ElectronPar[0] = -1.552;
- ElectronPar[1] = 1.748;
- ElectronPar[2] = 7.425;
- ElectronPar[3] = 0.980;
- ElectronPar[4] = 1.604;
-
- if (P < 0.3) {
- bayesAprioriCoefficients[0] = 0.28;
- bayesAprioriCoefficients[1] = 0.25;
- bayesAprioriCoefficients[2] = 0.26;
- bayesAprioriCoefficients[3] = 0.21;
- sigPi = sigmas[0][0];
- sigKa = sigmas[1][0];
- sigPr = sigmas[2][0];
- sigEl = sigmas[3][0];
- } else if ((P >= 0.3) && (P < 0.4)) {
- bayesAprioriCoefficients[0] = 0.91;
- bayesAprioriCoefficients[1] = 0.02;
- bayesAprioriCoefficients[2] = 0.06;
- bayesAprioriCoefficients[3] = 0.01;
- sigPi = sigmas[0][1];
- sigKa = sigmas[1][1];
- sigPr = sigmas[2][1];
- sigEl = sigmas[3][1];
- } else if ((P >= 0.4) && (P < 0.5)) {
- bayesAprioriCoefficients[0] = 0.70;
- bayesAprioriCoefficients[1] = 0.11;
- bayesAprioriCoefficients[2] = 0.13;
- bayesAprioriCoefficients[3] = 0.06;
- sigPi = sigmas[0][2];
- sigKa = sigmas[1][2];
- sigPr = sigmas[2][2];
- sigEl = sigmas[3][2];
- } else if ((P >= 0.5) && (P < 0.6)) {
- bayesAprioriCoefficients[0] = 0.39;
- bayesAprioriCoefficients[1] = 0.19;
- bayesAprioriCoefficients[2] = 0.32;
- bayesAprioriCoefficients[3] = 0.10;
- sigPi = sigmas[0][3];
- sigKa = sigmas[1][3];
- sigPr = sigmas[2][3];
- sigEl = sigmas[3][3];
- } else if ((P >= 0.6) && (P < 0.7)) {
- bayesAprioriCoefficients[0] = 0.41;
- bayesAprioriCoefficients[1] = 0.17;
- bayesAprioriCoefficients[2] = 0.37;
- bayesAprioriCoefficients[3] = 0.5;
- sigPi = sigmas[0][4];
- sigKa = sigmas[1][4];
- sigPr = sigmas[2][4];
- sigEl = sigmas[3][4];
- } else if ((P >= 0.7) && (P < 0.8)) {
- bayesAprioriCoefficients[0] = 0.43;
- bayesAprioriCoefficients[1] = 0.16;
- bayesAprioriCoefficients[2] = 0.31;
- bayesAprioriCoefficients[3] = 0.10;
- sigPi = sigmas[0][5];
- sigKa = sigmas[1][5];
- sigPr = sigmas[2][5];
- sigEl = sigmas[3][5];
- } else if ((P >= 0.8) && (P < 0.9)) {
- bayesAprioriCoefficients[0] = 0.44;
- bayesAprioriCoefficients[1] = 0.11;
- bayesAprioriCoefficients[2] = 0.43;
- bayesAprioriCoefficients[3] = 0.02;
- sigPi = sigmas[0][6];
- sigKa = sigmas[1][6];
- sigPr = sigmas[2][6];
- sigEl = sigmas[3][6];
- } else if ((P >= 0.9) && (P < 1.0)) {
- bayesAprioriCoefficients[0] = 0.36;
- bayesAprioriCoefficients[1] = 0.21;
- bayesAprioriCoefficients[2] = 0.36;
- bayesAprioriCoefficients[3] = 0.07;
- sigPi = sigmas[0][7];
- sigKa = sigmas[1][7];
- sigPr = sigmas[2][7];
- sigEl = sigmas[3][7];
- } else if ((P >= 1.0) && (P < 1.2)) {
- bayesAprioriCoefficients[0] = 0.32;
- bayesAprioriCoefficients[1] = 0.32;
- bayesAprioriCoefficients[2] = 0.32;
- bayesAprioriCoefficients[3] = 0.04;
- sigPi = sigmas[0][8];
- sigKa = sigmas[1][8];
- sigPr = sigmas[2][8];
- sigEl = sigmas[3][8];
- } else if (P >= 1.2) {
- bayesAprioriCoefficients[0] = 0.34;
- bayesAprioriCoefficients[1] = 0.27;
- bayesAprioriCoefficients[2] = 0.31;
- bayesAprioriCoefficients[3] = 0.08;
- sigPi = sigmas[0][9];
- sigKa = sigmas[1][9];
- sigPr = sigmas[2][9];
- sigEl = sigmas[3][9];
- }
-
- gausProb[0] = Gaus(dedx, BetheBlochFunction(P, PionPar), sigPi, kTRUE);
- gausProb[1] = Gaus(dedx, BetheBlochFunction(P, KaonPar), sigKa, kTRUE);
- gausProb[2] = Gaus(dedx, BetheBlochFunction(P, ProtonPar), sigPr, kTRUE);
- gausProb[3] = Gaus(dedx, BetheBlochFunction(P, ElectronPar), sigEl, kTRUE);
-
- sum = gausProb[0] + gausProb[1] + gausProb[2] + gausProb[3];
- if (sum == 0.) return 1;
- gausProb[0] /= sum;
- gausProb[1] /= sum;
- gausProb[2] /= sum;
- gausProb[3] /= sum;
- break;
-
- case 2: //bethe-bloch approximation with "standard" sigmas and different byesian coefficients
-
- //parameters were got from Bethe-Bloch approximation for 100000 BOX events
- /*AZ
- ProtonPar = new Float_t[4];
- PionPar = new Float_t[4];
- KaonPar = new Float_t[4];
- ElectronPar = new Float_t[4];
- */
- ProtonPar[0] = -3.957;
- ProtonPar[1] = 3.525;
- ProtonPar[2] = 16.468;
- ProtonPar[3] = 0.815;
- ProtonPar[4] = 5.247;
- PionPar[0] = -0.739;
- PionPar[1] = 7.373;
- PionPar[2] = 3904.790;
- PionPar[3] = 0.274;
- PionPar[4] = 5.497;
- KaonPar[0] = -2.590;
- KaonPar[1] = 4.918;
- KaonPar[2] = 79.722;
- KaonPar[3] = 0.357;
- KaonPar[4] = 4.511;
- ElectronPar[0] = -1.552;
- ElectronPar[1] = 1.748;
- ElectronPar[2] = 7.425;
- ElectronPar[3] = 0.980;
- ElectronPar[4] = 1.604;
-
- sigma = 0.07 * dedx;
- sigPi = sigPr = sigKa = sigEl = sigma;
-
- if (P < 0.3) {
- bayesAprioriCoefficients[0] = 0.28;
- bayesAprioriCoefficients[1] = 0.25;
- bayesAprioriCoefficients[2] = 0.26;
- bayesAprioriCoefficients[3] = 0.21;
- } else if ((P >= 0.3) && (P < 0.4)) {
- bayesAprioriCoefficients[0] = 0.91;
- bayesAprioriCoefficients[1] = 0.02;
- bayesAprioriCoefficients[2] = 0.06;
- bayesAprioriCoefficients[3] = 0.01;
- } else if ((P >= 0.4) && (P < 0.5)) {
- bayesAprioriCoefficients[0] = 0.70;
- bayesAprioriCoefficients[1] = 0.11;
- bayesAprioriCoefficients[2] = 0.13;
- bayesAprioriCoefficients[3] = 0.06;
- } else if ((P >= 0.5) && (P < 0.6)) {
- bayesAprioriCoefficients[0] = 0.39;
- bayesAprioriCoefficients[1] = 0.19;
- bayesAprioriCoefficients[2] = 0.32;
- bayesAprioriCoefficients[3] = 0.10;
- } else if ((P >= 0.6) && (P < 0.7)) {
- bayesAprioriCoefficients[0] = 0.41;
- bayesAprioriCoefficients[1] = 0.17;
- bayesAprioriCoefficients[2] = 0.37;
- bayesAprioriCoefficients[3] = 0.5;
- } else if ((P >= 0.7) && (P < 0.8)) {
- bayesAprioriCoefficients[0] = 0.43;
- bayesAprioriCoefficients[1] = 0.16;
- bayesAprioriCoefficients[2] = 0.31;
- bayesAprioriCoefficients[3] = 0.10;
- } else if ((P >= 0.8) && (P < 0.9)) {
- bayesAprioriCoefficients[0] = 0.44;
- bayesAprioriCoefficients[1] = 0.11;
- bayesAprioriCoefficients[2] = 0.43;
- bayesAprioriCoefficients[3] = 0.02;
- } else if ((P >= 0.9) && (P < 1.0)) {
- bayesAprioriCoefficients[0] = 0.36;
- bayesAprioriCoefficients[1] = 0.21;
- bayesAprioriCoefficients[2] = 0.36;
- bayesAprioriCoefficients[3] = 0.07;
- } else if ((P >= 1.0) && (P < 1.2)) {
- bayesAprioriCoefficients[0] = 0.32;
- bayesAprioriCoefficients[1] = 0.32;
- bayesAprioriCoefficients[2] = 0.32;
- bayesAprioriCoefficients[3] = 0.04;
- } else if (P >= 1.2) {
- bayesAprioriCoefficients[0] = 0.34;
- bayesAprioriCoefficients[1] = 0.27;
- bayesAprioriCoefficients[2] = 0.31;
- bayesAprioriCoefficients[3] = 0.08;
- }
-
- gausProb[0] = Gaus(dedx, BetheBlochFunction(P, PionPar), sigPi, kTRUE);
- gausProb[1] = Gaus(dedx, BetheBlochFunction(P, KaonPar), sigKa, kTRUE);
- gausProb[2] = Gaus(dedx, BetheBlochFunction(P, ProtonPar), sigPr, kTRUE);
- gausProb[3] = Gaus(dedx, BetheBlochFunction(P, ElectronPar), sigEl, kTRUE);
-
- sum = gausProb[0] + gausProb[1] + gausProb[2] + gausProb[3];
- if (sum == 0.) return 1;
- gausProb[0] /= sum;
- gausProb[1] /= sum;
- gausProb[2] /= sum;
- gausProb[3] /= sum;
- break;
-
- case 3: //parabolic approximation
-
- //parameters were got from parabolic approximation function for 2000 UrQMD events
- ProtonPar[0] = 0.031;
- ProtonPar[1] = -0.124;
- ProtonPar[2] = 1.412;
- PionPar[0] = 0.230;
- PionPar[1] = 0.088;
- PionPar[2] = 1.072;
- KaonPar[0] = 0.676;
- KaonPar[1] = 0.621;
- KaonPar[2] = 0.831;
- ElectronPar[0] = 0.000;
- ElectronPar[1] = -0.018;
- ElectronPar[2] = 2.055;
-
- Float_t invP = 1.0 / P;
-
- gausProb[0] = Gaus(dedx, ParabolicFunction(invP, PionPar), sigPi, kTRUE);
- gausProb[1] = Gaus(dedx, ParabolicFunction(invP, KaonPar), sigKa, kTRUE);
- gausProb[2] = Gaus(dedx, ParabolicFunction(invP, ProtonPar), sigPr, kTRUE);
- gausProb[3] = Gaus(dedx, ParabolicFunction(invP, ElectronPar), sigEl, kTRUE);
-
- sum = gausProb[0] + gausProb[1] + gausProb[2] + gausProb[3];
- if (sum == 0.) return 1;
- gausProb[0] /= sum;
- gausProb[1] /= sum;
- gausProb[2] /= sum;
- gausProb[3] /= sum;
-
- break;
- }
- //AZ Float_t *resultProb = new Float_t[Ntypes];
- BayesFunction(gausProb, bayesAprioriCoefficients, resultProb, Ntypes);
-
- Ppi = resultProb[0];
- Pk = resultProb[1];
- Pp = resultProb[2];
- Pe = resultProb[3];
-
- return 0;
-}
-
Int_t MpdParticleIdentification::GetTofProbs(Float_t P, Float_t beta, Float_t& Ppi, Float_t& Pk, Float_t& Pp, Float_t& Pe, Int_t method) {
const Float_t m2_proton = 0.938 * 0.938; //square of proton mass (GeV)
diff --git a/tpc/MpdParticleIdentification.h b/tpc/MpdParticleIdentification.h
index abdb60575a10c28d416df5a7e259c75b542ebd0b..c3d2fea38bc234c0badfc43b3c3968e2b942f014 100644
--- a/tpc/MpdParticleIdentification.h
+++ b/tpc/MpdParticleIdentification.h
@@ -19,15 +19,18 @@
// Base Class Headers ----------------
#include "TSystem.h"
+#include "MpdTPCpid.h"
-class MpdParticleIdentification {
+
+class MpdParticleIdentification : public MpdTPCpid {
public:
// Constructors/Destructors ---------
MpdParticleIdentification();
- ~MpdParticleIdentification();
+ virtual ~MpdParticleIdentification();
- Int_t GetTpcProbs(Float_t P, Float_t dedx, Int_t nHits, Float_t& Ppi, Float_t& PK, Float_t& Pp, Float_t& Pe, Int_t method);
+ // Int_t GetTpcProbs(Float_t P, Float_t dedx, Int_t nHits, Float_t& Ppi, Float_t& PK, Float_t& Pp, Float_t& Pe, Int_t method);
+ // The function is realized now in the MpdTPCpid-class with corrected probability coefficients
Int_t GetTofProbs(Float_t P, Float_t beta, Float_t& Ppi, Float_t& PK, Float_t& Pp, Float_t& Pe, Int_t method);
Int_t GetCombinedProbs(Float_t *tofProbs, Float_t *tpcProbs, Float_t *resultProbs, Int_t N);
diff --git a/tpc/MpdTPCpid.cxx b/tpc/MpdTPCpid.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..f7005ab25a9a86a4396f920e597756477d278e88
--- /dev/null
+++ b/tpc/MpdTPCpid.cxx
@@ -0,0 +1,285 @@
+//-----------------------------------------------------------
+// File and Version Information:
+// $Id$
+//
+// Description:
+// Implementation of class ParticleIdentification
+// see ParticleIdentification.h for details
+//
+// Environment:
+// Software developed for MPD at NICA.
+//
+// Author List:
+// Gyulnara Eyyubova
+//
+//-----------------------------------------------------------
+
+// C/C++ Headers ----------------------
+#include <iostream>
+#include <TMath.h>
+#include <TF1.h>
+#include <TH1.h>
+#include <TH2.h>
+// This Class' Header ------------------
+#include "MpdTPCpid.h"
+
+using namespace std;
+using namespace TMath;
+
+// Class Member definitions -----------
+
+MpdTPCpid::MpdTPCpid() :
+ProtonPar(),
+PionPar(),
+KaonPar(),
+ElectronPar(),
+sigmasPi(),
+sigmasPr(),
+sigmasKa(),
+sigmasEl(),
+fCoefficients(NULL) {
+ ReadTPCResponse();
+}
+
+MpdTPCpid::~MpdTPCpid() {
+ fCoefficients->Delete();
+}
+
+Int_t MpdTPCpid::GetTpcProbs(Float_t P, Float_t dedx, Int_t nHits, Float_t& Ppi, Float_t& Pk, Float_t& Pp, Float_t& Pe, Int_t method) {
+
+
+ const Int_t Ntypes = 4; //order: pion, kaon, proton, electron
+ Float_t resultProb[Ntypes];
+ const Int_t Nintervals = 10;
+
+ Double_t P_int[Nintervals + 1] = {0.1, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.2, 3};
+
+ //default sigmas
+ Float_t sigma = 0.07 * dedx;
+
+ Float_t sum = 0.0; // for normalizing
+
+ /*A priori coefficients for Bayes approach */
+ Float_t bayesAprioriCoefficients[Ntypes];
+
+ /*A "measured" probabilities */
+ Float_t gausProb[Ntypes];
+ for (Int_t i = 0; i < Ntypes; i++) gausProb[i] = 0;
+
+ Int_t i_p = 0; // momentum interval
+
+ switch (method) {
+
+ case 0: //bethe-bloch approximation equal bayesian coefficients
+
+ if (P > 3) i_p = 9;
+ for (Int_t k = 1; k < Nintervals; k++) if (P >= P_int[k] && P < P_int[k + 1]) i_p = k;
+
+
+ bayesAprioriCoefficients[0] = 1.0;
+ bayesAprioriCoefficients[1] = 1.0;
+ bayesAprioriCoefficients[2] = 1.0;
+ bayesAprioriCoefficients[3] = 1.0;
+
+ gausProb[0] = Gaus(dedx, BetheBlochFunction(P, PionPar), sigmasPi[i_p], kTRUE); // kTRUE = normilized by sqrt(2*Pi)*sigma.
+ gausProb[1] = Gaus(dedx, BetheBlochFunction(P, KaonPar), sigmasKa[i_p], kTRUE);
+ gausProb[2] = Gaus(dedx, BetheBlochFunction(P, ProtonPar), sigmasPr[i_p], kTRUE);
+ gausProb[3] = Gaus(dedx, BetheBlochFunction(P, ElectronPar), sigmasEl[i_p], kTRUE);
+
+ sum = gausProb[0] + gausProb[1] + gausProb[2] + gausProb[3];
+ if (sum == 0.) return 1;
+ gausProb[0] /= sum;
+ gausProb[1] /= sum;
+ gausProb[2] /= sum;
+ gausProb[3] /= sum;
+ break;
+
+ case 1: //bethe-bloch approximation with special byesian coefficients
+
+ if (P > 3) i_p = 9;
+ for (Int_t k = 1; k < Nintervals; k++) if (P >= P_int[k] && P < P_int[k + 1]) i_p = k;
+
+ if (P < 0.3) {
+ bayesAprioriCoefficients[0] = 0.28;
+ bayesAprioriCoefficients[1] = 0.25;
+ bayesAprioriCoefficients[2] = 0.26;
+ bayesAprioriCoefficients[3] = 0.21;
+ } else if ((P >= 0.3) && (P < 0.4)) {
+ bayesAprioriCoefficients[0] = 0.91;
+ bayesAprioriCoefficients[1] = 0.02;
+ bayesAprioriCoefficients[2] = 0.06;
+ bayesAprioriCoefficients[3] = 0.01;
+ } else if ((P >= 0.4) && (P < 0.5)) {
+ bayesAprioriCoefficients[0] = 0.70;
+ bayesAprioriCoefficients[1] = 0.11;
+ bayesAprioriCoefficients[2] = 0.13;
+ bayesAprioriCoefficients[3] = 0.06;
+ } else if ((P >= 0.5) && (P < 0.6)) {
+ bayesAprioriCoefficients[0] = 0.39;
+ bayesAprioriCoefficients[1] = 0.19;
+ bayesAprioriCoefficients[2] = 0.32;
+ bayesAprioriCoefficients[3] = 0.10;
+ } else if ((P >= 0.6) && (P < 0.7)) {
+ bayesAprioriCoefficients[0] = 0.41;
+ bayesAprioriCoefficients[1] = 0.17;
+ bayesAprioriCoefficients[2] = 0.37;
+ bayesAprioriCoefficients[3] = 0.5;
+ } else if ((P >= 0.7) && (P < 0.8)) {
+ bayesAprioriCoefficients[0] = 0.43;
+ bayesAprioriCoefficients[1] = 0.16;
+ bayesAprioriCoefficients[2] = 0.31;
+ bayesAprioriCoefficients[3] = 0.10;
+ } else if ((P >= 0.8) && (P < 0.9)) {
+ bayesAprioriCoefficients[0] = 0.44;
+ bayesAprioriCoefficients[1] = 0.11;
+ bayesAprioriCoefficients[2] = 0.43;
+ bayesAprioriCoefficients[3] = 0.02;
+ } else if ((P >= 0.9) && (P < 1.0)) {
+ bayesAprioriCoefficients[0] = 0.36;
+ bayesAprioriCoefficients[1] = 0.21;
+ bayesAprioriCoefficients[2] = 0.36;
+ bayesAprioriCoefficients[3] = 0.07;
+ } else if ((P >= 1.0) && (P < 1.2)) {
+ bayesAprioriCoefficients[0] = 0.32;
+ bayesAprioriCoefficients[1] = 0.32;
+ bayesAprioriCoefficients[2] = 0.32;
+ bayesAprioriCoefficients[3] = 0.04;
+ } else if (P >= 1.2) {
+ bayesAprioriCoefficients[0] = 0.34;
+ bayesAprioriCoefficients[1] = 0.27;
+ bayesAprioriCoefficients[2] = 0.31;
+ bayesAprioriCoefficients[3] = 0.08;
+ }
+
+ gausProb[0] = Gaus(dedx, BetheBlochFunction(P, PionPar), sigmasPi[i_p], kTRUE);
+ gausProb[1] = Gaus(dedx, BetheBlochFunction(P, KaonPar), sigmasKa[i_p], kTRUE);
+ gausProb[2] = Gaus(dedx, BetheBlochFunction(P, ProtonPar), sigmasPr[i_p], kTRUE);
+ gausProb[3] = Gaus(dedx, BetheBlochFunction(P, ElectronPar), sigmasEl[i_p], kTRUE);
+
+ sum = gausProb[0] + gausProb[1] + gausProb[2] + gausProb[3];
+ if (sum == 0.) return 1;
+ gausProb[0] /= sum;
+ gausProb[1] /= sum;
+ gausProb[2] /= sum;
+ gausProb[3] /= sum;
+ break;
+
+
+
+ }
+
+
+ BayesFunction(gausProb, bayesAprioriCoefficients, resultProb, Ntypes);
+
+ Ppi = resultProb[0];
+ Pk = resultProb[1];
+ Pp = resultProb[2];
+ Pe = resultProb[3];
+
+ return 0;
+}
+
+Float_t MpdTPCpid::BetheBlochFunction(Float_t x, Float_t *p) {
+ Float_t b = 1 / (x / Sqrt(1 + x * x));
+ return p[0] / Power(b, p[3]) * (p[1] - Power(b, p[3]) - Log(p[2] + Power(1 / x, p[4])));
+}
+
+/*Bayes function for probabilities calculation*/
+Int_t MpdTPCpid::BayesFunction(Float_t *measProb, Float_t *aprioriProb, Float_t *bayesProb, Int_t N) {
+
+ /*measProb - measured probabilities from TPC/TOF/etc*/
+ /*aprioriProb - a priori probabilities from some magic*/
+ /*bayesProb - result bayes probabilities */
+ /*N - number of types {pion, kaon, proton, electron} */
+
+ Float_t sumBayes = 0.0;
+ for (Int_t i = 0; i < N; ++i) {
+ sumBayes += measProb[i] * aprioriProb[i];
+ }
+
+ if (sumBayes == 0.) return 1;
+
+ for (Int_t i = 0; i < N; ++i) {
+ bayesProb[i] = measProb[i] * aprioriProb[i] / sumBayes;
+ }
+
+ return 0;
+}
+
+void MpdTPCpid::ReadTPCResponse() {
+ TString c;
+ /* Accessing the DeDx response, obtained with VHLEE generator */
+ fCoefficients = new TFile("$VMCWORKDIR/input/MPDTPCPidResponseVhelle.root");
+
+ if (fCoefficients->IsZombie()) {
+ printf("File MPDTPCPidResponseVhelle.root does not exist.\n");
+ }
+
+ TH2D *hdedx_pi = (TH2D*) fCoefficients->Get("dedx_pi");
+ TH2D *hdedx_pr = (TH2D*) fCoefficients->Get("dedx_pr");
+ TH2D *hdedx_ka = (TH2D*) fCoefficients->Get("dedx_ka");
+ TH2D *hdedx_el = (TH2D*) fCoefficients->Get("dedx_el");
+
+ TF1 *fBetheBlochPion = hdedx_pi->GetFunction("BetheBlochALEPH");
+ PionPar[0] = fBetheBlochPion->GetParameter(0);
+ PionPar[1] = fBetheBlochPion->GetParameter(1);
+ PionPar[2] = fBetheBlochPion->GetParameter(2);
+ PionPar[3] = fBetheBlochPion->GetParameter(3);
+ PionPar[4] = fBetheBlochPion->GetParameter(4);
+
+ TF1 *fBetheBlochProton = hdedx_pr->GetFunction("BetheBlochALEPH");
+ ProtonPar[0] = fBetheBlochProton->GetParameter(0);
+ ProtonPar[1] = fBetheBlochProton->GetParameter(1);
+ ProtonPar[2] = fBetheBlochProton->GetParameter(2);
+ ProtonPar[3] = fBetheBlochProton->GetParameter(3);
+ ProtonPar[4] = fBetheBlochProton->GetParameter(4);
+
+ TF1 *fBetheBlochKaon = hdedx_ka->GetFunction("BetheBlochALEPH");
+ KaonPar[0] = fBetheBlochKaon->GetParameter(0);
+ KaonPar[1] = fBetheBlochKaon->GetParameter(1);
+ KaonPar[2] = fBetheBlochKaon->GetParameter(2);
+ KaonPar[3] = fBetheBlochKaon->GetParameter(3);
+ KaonPar[4] = fBetheBlochKaon->GetParameter(4);
+
+ TF1 *fBetheBlochElectron = hdedx_el->GetFunction("BetheBlochALEPH");
+ ElectronPar[0] = fBetheBlochElectron->GetParameter(0);
+ ElectronPar[1] = fBetheBlochElectron->GetParameter(1);
+ ElectronPar[2] = fBetheBlochElectron->GetParameter(2);
+ ElectronPar[3] = fBetheBlochElectron->GetParameter(3);
+ ElectronPar[4] = fBetheBlochElectron->GetParameter(4);
+
+
+ const Int_t Nintervals = 10;
+ Double_t P_int[Nintervals + 1] = {0.1, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.2, 3};
+
+ for (Int_t k = 0; k < Nintervals; k++) {
+
+ Int_t p_low_bin = hdedx_pi->GetXaxis()->FindBin(P_int[k]);
+ Int_t p_up_bin = hdedx_pi->GetXaxis()->FindBin(P_int[k + 1]);
+
+ c = Form("hpi%i", k);
+ TH1D *hpi = (TH1D*) hdedx_pi->ProjectionY(c, p_low_bin, p_up_bin);
+ hpi->Fit("gaus");
+ sigmasPi[k] = hpi->GetFunction("gaus")->GetParameter(2);
+
+ c = Form("hpr%i", k);
+ TH1D *hpr = (TH1D*) hdedx_pr->ProjectionY("hpr", p_low_bin, p_up_bin);
+ hpr->Fit("gaus");
+ sigmasPr[k] = hpr->GetFunction("gaus")->GetParameter(2);
+
+ c = Form("hka%i", k);
+ TH1D *hka = (TH1D*) hdedx_ka->ProjectionY("hka", p_low_bin, p_up_bin);
+ hka->Fit("gaus");
+ sigmasKa[k] = hka->GetFunction("gaus")->GetParameter(2);
+
+ c = Form("hel%i", k);
+ TH1D *hel = (TH1D*) hdedx_el->ProjectionY("hel", p_low_bin, p_up_bin);
+ hel->Fit("gaus");
+ sigmasEl[k] = hel->GetFunction("gaus")->GetParameter(2);
+ }
+
+ // f->Close();
+
+}
+
+ClassImp(MpdTPCpid);
diff --git a/tpc/MpdTPCpid.h b/tpc/MpdTPCpid.h
new file mode 100644
index 0000000000000000000000000000000000000000..db0d3fde645d857f8fdf96470ec60cfc1eabd87f
--- /dev/null
+++ b/tpc/MpdTPCpid.h
@@ -0,0 +1,56 @@
+//-----------------------------------------------------------
+// File and Version Information:
+// $Id$
+//
+// Description:
+// class for particle identification
+//
+//
+// Environment:
+// Software developed for MPD at NICA.
+//
+// Author List:
+// Gyulnara Eyyubova
+//
+//-----------------------------------------------------------
+
+#ifndef MpdTPCpid_HH
+#define MpdTPCpid_HH
+
+// Base Class Headers ----------------
+#include "TSystem.h"
+#include <TGraph.h>
+#include <TFile.h>
+
+class MpdTPCpid {
+public:
+
+ // Constructors/Destructors ---------
+ MpdTPCpid();
+ virtual ~MpdTPCpid();
+
+ Int_t GetTpcProbs(Float_t P, Float_t dedx, Int_t nHits, Float_t& Ppi, Float_t& PK, Float_t& Pp, Float_t& Pe, Int_t method);
+ Float_t BetheBlochFunction(Float_t x, Float_t *p);
+ Int_t BayesFunction(Float_t *measProb, Float_t *aprioriProb, Float_t *bayesProb, Int_t N);
+ void ReadTPCResponse();
+
+private:
+
+ Float_t ProtonPar[5];
+ Float_t PionPar[5];
+ Float_t KaonPar[5];
+ Float_t ElectronPar[5];
+ //const Int_t Nintervals = 10;
+ Float_t sigmasPi[10];
+ Float_t sigmasPr[10];
+ Float_t sigmasKa[10];
+ Float_t sigmasEl[10];
+
+ TFile* fCoefficients;
+
+public:
+ ClassDef(MpdTPCpid, 1)
+
+};
+
+#endif
diff --git a/tpc/tpcLinkDef.h b/tpc/tpcLinkDef.h
index badcf2790a0e09ea5334586fca671b5e15f5e895..9c71ad5931529e76e22e447d0fed24ceb22c63c2 100644
--- a/tpc/tpcLinkDef.h
+++ b/tpc/tpcLinkDef.h
@@ -27,6 +27,7 @@
#pragma link C++ class MpdTpcSectorGeo+;
#pragma link C++ class MpdTpcDigitizerAZ+;
#pragma link C++ class MpdTpcClusterFinderAZ+;
+#pragma link C++ class MpdTPCpid+;
#endif