diff --git a/analysis/PWGHAD/hadron/CbmHadronAnalysis.cxx b/analysis/PWGHAD/hadron/CbmHadronAnalysis.cxx
index 09fba3f950c545c00385f91f0808a4f2976fe411..b0f5f48513a951dc5742bdd96d1a3a8a36775689 100644
--- a/analysis/PWGHAD/hadron/CbmHadronAnalysis.cxx
+++ b/analysis/PWGHAD/hadron/CbmHadronAnalysis.cxx
@@ -46,6 +46,7 @@ using namespace std;
 #include "FairMCEventHeader.h"
 #include "FairMCPoint.h"
 #include "FairRootManager.h"
+#include "FairRunAna.h"
 
 CbmDigiManager* fDigiMan;   // TOF Input Digis
 TClonesArray* fEventsColl;  // CBMEvents (time based)
@@ -68,6 +69,11 @@ Int_t nGlobTracks  = 0;
 Int_t NMASS        = 3;
 Float_t refMass[3] = {0.139, 0.494, 0.938};
 
+static TH1F* fhTofHitMul;
+static TH1F* fhStsHitMul;
+const Int_t fiTofHitMulMax = 200;
+static TFile* fHist;
+
 //___________________________________________________________________
 //
 // CbmHadronAnalysis
@@ -473,8 +479,23 @@ CbmHadronAnalysis::~CbmHadronAnalysis() {
 // ------------------------------------------------------------------
 void CbmHadronAnalysis::CreateHistogramms() {
   // Create histogramms
-  gROOT->cd();
-  LOG(info) << "CreateHistograms in " << gDirectory->GetName();
+  // gROOT->cd();
+  FairRunAna *fRun= FairRunAna::Instance();
+  fHist = fRun->GetOutputFile();
+  LOG(info) << "CreateHistograms in " << fHist->GetName();
+  // gSystem->cd(fHist->GetName());
+  fHist->ReOpen("Update");
+
+  fhTofHitMul      = new TH1F(Form("hTofHitMul"),
+                         Form("TofHit Multiplicity; M_{TofHit} "),
+                         fiTofHitMulMax,
+                         0.,
+                         (Double_t) fiTofHitMulMax);
+  fhStsHitMul      = new TH1F(Form("hStsHitMul"),
+                         Form("StsHit Multiplicity; M_{StsHit} "),
+                         fiTofHitMulMax,
+                         0.,
+                         (Double_t) fiTofHitMulMax);
 
   Float_t ymin   = -1.;
   Float_t ymax   = 4.;
@@ -3187,6 +3208,8 @@ void CbmHadronAnalysis::ExecEvent(Option_t*) {
   if (fGlobalTracks != NULL) nGlobTracks = fGlobalTracks->GetEntriesFast();
   if (fStsHits != NULL) nStsHits = fStsHits->GetEntriesFast();
 
+  fhTofHitMul->Fill((Double_t) nTofHits);
+  fhStsHitMul->Fill((Double_t) nStsHits);
 
   if (verbose > 0) {  //nh-debug
     LOG(debug) << "<D> HadronAnalysis::Exec starting with MCtrks " << nMCTracks
@@ -5752,7 +5775,7 @@ void CbmHadronAnalysis::Finish() {
 
   // Finish of the task execution
 
-  // WriteHistogramms();
+  WriteHistogramms();
 }
 // ------------------------------------------------------------------
 
@@ -5760,7 +5783,7 @@ void CbmHadronAnalysis::Finish() {
 // ------------------------------------------------------------------
 void CbmHadronAnalysis::WriteHistogramms() {
   // Write histogramms to the file
-  TFile* fHist = new TFile("data/auaumbias.hst.root", "RECREATE");
+  if (NULL != fHist)
   {
     TIter next(gDirectory->GetList());
     TH1* h;
@@ -5773,8 +5796,8 @@ void CbmHadronAnalysis::WriteHistogramms() {
       }
     }
   }
-  //fHist->ls();
-  fHist->Close();
+  fHist->ls();
+  //fHist->Close();
 }
 // ------------------------------------------------------------------
 static Int_t iCandEv = 0;
@@ -5782,8 +5805,6 @@ void CbmHadronAnalysis::ReconstructSecondaries() {
 #include "TLorentzVector.h"
 #include "TVector3.h"
 
-  static TH1F* fhTofHitMul;
-  static TH1F* fhStsHitMul;
   static TH1F* fhTofChi;
   static TH1F* fhDperp;
   static TH2F* fhdEdxMul;
@@ -5829,7 +5850,6 @@ void CbmHadronAnalysis::ReconstructSecondaries() {
   */
 
 
-  const Int_t fiTofHitMulMax = 200;
   const Int_t fiNMixClasses  = 10;
   const Double_t beamRotY    = -25.;
   const Double_t MLAM        = 1.1156;
@@ -5860,16 +5880,7 @@ void CbmHadronAnalysis::ReconstructSecondaries() {
   if (iCandEv == 0) {  //initialize
     // define some histograms
     Double_t MinvMin = secMass[0] + secMass[1];
-    fhTofHitMul      = new TH1F(Form("hTofHitMul"),
-                           Form("TofHit Multiplicity; M_{TofHit} "),
-                           fiTofHitMulMax,
-                           0.,
-                           (Double_t) fiTofHitMulMax);
-    fhStsHitMul      = new TH1F(Form("hStsHitMul"),
-                           Form("StsHit Multiplicity; M_{StsHit} "),
-                           fiTofHitMulMax,
-                           0.,
-                           (Double_t) fiTofHitMulMax);
+    cout << "Add secondary histos to " << fHist->GetName() << endl;
     fhTofChi         = new TH1F(
       Form("hTofChi"), Form("TofHit Merger; #chi "), 100, 0., dChiTofLim * 2.);
     fhDperp = new TH1F(Form("hDperp"),
@@ -6019,9 +6030,6 @@ void CbmHadronAnalysis::ReconstructSecondaries() {
     }
   }
 
-  fhTofHitMul->Fill((Double_t) nTofHits);
-  fhStsHitMul->Fill((Double_t) nStsHits);
-
   fvP.resize(fiNMixClasses);
   fvX.resize(fiNMixClasses);
   fvX0.resize(fiNMixClasses);
diff --git a/core/data/tof/CbmTofDetectorId.h b/core/data/tof/CbmTofDetectorId.h
index c0ca9f976a180a83c1ede422855576013b3a64d7..646a35594b595740878993f417020383480a6c7e 100644
--- a/core/data/tof/CbmTofDetectorId.h
+++ b/core/data/tof/CbmTofDetectorId.h
@@ -35,8 +35,8 @@ class CbmTofDetectorInfo
    : fDetectorSystem(ToIntegralType(detsystem)), fSMtype(smtype), fSModule(smodule),
     fCounter(counter), fGap(gap), fCell(cell) {};
     
-CbmTofDetectorInfo(ECbmModuleId detsystem, Int_t smtype, Int_t smodule, 
-            Int_t counter, Int_t counterType, Int_t gap, Int_t cell)
+ CbmTofDetectorInfo(ECbmModuleId detsystem, Int_t smtype, Int_t smodule,
+            Int_t counter, Int_t gap, Int_t cell, Int_t counterType)
    : fDetectorSystem(ToIntegralType(detsystem)), fSMtype(smtype), fSModule(smodule),
      fCounter(counter), fGap(gap), fCell(cell), fCounterType(counterType) {};
   
diff --git a/core/data/tof/CbmTofDetectorId_v21a.cxx b/core/data/tof/CbmTofDetectorId_v21a.cxx
index 349b8c7095c3862912e238c347dba25d7ec286dc..489318bd8d1bace48b6d22430f03e84c744ff64e 100644
--- a/core/data/tof/CbmTofDetectorId_v21a.cxx
+++ b/core/data/tof/CbmTofDetectorId_v21a.cxx
@@ -6,7 +6,7 @@
 
 #include "CbmTofDetectorId_v21a.h"
 
-const Int_t CbmTofDetectorId_v21a::shiftarray[] = {0,4,11,15,21,22,28	};
+const Int_t CbmTofDetectorId_v21a::shiftarray[] = {0,4,11,15,21,22,28};
 const Int_t CbmTofDetectorId_v21a::bitarray[]   = {4,7, 4, 6, 1, 6, 4};
 
 
@@ -28,6 +28,7 @@ CbmTofDetectorId_v21a::CbmTofDetectorId_v21a()
 		  (maskarray[5] << shiftarray[5]) |
 		  (maskarray[6] << shiftarray[6])
   );
+  std::cout << "<I> V21a module mask 0x" << std::hex << modulemask << std::endl;
 }
 
 CbmTofDetectorInfo CbmTofDetectorId_v21a::GetDetectorInfo(const Int_t detectorId)
@@ -118,6 +119,15 @@ Int_t CbmTofDetectorId_v21a::GetCellId(const Int_t detectorId)
 
 Int_t CbmTofDetectorId_v21a ::SetDetectorInfo(const CbmTofDetectorInfo detInfo)
 {
+  std::cout << "SetDetectorInfo for "
+		  << detInfo.fSMtype
+		  << detInfo.fSModule << " "
+		  << detInfo.fCounter
+		  << detInfo.fCounterType <<" "
+		  << detInfo.fGap
+		  << detInfo.fCell
+		  << std::endl;
+
   return ( (((detInfo.fDetectorSystem) & maskarray[0]) << shiftarray[0]) | 
            (((detInfo.fSMtype)         & maskarray[2]) << shiftarray[2]) | 
            (((detInfo.fSModule)        & maskarray[1]) << shiftarray[1]) | 
diff --git a/core/data/tof/CbmTofDetectorId_v21a.h b/core/data/tof/CbmTofDetectorId_v21a.h
index eff9be62b04da915f68e53ecaeb8679767086c01..6725a2007e16d800892ca005ec6c9fe14ece5416 100644
--- a/core/data/tof/CbmTofDetectorId_v21a.h
+++ b/core/data/tof/CbmTofDetectorId_v21a.h
@@ -7,6 +7,7 @@
 /** CbmTofDetectorId.h
  ** Defines unique detector identifier for all TOF modules.
  ** This class is the implementation for tof geometry version v21a 
+ ** nh, 28.07.2021 add counter type
  ** nh, 11.03.2014
  ** PAL, 23.09.2015: make the class common to both v14 and v15 geometries
  **                  Field 4 used as Side index (or fake Gap index in digitizer)
diff --git a/core/detectors/tof/CbmTofCreateDigiPar.cxx b/core/detectors/tof/CbmTofCreateDigiPar.cxx
index 671b265175571d240612de590325bb7b385d299e..c55304cec88f753e8171166b8f6ae807249682e7 100644
--- a/core/detectors/tof/CbmTofCreateDigiPar.cxx
+++ b/core/detectors/tof/CbmTofCreateDigiPar.cxx
@@ -110,7 +110,10 @@ InitStatus CbmTofCreateDigiPar::Init() {
   FairRootManager* ioman = FairRootManager::Instance();
   if (!ioman) LOG(fatal) << "No FairRootManager found";
 
-
+  if (k21a == geoVersion) {
+    LOG(info) << "Will now create digitization parameters for root geometry.";
+    FillCellMapRootGeometry();
+  }
   if (k14a == geoVersion) {
     LOG(info) << "Will now create digitization parameters for root geometry.";
     FillCellMapRootGeometry();
@@ -294,6 +297,7 @@ void CbmTofCreateDigiPar::FillCellMapRootGeometry() {
   // Example for full path to gap
   //   /cave_0/tof_v12b_0/module_0_0/gas_box_0/counter_0/Gap_0/Cell_1
 
+  //   /TOP_1/tof_v21b_mcbm_1/tof_v21b_mcbmStand_1/module_9_0/gas_box_0/counter_1/Gap_17/Cell_18
 
   /*  Int_t nrCells = 0;*/
   std::vector<CbmTofCell*> cellVector;
@@ -437,16 +441,16 @@ void CbmTofCreateDigiPar::FillCellInfoFromGeoHandler(TString FullPath) {
   fZ = fGeoHandler->GetZ(FullPath);
 
   LOG(debug2) << "FCI: " << FullPath.Data();
-  LOG(debug2) << "FCI: X: " << fX;
-  LOG(debug2) << " Y: " << fY;
-  LOG(debug2) << " Z: " << fZ;
-  LOG(debug2) << " SizeX: " << fSizex;
-  LOG(debug2) << " SizeY: " << fSizey;
-  LOG(debug2) << Form(" DetID: 0x%08x", fDetID);
-  LOG(debug2) << " Region: " << fRegion;
-  LOG(debug2) << " Module: " << fCounter;
-  LOG(debug2) << " Gap: " << fGap;
-  LOG(debug2) << " Cell: " << fCell;
+  LOG(debug2) << "FCI: X: " << fX
+  	  	  	  << " Y: " << fY
+  	  	  	  << " Z: " << fZ
+  	  	  	  << " SizeX: " << fSizex
+  	  	  	  << " SizeY: " << fSizey;
+  LOG(debug2) << Form(" DetID: 0x%08x", fDetID)
+  	  	  	  << " Region: " << fRegion
+			  << " Counter: " << fCounter
+			  << " Gap: " << fGap
+			  << " Cell: " << fCell;
 
   fCellID = fGeoHandler->GetCellId(fDetID);
 
@@ -459,12 +463,12 @@ void CbmTofCreateDigiPar::FillCellInfoFromGeoHandler(TString FullPath) {
 
   LOG(debug2) << "FCI: Cell ID: " << Form("0x%08x", fCellID) << " detId "
               << Form("0x%08x", fDetID);
-  LOG(debug2) << " Region:  " << fGeoHandler->GetRegion(fCellID);
-  LOG(debug2) << " SMTYP:   " << fGeoHandler->GetSMType(fCellID);
-  LOG(debug2) << " SModule: " << fGeoHandler->GetSModule(fCellID);
-  LOG(debug2) << " Module:  " << fGeoHandler->GetCounter(fCellID);
-  LOG(debug2) << " Gap:     " << fGeoHandler->GetGap(fCellID);
-  LOG(debug2) << " Cell: " << fGeoHandler->GetCell(fCellID);
+  LOG(debug2) << " Region:  " << fGeoHandler->GetRegion(fCellID)
+			  << " SMTYP:   " << fGeoHandler->GetSMType(fCellID)
+			  << " SModule: " << fGeoHandler->GetSModule(fCellID)
+			  << " Module:  " << fGeoHandler->GetCounter(fCellID)
+			  << " Gap:     " << fGeoHandler->GetGap(fCellID)
+			  << " Cell: " << fGeoHandler->GetCell(fCellID);
 }
 
 
diff --git a/core/detectors/tof/CbmTofGeoHandler.cxx b/core/detectors/tof/CbmTofGeoHandler.cxx
index f133d8700fb59a073465451a10bfc226075a26c3..d80311bb59f633bbbf808e672d0664f929494be2 100644
--- a/core/detectors/tof/CbmTofGeoHandler.cxx
+++ b/core/detectors/tof/CbmTofGeoHandler.cxx
@@ -162,6 +162,7 @@ Int_t CbmTofGeoHandler::GetUniqueDetectorId() {
   Int_t smtype  = 0;
   Int_t smodule = 0;
   Int_t counter = 0;
+  Int_t countertype = 0;
   Int_t gap     = 0;
   Int_t cell    = 0;
   TString Volname;
@@ -195,20 +196,40 @@ Int_t CbmTofGeoHandler::GetUniqueDetectorId() {
     //    smodule=smtype;   // for test beam setup
     gap = 0;
     cell--;
+  } else if (fGeoVersion == k21a) {  // test beam
+    if (fUseNodeName) {
+      Volname = CurrentNodeOffName(4);
+    } else {
+      Volname = CurrentVolOffName(4);
+    }
+    smtype = Volname[7] - '0';
+    CurrentVolOffID(4, smodule);
+    CurrentVolOffID(2, counter);
+    CurrentVolOffID(1, gap);
+    CurrentVolID(cell);
+    //    counter=smodule;  // necessary for plastics
+    //    smodule=smtype;   // for test beam setup
+    gap = 0;
+    cell--;
   }
 
-  LOG(debug1) << "GeoHand: ";
   LOG(debug1) << " Volname: " << Volname << ", " << CurrentVolOffName(3) << ", "
               << CurrentVolOffName(2) << ", " << CurrentVolOffName(1) << ", "
               << CurrentVolOffName(0);
-  LOG(debug1) << " SMtype: " << smtype;
-  LOG(debug1) << " SModule: " << smodule;
-  LOG(debug1) << " Counter: " << counter;
-  LOG(debug1) << " Gap: " << gap;
-  LOG(debug1) << " Cell: " << cell;
+
+  TString cTemp=CurrentVolOffName(2);
+  TString cType=cTemp(8,2);  // 1 character only
+  countertype=cType.Atoi();
+
+  LOG(debug1) << " SMtype: " << smtype
+		  	  << " SModule: " << smodule
+			  << " CounterType: " << countertype
+			  << " Counter: " << counter
+			  << " Gap: " << gap
+			  << " Strip: " << cell;
 
   CbmTofDetectorInfo detInfo(
-    ECbmModuleId::kTof, smtype, smodule, counter, gap, cell);
+    ECbmModuleId::kTof, smtype, smodule, counter, gap, cell, countertype);
 
   Int_t result = fTofId->SetDetectorInfo(detInfo);
   LOG(debug1) << " Unique ID: " << Form("0x%08x", result);
@@ -221,6 +242,7 @@ Int_t CbmTofGeoHandler::GetUniqueCounterId() {
 
   Int_t smtype  = 0;
   Int_t smodule = 0;
+  Int_t countertype = 0;
   Int_t counter = 0;
   Int_t gap     = 0;
   Int_t cell    = 0;
@@ -241,39 +263,54 @@ Int_t CbmTofGeoHandler::GetUniqueCounterId() {
     CurrentVolOffID(1, gap);
     CurrentVolID(cell);
   } else if (fGeoVersion == k14a) {  // test beam
-    if (fUseNodeName) {
-      Volname = CurrentNodeOffName(4);
-    } else {
-      Volname = CurrentVolOffName(4);
-    }
-    smtype = Volname[7] - '0';
-    CurrentVolOffID(4, smodule);
-    CurrentVolOffID(2, counter);
-    CurrentVolOffID(1, gap);
-    CurrentVolID(cell);
-    //    counter=smodule;  // necessary for plastics
-    //    smodule=smtype;   // for test beam setup
+	    if (fUseNodeName) {
+	      Volname = CurrentNodeOffName(4);
+	    } else {
+	      Volname = CurrentVolOffName(4);
+	    }
+	    smtype = Volname[7] - '0';
+	    CurrentVolOffID(4, smodule);
+	    CurrentVolOffID(2, counter);
+	    CurrentVolOffID(1, gap);
+	    CurrentVolID(cell);
+	    //    counter=smodule;  // necessary for plastics
+	    //    smodule=smtype;   // for test beam setup
+  }	else if (fGeoVersion == k21a) {  // test beam
+		if (fUseNodeName) {
+		  Volname = CurrentNodeOffName(4);
+		} else {
+		  Volname = CurrentVolOffName(4);
+		}
+		smtype = Volname[7] - '0';
+		CurrentVolOffID(4, smodule);
+		TString cTemp=CurrentVolOffName(2);
+		TString cType=cTemp(8,2);  // 1 character only
+		countertype=cType.Atoi();
+		CurrentVolOffID(2, counter);
+		CurrentVolOffID(1, gap);
+		CurrentVolID(cell);
   }
 
   cell = 0;
 
   fDetectorInfoArray =
-    CbmTofDetectorInfo(ECbmModuleId::kTof, smtype, smodule, counter, gap, cell);
+    CbmTofDetectorInfo(ECbmModuleId::kTof, smtype, smodule, counter, gap, cell, countertype);
 
   gap = 0;
 
-  LOG(debug1) << "GeoHand: ";
+  LOG(debug1) << "GetUniqueCounterId: ";
   LOG(debug1) << " Volname: " << Volname << ", " << CurrentVolOffName(3) << ", "
               << CurrentVolOffName(2) << ", " << CurrentVolOffName(1) << ", "
               << CurrentVolOffName(0);
-  LOG(debug1) << " SMtype: " << smtype;
-  LOG(debug1) << " SModule: " << smodule;
-  LOG(debug1) << " Counter: " << counter;
-  LOG(debug1) << " Gap: " << gap;
-  LOG(debug1) << " Cell: " << cell;
+  LOG(debug1) << " SMtype: " << smtype
+		  	  << " SModule: " << smodule
+			  << " CounterType: " << countertype
+			  << " Counter: " << counter
+			  << " Gap: " << gap
+			  << " Cell: " << cell;
 
   CbmTofDetectorInfo detInfo(
-    ECbmModuleId::kTof, smtype, smodule, counter, gap, cell);
+    ECbmModuleId::kTof, smtype, smodule, counter, gap, cell, countertype);
 
   Int_t result        = fTofId->SetDetectorInfo(detInfo);
   fLastUsedDetectorID = result;