diff --git a/macro/beamtime/hd2020/ana_digi_cal.C b/macro/beamtime/hd2020/ana_digi_cal.C
new file mode 100644
index 0000000000000000000000000000000000000000..b348c8d074c117818db1f1143392d78a99c9ffc3
--- /dev/null
+++ b/macro/beamtime/hd2020/ana_digi_cal.C
@@ -0,0 +1,461 @@
+void ana_digi_cal(Int_t nEvents      = 10000000,
+                  Int_t calMode      = 53,
+                  Int_t calSel       = 0,
+                  Int_t calSm        = 900,
+                  Int_t RefSel       = 1,
+                  TString cFileId    = "Test",
+                  Int_t iCalSet      = 910920900,
+                  Bool_t bOut        = 0,
+                  Int_t iSel2        = 0,
+                  Double_t dDeadtime = 50,
+                  TString cCalId     = "XXX",
+                  Int_t iPlot        = 1) {
+  Int_t iVerbose = 1;
+  Int_t iBugCor  = 0;
+  //Specify log level (INFO, DEBUG, DEBUG1, ...)
+  //TString logLevel = "FATAL";
+  //TString logLevel = "ERROR";
+  TString logLevel = "INFO";
+  //TString logLevel = "DEBUG";
+  //TString logLevel = "DEBUG1";
+  //TString logLevel = "DEBUG2";
+  //TString logLevel = "DEBUG3";
+  FairLogger::GetLogger();
+  gLogger->SetLogScreenLevel(logLevel);
+  gLogger->SetLogVerbosityLevel("VERYHIGH");
+
+  TString workDir = gSystem->Getenv("VMCWORKDIR");
+  /*
+   TString workDir    = (TString)gInterpreter->ProcessLine(".! pwd");
+   cout << "workdir = "<< workDir.Data() << endl;
+   return;
+  */
+  TString paramDir = workDir + "/macro/beamtime/hd2020/";
+  //   TString ParFile    = paramDir + "data/unp_tof_" + cFileId + ".params.root";
+  // TString InputFile  = paramDir + "data/unp_tof_" + cFileId + ".root";
+  TString ParFile   = paramDir + "data/" + cFileId + ".params.root";
+  TString InputFile = paramDir + "data/" + cFileId + ".root";
+  //TString InputFile  = "/d/cbm/HD2020/" + cFileId + ".root";
+  TString OutputFile =
+    paramDir + "data/digidev_" + cFileId
+    + Form("_%09d_%03d_%02.0f_Cal", iCalSet, iSel2, dDeadtime) + cCalId
+    + ".out.root";
+
+  TString shcmd = "rm -v " + ParFile;
+  gSystem->Exec(shcmd.Data());
+
+  TList* parFileList = new TList();
+
+  TString FId = cFileId;
+  //   TString TofGeo="v18n_cosmicHD";   // for Buc2018 on top
+  //   TString TofGeo="v18m_cosmicHD";  // for Buc2015 measurement
+  //   TString TofGeo="v18o_cosmicHD";        // for Buc2018 sandwiched
+  TString TofGeo = "v20b_cosmicHD";  // for Buc2020 sandwiched
+  cout << "Geometry version " << TofGeo << endl;
+
+  TObjString* tofDigiFile = new TObjString(
+    workDir + "/parameters/tof/tof_" + TofGeo + ".digi.par");  // TOF digi file
+  // parFileList->Add(tofDigiFile);
+
+  //   TObjString tofDigiBdfFile = new TObjString( paramDir + "/tof." + FPar + "digibdf.par");
+  TObjString* tofDigiBdfFile =
+    new TObjString(workDir + "/parameters/tof/tof_" + TofGeo + ".digibdf.par");
+  parFileList->Add(tofDigiBdfFile);
+
+  TString geoDir  = gSystem->Getenv("VMCWORKDIR");
+  TString geoFile = geoDir + "/geometry/tof/geofile_tof_" + TofGeo + ".root";
+  TFile* fgeo     = new TFile(geoFile);
+  TGeoManager* geoMan = (TGeoManager*) fgeo->Get("FAIRGeom");
+  if (NULL == geoMan) {
+    cout << "<E> FAIRGeom not found in geoFile" << endl;
+    return;
+  }
+
+  // -----   Reconstruction run   -------------------------------------------
+  FairRunAna* run = new FairRunAna();
+  run->SetInputFile(InputFile.Data());
+  //run->AddFriend(InputFile.Data());
+  // run->SetOutputFile(OutputFile);
+  //run->SetSink( new FairRootFileSink( OutputFile.Data() ) );
+  run->SetUserOutputFileName(OutputFile.Data());
+  run->SetSink(new FairRootFileSink(run->GetUserOutputFileName()));
+  CbmTofEventClusterizer* tofClust =
+    new CbmTofEventClusterizer("TOF Event Clusterizer", iVerbose, bOut);
+
+  tofClust->SetCalMode(calMode);
+  tofClust->SetCalSel(calSel);
+  tofClust->SetCaldXdYMax(30.);  // geometrical matching window in cm
+  tofClust->SetCalCluMulMax(
+    3.);  // Max Counter Cluster Multiplicity for filling calib histos
+  tofClust->SetCalRpc(calSm);   // select detector for calibration update
+  tofClust->SetTRefId(RefSel);  // reference trigger for offset calculation
+  tofClust->SetTotMax(20.);     // Tot upper limit for walk corection
+  tofClust->SetTotMin(
+    0.01);  //(12000.);  // Tot lower limit for walk correction
+  tofClust->SetTotPreRange(
+    5.);  // effective lower Tot limit  in ns from peak position
+  tofClust->SetTotMean(5.);       // Tot calibration target value in ns
+  tofClust->SetMaxTimeDist(0.5);  // default cluster range in ns
+  //tofClust->SetMaxTimeDist(0.);       //Deb// default cluster range in ns
+  tofClust->SetDelTofMax(
+    4.);  // acceptance range for cluster distance in ns (!)
+  tofClust->SetSel2MulMax(3);  // limit Multiplicity in 2nd selector
+  tofClust->SetChannelDeadtime(dDeadtime);  // artificial deadtime in ns
+  tofClust->SetEnableAvWalk(kFALSE);        // for low statistics Cosmic
+  tofClust->SetEnableMatchPosScaling(
+    kFALSE);  // turn off projection to nominal target
+  tofClust->SetYFitMin(1.E4);
+  // tofClust->SetTimePeriod(25600.);       // ignore coarse time
+  // tofClust->SetCorMode(iBugCor);         // correct missing hits
+  tofClust->SetIdMode(0);  // 1 = calibrate on module level
+  //   tofClust->SetDeadStrips(15,23);   // declare dead strip for T0M3,Rpc0,Strip 23
+  //   tofClust->SetDeadStrips(27,30);   // declare dead strip for T6M0,Rpc1,Strip 30
+  //   tofClust->SetDeadStrips(28,19);   // declare dead strip for T7M0,Rpc0,Strip 19
+  //tofClust->SetDeadStrips(25,16);   // declare non-existant diamond strip (#5) dead
+
+  Int_t calSelRead = calSel;
+  if (calSel < 0) calSelRead = 0;
+  TString cFname = Form("%s_set%09d_%02d_%01dtofClust.hst.root",
+                        cFileId.Data(),
+                        iCalSet,
+                        calMode,
+                        calSelRead);
+  if (cCalId != "XXX")
+    cFname = Form("%s_set%09d_%02d_%01dtofClust.hst.root",
+                  cCalId.Data(),
+                  iCalSet,
+                  calMode,
+                  calSelRead);
+  tofClust->SetCalParFileName(cFname);
+  TString cOutFname =
+    Form("tofClust_%s_set%09d.hst.root", cFileId.Data(), iCalSet);
+  tofClust->SetOutHstFileName(cOutFname);
+
+  TString cAnaFile =
+    Form("%s_%09d%03d_tofAna.hst.root", cFileId.Data(), iCalSet, iSel2);
+
+  switch (calMode) {
+    case -1:                      // initial check of raw data
+      tofClust->SetTotMax(256.);  // range in bin number
+      tofClust->SetTotPreRange(256.);
+      //tofClust->SetTotMin(1.);
+      tofClust->SetTRefDifMax(26000.);  // in ns
+      tofClust->PosYMaxScal(10000.);    // in % of length
+      tofClust->SetMaxTimeDist(0.);     // no cluster building
+      //tofClust->SetTimePeriod(25600.);       // inspect coarse time
+      break;
+    case 0:                       // initial calibration
+      tofClust->SetTotMax(256.);  // range in bin number
+      tofClust->SetTotPreRange(256.);
+      //tofClust->SetTotMin(1.);
+      tofClust->SetTRefDifMax(1000.);  // in ns
+      tofClust->PosYMaxScal(100.);     // in % of length
+      tofClust->SetMaxTimeDist(0.);    // no cluster building
+      break;
+    case 1:                       // save offsets, update walks, for diamonds
+      tofClust->SetTotMax(256.);  // range in bin number
+      tofClust->SetTotPreRange(256.);
+      tofClust->SetTRefDifMax(6.25);  // in ns
+      //tofClust->SetTimePeriod(6.25);       // inspect coarse time
+      tofClust->PosYMaxScal(10.);  // in % of length
+      break;
+    case 11:
+      tofClust->SetTRefDifMax(5.);  // in ns
+      tofClust->PosYMaxScal(3.0);   // in % of length
+      break;
+    case 21:
+      tofClust->SetTRefDifMax(2.5);  // in ns
+      tofClust->PosYMaxScal(2.0);    // in % of length
+      break;
+    case 31:
+      tofClust->SetTRefDifMax(2.);  // in ns
+      tofClust->PosYMaxScal(1.5);   // in % of length
+      break;
+    case 41:
+      tofClust->SetTRefDifMax(1.5);  // in ns
+      tofClust->PosYMaxScal(0.8);    // in % of length
+      break;
+    case 51:
+      tofClust->SetTRefDifMax(1.2);  // in ns
+      tofClust->PosYMaxScal(0.7);    // in % of length
+      break;
+    case 61:
+      tofClust->SetTRefDifMax(1.0);  // in ns
+      tofClust->PosYMaxScal(0.75);   // in % of length
+      break;
+    case 71:
+      tofClust->SetTRefDifMax(0.8);  // in ns
+      tofClust->PosYMaxScal(0.6);    // in % of length
+      break;
+
+    case 2:                           // time difference calibration
+      tofClust->SetTRefDifMax(300.);  // in ns
+      tofClust->PosYMaxScal(1000.);   //in % of length
+      break;
+
+    case 3:                           // time offsets
+      tofClust->SetTRefDifMax(200.);  // in ns
+      tofClust->PosYMaxScal(10.);     //in % of length
+      tofClust->SetMaxTimeDist(0.);   // no cluster building
+      break;
+    case 12:
+    case 13:
+      tofClust->SetTRefDifMax(100.);  // in ns
+      tofClust->PosYMaxScal(10.);     //in % of length
+      break;
+    case 22:
+    case 23:
+      tofClust->SetTRefDifMax(40.);  // in ns
+      tofClust->PosYMaxScal(5.);     //in % of length
+      break;
+    case 32:
+    case 33:
+      tofClust->SetTRefDifMax(20.);  // in ns
+      tofClust->PosYMaxScal(4.);     //in % of length
+      break;
+    case 42:
+    case 43:
+      tofClust->SetTRefDifMax(10.);  // in ns
+      tofClust->PosYMaxScal(2.);     //in % of length
+      break;
+    case 52:
+    case 53:
+      tofClust->SetTRefDifMax(7.);  // in ns
+      tofClust->PosYMaxScal(1.5);   //in % of length
+      break;
+    case 62:
+    case 63:
+      tofClust->SetTRefDifMax(4.);  // in ns
+      tofClust->PosYMaxScal(1.);    //in % of length
+      break;
+    case 72:
+    case 73:
+      tofClust->SetTRefDifMax(2.8);  // in ns
+      tofClust->PosYMaxScal(0.9);    //in % of length
+      break;
+    case 82:
+    case 83:
+      tofClust->SetTRefDifMax(1.8);  // in ns
+      tofClust->PosYMaxScal(0.8);    //in % of length
+      break;
+    case 92:
+    case 93:
+      tofClust->SetTRefDifMax(1.0);  // in ns
+      tofClust->PosYMaxScal(0.75);   //in % of length
+      break;
+
+    case 4:  // velocity dependence (DelTOF)
+    case 14:
+      tofClust->SetTRefDifMax(25.);  // in ns
+      tofClust->PosYMaxScal(2.0);    //in % of length
+      break;
+    case 24:
+      tofClust->SetTRefDifMax(5.);  // in ns
+      tofClust->PosYMaxScal(1.5);   //in % of length
+      break;
+    case 34:
+      tofClust->SetTRefDifMax(5.);  // in ns
+      tofClust->PosYMaxScal(1.2);   //in % of length
+      break;
+    case 44:
+      tofClust->SetTRefDifMax(3.);  // in ns
+      tofClust->PosYMaxScal(1.0);   //in % of length
+      break;
+    case 54:
+      tofClust->SetTRefDifMax(2.5);  // in ns
+      tofClust->PosYMaxScal(0.9);    //in % of length
+      break;
+    case 64:
+      tofClust->SetTRefDifMax(2.5);  // in ns
+      tofClust->PosYMaxScal(0.8);    //in % of length
+      break;
+    case 74:
+      tofClust->SetTRefDifMax(1.5);  // in ns
+      tofClust->PosYMaxScal(0.7);    //in % of length
+      break;
+    default:
+      cout << "<E> Calib mode not implemented! stop execution of script"
+           << endl;
+      return;
+  }
+
+  run->AddTask(tofClust);
+
+  Int_t iBRef    = iCalSet % 1000;
+  Int_t iSet     = (iCalSet - iBRef) / 1000;
+  Int_t iRSel    = 0;
+  Int_t iRSelTyp = 0;
+  Int_t iRSelSm  = 0;
+  Int_t iRSelRpc = 0;
+  iRSel          = iBRef;  // use diamond
+  if (iSel2 == 0) {
+    // iSel2=iBRef;
+  } else {
+    if (iSel2 < 0) iSel2 = -iSel2;
+  }
+
+  Int_t iRSelin = iRSel;
+  iRSelRpc      = iRSel % 10;
+  iRSelTyp      = (iRSel - iRSelRpc) / 10;
+  iRSelSm       = iRSelTyp % 10;
+  iRSelTyp      = (iRSelTyp - iRSelSm) / 10;
+
+  tofClust->SetBeamRefId(iRSelTyp);  // define Beam reference counter
+  tofClust->SetBeamRefSm(iRSelSm);
+  tofClust->SetBeamRefDet(iRSelRpc);
+  tofClust->SetBeamAddRefMul(-1);
+  tofClust->SetBeamRefMulMax(3);
+
+  Int_t iSel2in  = iSel2;
+  Int_t iSel2Rpc = iSel2 % 10;
+  iSel2          = (iSel2 - iSel2Rpc) / 10;
+  Int_t iSel2Sm  = iSel2 % 10;
+  iSel2          = (iSel2 - iSel2Sm) / 10;
+  if (iSel2 > 0) {
+    tofClust->SetSel2Id(iSel2);
+    tofClust->SetSel2Sm(iSel2Sm);
+    tofClust->SetSel2Rpc(iSel2Rpc);
+  }
+
+  Int_t iRef    = iSet % 1000;
+  Int_t iDut    = (iSet - iRef) / 1000;
+  Int_t iDutRpc = iDut % 10;
+  iDut          = (iDut - iDutRpc) / 10;
+  Int_t iDutSm  = iDut % 10;
+  iDut          = (iDut - iDutSm) / 10;
+
+  tofClust->SetDutId(iDut);
+  tofClust->SetDutSm(iDutSm);
+  tofClust->SetDutRpc(iDutRpc);
+
+  Int_t iRefRpc = iRef % 10;
+  iRef          = (iRef - iRefRpc) / 10;
+  Int_t iRefSm  = iRef % 10;
+  iRef          = (iRef - iRefSm) / 10;
+
+  tofClust->SetSelId(iRef);
+  tofClust->SetSelSm(iRefSm);
+  tofClust->SetSelRpc(iRefRpc);
+
+  cout << "Run with iRSel = " << iRSel << ", iSel2 = " << iSel2in << endl;
+
+
+  // -----  Parameter database   --------------------------------------------
+
+  FairRuntimeDb* rtdb       = run->GetRuntimeDb();
+  Bool_t kParameterMerged   = kTRUE;
+  FairParRootFileIo* parIo2 = new FairParRootFileIo(kParameterMerged);
+  parIo2->open(ParFile.Data(), "UPDATE");
+  parIo2->print();
+  rtdb->setFirstInput(parIo2);
+
+  FairParAsciiFileIo* parIo1 = new FairParAsciiFileIo();
+  parIo1->open(parFileList, "in");
+  parIo1->print();
+  rtdb->setSecondInput(parIo1);
+  rtdb->print();
+  rtdb->printParamContexts();
+
+  //  FairParRootFileIo* parInput1 = new FairParRootFileIo();
+  //  parInput1->open(ParFile.Data());
+  //  rtdb->setFirstInput(parInput1);
+
+  // -----   Intialise and run   --------------------------------------------
+  run->Init();
+  cout << "Starting run" << endl;
+  run->Run(0, nEvents);
+  //tofClust->Finish();
+  // ------------------------------------------------------------------------
+  // default display
+  /*
+  TString Display_Status = "pl_over_Mat04D4best.C";
+  TString Display_Funct = "pl_over_Mat04D4best()";  
+  gROOT->LoadMacro(Display_Status);
+  */
+
+  gROOT->LoadMacro("save_hst.C");
+
+  gROOT->LoadMacro("fit_ybox.h");
+  gROOT->LoadMacro("pl_all_CluMul.C");
+  gROOT->LoadMacro("pl_all_CluRate.C");
+  gROOT->LoadMacro("pl_all_CluPosEvol.C");
+  gROOT->LoadMacro("pl_all_CluTimeEvol.C");
+  gROOT->LoadMacro("pl_over_cluSel.C");
+  gROOT->LoadMacro("pl_over_clu.C");
+  gROOT->LoadMacro("pl_over_Walk2.C");
+  gROOT->LoadMacro("pl_all_dTSel.C");
+  gROOT->LoadMacro("pl_over_MatD4sel.C");
+  gROOT->LoadMacro("pl_all_Sel2D.C");
+  gROOT->LoadMacro("pl_all_2D.C");
+
+  TString FSave = Form("save_hst(\"CluStatus%d_%d_Cal_%s.hst.root\")",
+                       iCalSet,
+                       iSel2in,
+                       cCalId.Data());
+  gInterpreter->ProcessLine(FSave.Data());
+
+  if (iPlot) {
+
+    switch (iSet) {
+      default:
+        for (Int_t iOpt = 0; iOpt < 7; iOpt++) {
+          for (Int_t iSel = 0; iSel < 2; iSel++) {
+            gInterpreter->ProcessLine(Form("pl_all_Sel2D(%d,%d)", iOpt, iSel));
+          }
+        }
+
+        for (Int_t iOpt = 0; iOpt < 12; iOpt++) {
+          gInterpreter->ProcessLine(Form("pl_all_2D(%d)", iOpt));
+        }
+        /*
+	gInterpreter->ProcessLine("pl_over_clu(0,0,0)");
+	gInterpreter->ProcessLine("pl_over_clu(0,0,1)");
+	gInterpreter->ProcessLine("pl_over_clu(0,0,2)");
+	gInterpreter->ProcessLine("pl_over_clu(0,0,3)");
+	gInterpreter->ProcessLine("pl_over_clu(0,0,4)");
+	gInterpreter->ProcessLine("pl_over_clu(0,1,0)");
+	gInterpreter->ProcessLine("pl_over_clu(0,1,1)");
+	gInterpreter->ProcessLine("pl_over_clu(0,1,2)");
+	gInterpreter->ProcessLine("pl_over_clu(0,1,3)");
+	gInterpreter->ProcessLine("pl_over_clu(0,1,4)");
+	gInterpreter->ProcessLine("pl_over_clu(0,2,0)");
+	gInterpreter->ProcessLine("pl_over_clu(0,2,1)");
+	gInterpreter->ProcessLine("pl_over_clu(0,2,2)");
+	gInterpreter->ProcessLine("pl_over_clu(0,2,3)");
+	gInterpreter->ProcessLine("pl_over_clu(0,2,4)");
+	gInterpreter->ProcessLine("pl_over_clu(0,3,0)");
+	gInterpreter->ProcessLine("pl_over_clu(0,3,1)");
+	gInterpreter->ProcessLine("pl_over_clu(0,3,2)");
+	gInterpreter->ProcessLine("pl_over_clu(0,3,3)");
+	gInterpreter->ProcessLine("pl_over_clu(0,3,4)");
+	gInterpreter->ProcessLine("pl_over_clu(0,4,0)");
+	gInterpreter->ProcessLine("pl_over_clu(0,4,1)");
+	gInterpreter->ProcessLine("pl_over_clu(0,4,2)");
+	gInterpreter->ProcessLine("pl_over_clu(0,4,3)");
+	gInterpreter->ProcessLine("pl_over_clu(0,4,4)");
+	
+	gInterpreter->ProcessLine("pl_over_clu(5,0,0)");
+	gInterpreter->ProcessLine("pl_over_cluSel(0,5,0,0)");
+	gInterpreter->ProcessLine("pl_over_cluSel(1,5,0,0)");
+	
+	for(Int_t iSm=0; iSm<3; iSm++)
+	for (Int_t iRpc=0; iRpc<5; iRpc++)
+	for (Int_t iSel=0; iSel<2; iSel++){
+	gInterpreter->ProcessLine(Form("pl_over_cluSel(%d,0,%d,%d)",iSel,iSm,iRpc));
+	gInterpreter->ProcessLine(Form("pl_over_Walk2(%d,0,%d,%d)",iSel,iSm,iRpc));
+	}
+      */
+        gInterpreter->ProcessLine("pl_all_CluMul()");
+        gInterpreter->ProcessLine("pl_all_CluRate()");
+        gInterpreter->ProcessLine("pl_all_CluRate(5,1)");
+        gInterpreter->ProcessLine("pl_all_CluPosEvol()");
+        gInterpreter->ProcessLine("pl_all_CluTimeEvol()");
+        gInterpreter->ProcessLine("pl_all_dTSel()");
+
+        //  gInterpreter->ProcessLine("pl_over_MatD4sel()");
+        //  gInterpreter->ProcessLine(Display_Funct.Data());
+        break;
+        ;
+    }
+  }
+}
diff --git a/macro/beamtime/hd2020/ana_trks.C b/macro/beamtime/hd2020/ana_trks.C
new file mode 100644
index 0000000000000000000000000000000000000000..12163e4512e03761063b8df4400addfcec47f497
--- /dev/null
+++ b/macro/beamtime/hd2020/ana_trks.C
@@ -0,0 +1,636 @@
+void ana_trks(Int_t nEvents        = 10000,
+              Int_t iSel           = 910920,
+              Int_t iGenCor        = 1,
+              TString cFileId      = "HD_cosmic_2020-02-08_18:50:47.50.3.0.0",
+              TString cSet         = "900920910",
+              Int_t iSel2          = 911,
+              Int_t iTrackingSetup = 1,
+              Double_t dScalFac    = 1.,
+              Double_t dChi2Lim2   = 50.,
+              Double_t dDeadtime   = 50,
+              TString cCalId       = "",
+              Int_t iAnaCor        = 1,
+              Bool_t bUseSigCalib  = kFALSE,
+	          Int_t iCalSet        = 900920910,
+              Int_t iCalOpt        = 1,
+              Int_t iMc            = 0) {
+  Int_t iVerbose = 1;
+  if (cCalId == "") cCalId = cFileId;
+  // Specify log level (INFO, DEBUG, DEBUG1, ...)
+  //TString logLevel = "FATAL";
+  //TString logLevel = "ERROR";
+  TString logLevel = "INFO";
+  //TString logLevel = "DEBUG";
+  //TString logLevel = "DEBUG1";
+  //TString logLevel = "DEBUG2";
+  //TString logLevel = "DEBUG3";
+  TString workDir = gSystem->Getenv("VMCWORKDIR");
+  //TString workDir          = "../../..";
+  TString paramDir      = workDir + "/macro/beamtime/hd2020";
+  TString ParFile       = paramDir + "/data/" + cFileId.Data() + ".params.root";
+  TString InputFile     = paramDir + "/data/" + cFileId.Data() + ".root";
+  TString InputDigiFile = paramDir + "/data/digidev_" + cFileId.Data()
+                          + Form("_%s_%02.0f_Cal", cSet.Data(), dDeadtime)
+                          + cCalId + ".out.root";
+  if (iMc == 1) {
+    InputFile     = paramDir + "/data/" + cFileId.Data() + ".raw.root";
+    InputDigiFile = paramDir + "/data/" + cFileId.Data() + ".rec.root";
+  }
+  TString OutputFile = paramDir + "/data/hits_" + cFileId.Data()
+                       + Form("_%s_%06d_%03d", cSet.Data(), iSel, iSel2)
+                       + ".out.root";
+  TString cHstFile =
+    paramDir
+    + Form(
+        "/hst/%s_%03.0f_%s_%06d_%03d_%03.1f_%03.1f_trk%03d_Cal%s_Ana.hst.root",
+        cFileId.Data(),
+        dDeadtime,
+        cSet.Data(),
+        iSel,
+        iSel2,
+        dScalFac,
+        dChi2Lim2,
+        iTrackingSetup,
+        cCalId.Data());
+  TString cTrkFile = Form("%s_tofFindTracks.hst.root", cCalId.Data());
+  TString cAnaFile = Form("%s_TrkAnaTestBeam.hst.root", cFileId.Data());
+
+  cout << " InputDigiFile = " << InputDigiFile << endl;
+
+  TString shcmd = "rm -v " + ParFile;
+  gSystem->Exec(shcmd.Data());
+
+  TList* parFileList = new TList();
+
+  //TString TofGeo="v18m_cosmicHD";  //Buc2015 setup, r0078
+  //TString TofGeo = "v19a_cosmicHD";  //Buc2018 sandwiched, r0079
+  TString TofGeo = "v20b_cosmicHD";  //Buc2020 sandwiched
+
+  cout << "Geometry version " << TofGeo << endl;
+
+  TObjString* tofDigiFile = new TObjString(
+    workDir + "/parameters/tof/tof_" + TofGeo + ".digi.par");  // TOF digi file
+  parFileList->Add(tofDigiFile);
+
+  // TObjString tofDigiBdfFile =  paramDir + "/tof.digibdf.par";
+  // TObjString tofDigiBdfFile =  paramDir + "/tof." + FPar + "digibdf.par";
+  TObjString* tofDigiBdfFile =
+    new TObjString(workDir + "/parameters/tof/tof_" + TofGeo + ".digibdf.par");
+  parFileList->Add(tofDigiBdfFile);
+
+  TString geoDir  = gSystem->Getenv("VMCWORKDIR");
+  TString geoFile = geoDir + "/geometry/tof/geofile_tof_" + TofGeo + ".root";
+  TFile* fgeo     = new TFile(geoFile);
+  TGeoManager* geoMan = (TGeoManager*) fgeo->Get("FAIRGeom");
+  if (NULL == geoMan) {
+    cout << "<E> FAIRGeom not found in geoFile" << endl;
+    return;
+  }
+
+  // -----   Reconstruction run   -------------------------------------------
+  FairRunAna* run = new FairRunAna();
+  cout << "InputFile:     " << InputFile.Data() << endl;
+  cout << "InputDigiFile: " << InputDigiFile.Data() << endl;
+
+  //run->SetInputFile(InputFile.Data());
+  //run->AddFriend(InputDigiFile.Data());
+  run->SetInputFile(InputDigiFile.Data());
+  //run->AddFriend(InputFile.Data());
+  //run->SetOutputFile(OutputFile);
+  run->SetUserOutputFileName(OutputFile.Data());
+  run->SetSink(new FairRootFileSink(run->GetUserOutputFileName()));
+
+  FairLogger::GetLogger()->SetLogScreenLevel(logLevel.Data());
+  FairLogger::GetLogger()->SetLogVerbosityLevel("MEDIUM");
+
+  // -----   Local selection variables  -------------------------------------------
+
+  Int_t iRef    = iSel % 1000;
+  Int_t iDut    = (iSel - iRef) / 1000;
+  Int_t iDutRpc = iDut % 10;
+  iDut          = (iDut - iDutRpc) / 10;
+  Int_t iDutSm  = iDut % 10;
+  iDut          = (iDut - iDutSm) / 10;
+  Int_t iRefRpc = iRef % 10;
+  iRef          = (iRef - iRefRpc) / 10;
+  Int_t iRefSm  = iRef % 10;
+  iRef          = (iRef - iRefSm) / 10;
+
+  Int_t iSel2in  = iSel2;
+  Int_t iSel2Rpc = iSel2 % 10;
+  iSel2          = (iSel2 - iSel2Rpc) / 10;
+  Int_t iSel2Sm  = iSel2 % 10;
+  iSel2          = (iSel2 - iSel2Sm) / 10;
+
+
+  Int_t calMode = 93;
+  Int_t calSel  = 1;
+  Bool_t bOut   = kFALSE;
+
+  CbmTofEventClusterizer* tofClust =
+    new CbmTofEventClusterizer("TOF Event Clusterizer", iVerbose, bOut);
+  Int_t calSelRead = calSel;
+  if (calSel < 0) calSelRead = 0;
+  TString cFname = Form("%s_set%09d_%02d_%01dtofClust.hst.root",
+                        cFileId.Data(),
+                        iCalSet,
+                        calMode,
+                        calSelRead);
+  if (cCalId != "XXX")
+    cFname = Form("%s_set%09d_%02d_%01dtofClust.hst.root",
+                  cCalId.Data(),
+                  iCalSet,
+                  calMode,
+                  calSelRead);
+  tofClust->SetCalParFileName(cFname);
+  TString cOutFname =
+    Form("tofClust_%s_set%09d.hst.root", cFileId.Data(), iCalSet);
+  tofClust->SetOutHstFileName(cOutFname);
+
+  // =========================================================================
+  // ===                       Tracking                                    ===
+  // =========================================================================
+
+  CbmTofTrackFinder* tofTrackFinder = new CbmTofTrackFinderNN();
+  tofTrackFinder->SetMaxTofTimeDifference(0.2);  // in ns/cm
+  tofTrackFinder->SetTxLIM(1.6);                 // max slope dx/dz
+  tofTrackFinder->SetTyLIM(1.6);                // max dev from mean slope dy/dz
+  tofTrackFinder->SetTyMean(0.);                // mean slope dy/dz
+  tofTrackFinder->SetMaxTofTimeDifference(1.);  // in ns/cm
+
+  CbmTofTrackFitter* tofTrackFitter = new CbmTofTrackFitterKF(0, 211);
+  TFitter* MyFit                    = new TFitter(1);  // initialize Minuit
+  tofTrackFinder->SetFitter(tofTrackFitter);
+  CbmTofFindTracks* tofFindTracks = new CbmTofFindTracks("TOF Track Finder");
+  tofFindTracks->UseFinder(tofTrackFinder);
+  tofFindTracks->UseFitter(tofTrackFitter);
+  tofFindTracks->SetCalOpt(
+    iCalOpt);  // 1 - update offsets, 2 - update walk, 0 - bypass
+  tofFindTracks->SetCorMode(iGenCor);  // valid options: 0,1,2,3,4,5,6, 10 - 19
+  //tofFindTracks->SetTtTarg(0.057);            // target value for inverse velocity, > 0.033 ns/cm!
+  tofFindTracks->SetTtTarg(
+    0.035);  // target value for inverse velocity, > 0.033 ns/cm!
+  tofFindTracks->SetCalParFileName(
+    cTrkFile);  // Tracker parameter value file name
+  //   tofFindTracks->SetBeamCounter(5,0,0);         // default beam counter
+  tofFindTracks->SetStationMaxHMul(
+    30);  // Max Hit Multiplicity in any used station
+
+  tofFindTracks->SetT0MAX(dScalFac);  // in ns
+  tofFindTracks->SetSIGT(0.08);       // default in ns
+  tofFindTracks->SetSIGX(0.3);        // default in cm
+  tofFindTracks->SetSIGY(0.6);        // default in cm
+  tofFindTracks->SetSIGZ(0.05);       // default in cm
+  tofFindTracks->SetUseSigCalib(
+    bUseSigCalib);  // ignore resolutions in CalPar file
+  tofFindTracks->SetRefVelMean(
+    30.);  // set mean velocity value for acceptable tracklets
+  //   tofFindTracks->SetRefDVel(5.);         // set velocity window width
+  tofTrackFinder->SetSIGLIM(dChi2Lim2
+                            * 2.);  // matching window in multiples of chi2
+  tofTrackFinder->SetChiMaxAccept(dChi2Lim2);  // max tracklet chi2
+
+  Int_t iMinNofHits   = -1;
+  Int_t iNStations    = 0;
+  Int_t iNReqStations = 3;
+
+  switch (iTrackingSetup) {
+    case 0:  // bypass mode
+      iMinNofHits = -1;
+      iNStations  = 1;
+      tofFindTracks->SetStation(0, 5, 0, 0);  // Diamond
+      break;
+
+    case 1:  // for calibration mode of full setup
+      iMinNofHits   = 3;
+      iNStations    = 8;
+      iNReqStations = 4;
+      tofFindTracks->SetStation(0, 9, 0, 0);
+      tofFindTracks->SetStation(1, 9, 2, 1);
+      tofFindTracks->SetStation(2, 9, 0, 1);
+      tofFindTracks->SetStation(3, 9, 2, 0);
+      tofFindTracks->SetStation(4, 9, 1, 0);
+      tofFindTracks->SetStation(5, 9, 1, 1);
+      tofFindTracks->SetStation(6, 6, 0, 0);
+      tofFindTracks->SetStation(7, 6, 0, 1);
+      //tofFindTracks->SetStation(8, 6, 1, 0);
+      //tofFindTracks->SetStation(9, 6, 1, 1);
+      break;
+      
+    case 11:  // for calibration mode of full setup
+      iMinNofHits   = 3;
+      iNStations    = 6;
+      iNReqStations = 4;
+      tofFindTracks->SetStation(0, 9, 2, 1);
+      tofFindTracks->SetStation(1, 9, 0, 0);
+      tofFindTracks->SetStation(2, 9, 1, 1);
+      tofFindTracks->SetStation(3, 9, 0, 1);
+      tofFindTracks->SetStation(4, 6, 0, 0);
+      tofFindTracks->SetStation(5, 6, 0, 1);
+      break;      
+
+    case 2:
+      iMinNofHits   = 6;
+      iNStations    = 7;
+      iNReqStations = 7;
+      tofFindTracks->SetStation(0, 9, 0, 0);
+      tofFindTracks->SetStation(1, 9, 2, 1);
+      tofFindTracks->SetStation(2, 9, 0, 1);
+      tofFindTracks->SetStation(3, 9, 2, 0);
+      tofFindTracks->SetStation(4, 9, 1, 0);
+      tofFindTracks->SetStation(5, 9, 1, 1);
+      tofFindTracks->SetStation(6, iDut, iDutSm, iDutRpc);
+      //     tofFindTracks->SetStation(4, 9, 0, 0);
+      break;
+
+    case 21:
+      iMinNofHits   = 4;
+      iNStations    = 5;
+      iNReqStations = 5;
+      tofFindTracks->SetStation(0, 9, 0, 0);
+      tofFindTracks->SetStation(1, 9, 2, 1);
+      tofFindTracks->SetStation(2, 9, 1, 1);
+      tofFindTracks->SetStation(3, 9, 0, 1);
+      //tofFindTracks->SetStation(4, 9, 0, 0);
+      //tofFindTracks->SetStation(4, 9, 0, 1);
+      tofFindTracks->SetStation(4, iDut, iDutSm, iDutRpc);
+      break;
+
+    case 22:
+      iMinNofHits   = 5;
+      iNStations    = 7;
+      iNReqStations = 7;
+      tofFindTracks->SetStation(0, 9, 1, 0);
+      tofFindTracks->SetStation(1, 9, 2, 0);
+      tofFindTracks->SetStation(2, 9, 1, 1);
+      tofFindTracks->SetStation(3, 9, 2, 1);
+      //tofFindTracks->SetStation(4, 9, 0, 0);
+      tofFindTracks->SetStation(4, 9, 0, 1);
+      tofFindTracks->SetStation(5, 1, 0, 2);
+      tofFindTracks->SetStation(6, iDut, iDutSm, iDutRpc);
+      //     tofFindTracks->SetStation(4, 9, 0, 0);
+      break;
+
+    case 3:
+      iMinNofHits   = 3;
+      iNStations    = 4;
+      iNReqStations = 4;
+      tofFindTracks->SetStation(0, 9, 0, 0);
+      tofFindTracks->SetStation(1, 9, 2, 1);
+      tofFindTracks->SetStation(2, 9, 0, 1);
+      tofFindTracks->SetStation(3, iDut, iDutSm, iDutRpc);
+      break;
+
+    case 4:
+      iMinNofHits   = 6;
+      iNStations    = 7;
+      iNReqStations = 7;
+      tofFindTracks->SetStation(0, 9, 1, 0);
+      tofFindTracks->SetStation(1, 9, 2, 0);
+      tofFindTracks->SetStation(2, 9, 1, 1);
+      tofFindTracks->SetStation(3, 9, 2, 1);
+      tofFindTracks->SetStation(4, 9, 0, 1);
+      tofFindTracks->SetStation(5, 9, 0, 0);
+      tofFindTracks->SetStation(6, iDut, iDutSm, iDutRpc);
+      break;
+
+    case 5:  // for calibration of 2-stack and add-on counters (STAR2, CERN)
+      iMinNofHits   = 7;
+      iNStations    = 8;
+      iNReqStations = 8;
+      tofFindTracks->SetStation(0, 9, 1, 0);
+      tofFindTracks->SetStation(1, 9, 2, 0);
+      tofFindTracks->SetStation(2, 9, 1, 1);
+      tofFindTracks->SetStation(3, 9, 2, 1);
+      tofFindTracks->SetStation(4, 9, 0, 0);
+      tofFindTracks->SetStation(5, 9, 0, 1);
+      tofFindTracks->SetStation(6, 6, 0, 0);
+      tofFindTracks->SetStation(7, 6, 0, 1);
+      break;
+
+    case 6:
+      iMinNofHits   = 5;
+      iNStations    = 6;
+      iNReqStations = 6;
+      tofFindTracks->SetStation(0, 9, 0, 0);
+      tofFindTracks->SetStation(1, 9, 2, 0);
+      tofFindTracks->SetStation(2, 9, 0, 1);
+      tofFindTracks->SetStation(3, 9, 2, 1);
+      tofFindTracks->SetStation(4, 9, 1, 0);
+      tofFindTracks->SetStation(5, 9, 1, 1);
+      break;
+
+    case 7:  // for calibration of 2-stack and add-on counters (BUC)
+      iMinNofHits   = 3;
+      iNStations    = 5;
+      iNReqStations = 4;
+      tofFindTracks->SetStation(0, 5, 0, 0);
+      tofFindTracks->SetStation(1, 0, 4, 3);
+      tofFindTracks->SetStation(2, 0, 3, 3);
+      tofFindTracks->SetStation(3, 6, 0, 0);
+      tofFindTracks->SetStation(4, 6, 0, 1);
+      break;
+
+    case 8:  // evaluation of add-on counters (BUC)
+      iMinNofHits   = 3;
+      iNStations    = 4;
+      iNReqStations = 4;
+      tofFindTracks->SetStation(0, 5, 0, 0);
+      tofFindTracks->SetStation(1, 0, 4, 3);
+      tofFindTracks->SetStation(2, 0, 3, 3);
+      tofFindTracks->SetStation(3, iDut, iDutSm, iDutRpc);
+      break;
+
+    case 10:
+      iMinNofHits   = 3;
+      iNStations    = 4;
+      iNReqStations = 4;
+      tofFindTracks->SetStation(0, 5, 0, 0);
+      tofFindTracks->SetStation(1, 0, 1, 2);
+      tofFindTracks->SetStation(2, 0, 0, 2);
+      tofFindTracks->SetStation(3, 0, 2, 2);
+
+    default:
+      cout << "Tracking setup " << iTrackingSetup << " not implemented "
+           << endl;
+      return;
+      ;
+  }
+  tofFindTracks->SetMinNofHits(iMinNofHits);
+  tofFindTracks->SetNStations(iNStations);
+  tofFindTracks->SetNReqStations(iNReqStations);
+  tofFindTracks->PrintSetup();
+  run->AddTask(tofFindTracks);
+
+  // =========================================================================
+  // ===                       Analysis                                    ===
+  // =========================================================================
+  CbmTofAnaTestbeam* tofAnaTestbeam =
+    new CbmTofAnaTestbeam("TOF TestBeam Analysis", iVerbose);
+  tofAnaTestbeam->SetCorMode(iAnaCor);  // 1 - DTD4, 2 - X4, 3 - Y4, 4 - Texp
+  tofAnaTestbeam->SetHitDistMin(30.);   // initialization
+  //tofAnaTestbeam->SetEnableMatchPosScaling(kTRUE);
+  tofAnaTestbeam->SetEnableMatchPosScaling(kFALSE);
+  tofAnaTestbeam->SetSpillDuration(30000.);
+  //CbmTofAnaTestbeam defaults
+  //   tofAnaTestbeam->SetR0LimFit(20.);  // limit distance of fitted track to nominal vertex
+  tofAnaTestbeam->SetDXMean(0.);
+  tofAnaTestbeam->SetDYMean(0.);
+  tofAnaTestbeam->SetDTMean(0.);  // in ns
+  tofAnaTestbeam->SetDXWidth(0.5);
+  tofAnaTestbeam->SetDYWidth(1.0);
+  tofAnaTestbeam->SetDTWidth(0.1);  // in ns
+  tofAnaTestbeam->SetCalParFileName(cAnaFile);
+  Double_t dScalFacA = 0.9;  // dScalFac is used for tracking
+  tofAnaTestbeam->SetPosY4Sel(
+    0.5 * dScalFacA);  // Y Position selection in fraction of strip length
+  tofAnaTestbeam->SetDTDia(0.);    // Time difference to additional diamond
+  tofAnaTestbeam->SetMul0Max(20);  // Max Multiplicity in dut
+  tofAnaTestbeam->SetMul4Max(30);  // Max Multiplicity in Ref - RPC
+  tofAnaTestbeam->SetMulDMax(3);   // Max Multiplicity in Diamond / BeamRef
+  tofAnaTestbeam->SetTOffD4(14.);  // initialization
+  tofAnaTestbeam->SetDTD4MAX(
+    6.);  // initialization of Max time difference Ref - BRef
+
+  //tofAnaTestbeam->SetTShift(-28000.);// initialization
+  tofAnaTestbeam->SetPosYS2Sel(
+    0.55);  // Y Position selection in fraction of strip length
+  tofAnaTestbeam->SetChS2Sel(0.);     // Center of channel selection window
+  tofAnaTestbeam->SetDChS2Sel(100.);  // Width  of channel selection window
+  tofAnaTestbeam->SetSel2TOff(0.);    // Shift Sel2 time peak to 0
+  tofAnaTestbeam->SetChi2Lim(5.);     // initialization of Chi2 selection limit
+  tofAnaTestbeam->SetChi2Lim2(
+    3.);  // initialization of Chi2 selection limit for Mref-Sel2 pair
+  tofAnaTestbeam->SetDutDX(
+    15.);  // limit inspection of tracklets to selected region
+  tofAnaTestbeam->SetDutDY(
+    15.);  // limit inspection of tracklets to selected region
+  tofAnaTestbeam->SetSIGLIM(3.);  // max matching chi2
+  tofAnaTestbeam->SetSIGT(0.08);  // in ns
+  tofAnaTestbeam->SetSIGX(0.3);   // in cm
+  tofAnaTestbeam->SetSIGY(0.6);   // in cm
+
+  Int_t iRSel    = 921;
+  Int_t iRSelTyp = 9;
+  Int_t iRSelSm  = 2;
+  Int_t iRSelRpc = 1;
+  Int_t iRSelin  = iRSel;
+
+  tofAnaTestbeam->SetBeamRefSmType(iRSelTyp);  // common reaction reference
+  tofAnaTestbeam->SetBeamRefSmId(iRSelSm);
+  tofAnaTestbeam->SetBeamRefRpc(iRSelRpc);
+
+  if (iSel2 >= 0) {
+    tofAnaTestbeam->SetMrpcSel2(
+      iSel2);  // initialization of second selector Mrpc Type
+    tofAnaTestbeam->SetMrpcSel2Sm(
+      iSel2Sm);  // initialization of second selector Mrpc SmId
+    tofAnaTestbeam->SetMrpcSel2Rpc(
+      iSel2Rpc);  // initialization of second selector Mrpc RpcId
+  }
+
+  tofAnaTestbeam->SetDut(iDut);            // Device under test
+  tofAnaTestbeam->SetDutSm(iDutSm);        // Device under test
+  tofAnaTestbeam->SetDutRpc(iDutRpc);      // Device under test
+  tofAnaTestbeam->SetMrpcRef(iRef);        // Reference RPC
+  tofAnaTestbeam->SetMrpcRefSm(iRefSm);    // Reference RPC
+  tofAnaTestbeam->SetMrpcRefRpc(iRefRpc);  // Reference RPC
+
+  cout << "dispatch iSel = " << iSel << ", iSel2in = " << iSel2in
+       << ", iRSelin = " << iRSelin << ", iRSel = " << iRSel << endl;
+  if (1) {
+    switch (iSel) {
+
+      case 600921:
+      case 601921:
+      case 700921:
+      case 911921:
+      case 920921:
+        switch (iRSelin) {
+          case 920:
+            tofAnaTestbeam->SetTShift(0.6);  // Shift DTD4 to 0
+            tofAnaTestbeam->SetTOffD4(18.);  // Shift DTD4 to physical value
+            switch (iSel2in) {
+              case 910:
+                tofAnaTestbeam->SetSel2TOff(0.);  // Shift Sel2 time peak to 0
+                break;
+              default:;
+            }
+            break;
+
+          case 910:
+            tofAnaTestbeam->SetTShift(-1.);  // Shift DTD4 to 0
+            tofAnaTestbeam->SetTOffD4(18.);  // Shift DTD4 to physical value
+            switch (iSel2in) {
+              case 910:
+                tofAnaTestbeam->SetSel2TOff(0.);  // Shift Sel2 time peak to 0
+                break;
+              default:;
+            }
+            break;
+
+          default:;
+        }
+        break;
+
+      case 600910:
+      case 601910:
+      case 700910:
+      case 900910:
+      case 901910:
+        switch (iRSelin) {
+          case 920:
+            tofAnaTestbeam->SetTShift(0.4);  // Shift DTD4 to 0
+            tofAnaTestbeam->SetTOffD4(18.);  // Shift DTD4 to physical value
+
+            switch (iSel2in) {
+              case 920:
+                tofAnaTestbeam->SetSel2TOff(0.5);  // Shift Sel2 time peak to 0
+                break;
+              case 921:
+                tofAnaTestbeam->SetSel2TOff(0.41);  // Shift Sel2 time peak to 0
+                break;
+
+              default:;
+            }
+            break;
+          default:;
+        }
+        break;
+
+      case 101900:
+      case 600900:
+      case 601900:
+        switch (iRSelin) {
+          case 920:
+            tofAnaTestbeam->SetTShift(0.0);  // Shift DTD4 to 0
+            tofAnaTestbeam->SetTOffD4(11.);  // Shift DTD4 to physical value
+
+            switch (iSel2in) {
+              case 910:
+                tofAnaTestbeam->SetSel2TOff(0.);  // Shift Sel2 time peak to 0
+                break;
+              case 911:
+                tofAnaTestbeam->SetSel2TOff(0.);  // Shift Sel2 time peak to 0
+                break;
+
+              default:;
+            }
+            break;
+          default:;
+        }
+        break;
+
+      case 600043:
+      case 601043:
+        switch (iRSelin) {
+          case 500:
+            tofAnaTestbeam->SetTShift(5.3);  // Shift DTD4 to 0
+            tofAnaTestbeam->SetTOffD4(11.);  // Shift DTD4 to physical value
+
+            switch (iSel2in) {
+              case 33:
+                tofAnaTestbeam->SetSel2TOff(
+                  -0.55);  // Shift Sel2 time peak to 0
+                break;
+
+              default:;
+            }
+            break;
+          default:;
+        }
+        break;
+
+      default:
+        cout << "Better to define analysis setup! Running with default offset "
+                "parameter... "
+             << endl;
+        // return;
+    }  // end of different subsets
+
+    cout << " Initialize TSHIFT to " << tofAnaTestbeam->GetTShift() << endl;
+    run->AddTask(tofAnaTestbeam);
+  }
+
+  // -----  Parameter database   --------------------------------------------
+  FairRuntimeDb* rtdb       = run->GetRuntimeDb();
+  Bool_t kParameterMerged   = kTRUE;
+  FairParRootFileIo* parIo2 = new FairParRootFileIo(kParameterMerged);
+  parIo2->open(ParFile.Data(), "UPDATE");
+  parIo2->print();
+  rtdb->setFirstInput(parIo2);
+
+  FairParAsciiFileIo* parIo1 = new FairParAsciiFileIo();
+  parIo1->open(parFileList, "in");
+  parIo1->print();
+  rtdb->setSecondInput(parIo1);
+  rtdb->print();
+  rtdb->printParamContexts();
+
+  //  FairParRootFileIo* parInput1 = new FairParRootFileIo();
+  //  parInput1->open(ParFile.Data());
+  //  rtdb->setFirstInput(parInput1);
+
+  // -----   Intialise and run   --------------------------------------------
+  run->Init();
+  cout << "Starting run" << endl;
+  run->Run(0, nEvents);
+  //run->Run(nEvents-1, nEvents); //debugging single events for memory leak
+  // ------------------------------------------------------------------------
+  TString SaveToHstFile = "save_hst(\"" + cHstFile + "\")";
+  gROOT->LoadMacro("save_hst.C");
+  gInterpreter->ProcessLine(SaveToHstFile);
+
+  // default displays, plot results
+  TString Display_Status = "pl_over_Mat04D4best.C";
+  TString Display_Funct;
+  if (iGenCor < 0) {
+    Display_Funct = "pl_over_Mat04D4best(1)";
+  } else {
+    Display_Funct = "pl_over_Mat04D4best(0)";
+  }
+  gROOT->LoadMacro(Display_Status);
+
+  cout << "Exec " << Display_Funct.Data() << endl;
+  gInterpreter->ProcessLine(Display_Funct);
+
+  gROOT->LoadMacro("pl_over_MatD4sel.C");
+  gROOT->LoadMacro("pl_eff_XY.C");
+  gROOT->LoadMacro("pl_over_trk.C");
+  gROOT->LoadMacro("pl_calib_trk.C");
+  gROOT->LoadMacro("pl_XY_trk.C");
+  //gROOT->LoadMacro("pl_vert_trk.C");
+  gROOT->LoadMacro("pl_pull_trk.C");
+  gROOT->LoadMacro("pl_all_Track2D.C");
+  gROOT->LoadMacro("pl_TIS.C");
+  gROOT->LoadMacro("pl_TIR.C");
+  gROOT->LoadMacro("pl_Eff_XY.C");
+  gROOT->LoadMacro("pl_Eff_DTLH.C");
+  gROOT->LoadMacro("pl_Eff_TIS.C");
+  gROOT->LoadMacro("pl_Dut_Res.C");
+
+  gInterpreter->ProcessLine("pl_over_MatD4sel()");
+  gInterpreter->ProcessLine("pl_TIS()");
+  gInterpreter->ProcessLine("pl_TIR()");
+  //gInterpreter->ProcessLine("pl_eff_XY()");
+  gInterpreter->ProcessLine("pl_calib_trk()");
+
+  gInterpreter->ProcessLine("pl_all_Track2D(1)");
+  gInterpreter->ProcessLine("pl_all_Track2D(2)");
+  gInterpreter->ProcessLine("pl_all_Track2D(4)");
+  
+  TString over_trk = "pl_over_trk(" + (TString)(Form("%d", iNStations)) + ")";
+  gInterpreter->ProcessLine(over_trk);
+
+  TString XY_trk = "pl_XY_trk(" + (TString)(Form("%d", iNStations)) + ")";
+  gInterpreter->ProcessLine(XY_trk);
+
+  TString Pull0 = (TString)(Form("pl_pull_trk(%d,%d,1)", iNStations, 0));
+  gInterpreter->ProcessLine(Pull0);
+  TString Pull1 = (TString)(Form("pl_pull_trk(%d,%d,1)", iNStations, 1));
+  gInterpreter->ProcessLine(Pull1);
+  TString Pull3 = (TString)(Form("pl_pull_trk(%d,%d,1)", iNStations, 3));
+  gInterpreter->ProcessLine(Pull3);
+  TString Pull4 = (TString)(Form("pl_pull_trk(%d,%d,1)", iNStations, 4));
+  gInterpreter->ProcessLine(Pull4);
+}
diff --git a/macro/beamtime/hd2020/dis_digi.C b/macro/beamtime/hd2020/dis_digi.C
new file mode 100644
index 0000000000000000000000000000000000000000..c70f557b9c09dc59a251ba8b37b4af1df86b3c00
--- /dev/null
+++ b/macro/beamtime/hd2020/dis_digi.C
@@ -0,0 +1,716 @@
+void dis_digi(Int_t nEvents = 100, Int_t calMode=93, Int_t calSel=1, Int_t calSm=0, Int_t RefSel=1, 
+TString cFileId="68.50.7.1", Int_t iCalSet=10500, Bool_t bOut=0, Int_t iSel2=20, Double_t dDeadtime=50, 
+Int_t iGenCor=1, Int_t iTrackingSetup=1, Double_t dScalFac=5., Double_t dChi2Lim2=10., TString cCalId="XXX", Bool_t bUseSigCalib=kFALSE) 
+{
+  Int_t iVerbose = 1;
+  //Specify log level (INFO, DEBUG, DEBUG1, ...)
+  //TString logLevel = "FATAL";
+  //TString logLevel = "ERROR";
+  //TString logLevel = "INFO";
+  //TString logLevel = "DEBUG";
+  TString logLevel = "DEBUG1";
+  //TString logLevel = "DEBUG2";
+  //TString logLevel = "DEBUG3";
+  FairLogger::GetLogger();
+  gLogger->SetLogScreenLevel(logLevel);
+  gLogger->SetLogVerbosityLevel("MEDIUM");
+  
+  TString workDir    = gSystem->Getenv("VMCWORKDIR");
+  /*
+   TString workDir    = (TString)gInterpreter->ProcessLine(".! pwd");
+   cout << "workdir = "<< workDir.Data() << endl;
+   return;
+  */
+   TString paramDir   = workDir + "/macro/beamtime/hd2020/";
+   TString ParFile    = paramDir + "data/" + cFileId + ".params.root";
+   TString InputFile  = paramDir + "data/" + cFileId + ".root";
+   TString OutputFile = paramDir + "data/disdigi_"   + cFileId + Form("_%09d%03d",iCalSet,iSel2) + ".out.root";
+
+   TString cTrkFile=Form("%s_tofFindTracks.hst.root",cFileId.Data());
+
+   TList *parFileList = new TList();
+
+   TString FId=cFileId;
+   //TString TofGeo="v18o_cosmicHD";
+   TString TofGeo="v20b_cosmicHD"; // for Buc2020 sandwiched
+    cout << "Geometry version "<<TofGeo<<endl;
+
+   TObjString *tofDigiFile = new TObjString(workDir + "/parameters/tof/tof_" + TofGeo + ".digi.par"); // TOF digi file
+   parFileList->Add(tofDigiFile);   
+
+   //   TObjString tofDigiBdfFile = new TObjString( paramDir + "/tof." + FPar + "digibdf.par");
+   TObjString *tofDigiBdfFile = new TObjString( workDir  + "/parameters/tof/tof_" + TofGeo +".digibdf.par");
+   parFileList->Add(tofDigiBdfFile);
+
+   TString geoDir  = gSystem->Getenv("VMCWORKDIR");
+   TString geoFile = geoDir + "/geometry/tof/geofile_tof_" + TofGeo + ".root";
+   TFile* fgeo = new TFile(geoFile);
+   TGeoManager *geoMan = (TGeoManager*) fgeo->Get("FAIRGeom");
+   if (NULL == geoMan){
+     cout << "<E> FAIRGeom not found in geoFile"<<endl;
+     return;
+   }
+   
+   if(0){
+   TGeoVolume* master=geoMan->GetTopVolume();
+   master->SetVisContainers(1); 
+   master->Draw("ogl"); 
+   }
+
+  // -----   Reconstruction run   -------------------------------------------
+  FairRunAna *run= new FairRunAna();
+  run->SetInputFile(InputFile.Data());
+  //run->SetOutputFile(OutputFile);
+   run->SetUserOutputFileName(OutputFile.Data());
+   run->SetSink(new FairRootFileSink(run->GetUserOutputFileName()));
+
+   CbmTofEventClusterizer* tofClust = new CbmTofEventClusterizer("TOF Event Clusterizer",iVerbose, bOut);
+
+   tofClust->SetCalMode(calMode);
+   tofClust->SetCalSel(calSel);
+   tofClust->SetCaldXdYMax(50.);         // geometrical matching window in cm 
+   tofClust->SetCalCluMulMax(4.);        // Max Counter Cluster Multiplicity for filling calib histos  
+   tofClust->SetCalRpc(calSm);           // select detector for calibration update  
+   tofClust->SetTRefId(RefSel);          // reference trigger for offset calculation 
+   tofClust->SetTotMax(20.);             // Tot upper limit for walk corection
+   tofClust->SetTotMin(0.01);            //(12000.);  // Tot lower limit for walk correction
+   tofClust->SetTotPreRange(2.);         // effective lower Tot limit  in ns from peak position
+   tofClust->SetTotMean(2.);             // Tot calibration target value in ns 
+   tofClust->SetMaxTimeDist(1.0);        // default cluster range in ns 
+   //tofClust->SetMaxTimeDist(0.);       //Deb// default cluster range in ns 
+   tofClust->SetDelTofMax(60.);          // acceptance range for cluster correlation in cm (!)  
+   tofClust->SetSel2MulMax(4);           // limit Multiplicity in 2nd selector
+   tofClust->SetChannelDeadtime(dDeadtime);    // artificial deadtime in ns 
+   tofClust->SetEnableAvWalk(kTRUE);
+   tofClust->SetYFitMin(1.E4);
+   //tofClust->SetTimePeriod(6.25);       // ignore coarse time 
+   //tofClust->SetCorMode(2);              // correct missing hits
+
+   Int_t calSelRead = calSel;
+   if (calSel<0) calSelRead=0;
+   TString cFname=Form("%s_set%09d_%02d_%01dtofClust.hst.root",cFileId.Data(),iCalSet,calMode,calSelRead);
+   tofClust->SetCalParFileName(cFname);
+   TString cOutFname=Form("tofClust_%s_set%09d.hst.root",cFileId.Data(),iCalSet);
+   tofClust->SetOutHstFileName(cOutFname);
+
+   TString cAnaFile=Form("%s_%09d%03d_tofAnaCosmic.hst.root",cFileId.Data(),iCalSet,iSel2);
+      
+   switch (calMode) {
+   case -1:                                    // initial calibration 
+     tofClust->SetTotMax(256.);         // range in bin number 
+     tofClust->SetTotPreRange(256.);
+     //tofClust->SetTotMin(1.);
+     tofClust->SetTRefDifMax(50000.);   // in ns 
+     tofClust->PosYMaxScal(10000.);     // in % of length 
+     tofClust->SetMaxTimeDist(0.);      // no cluster building  
+     //tofClust->SetTimePeriod(0.);       // inspect coarse time 
+     break;
+    case 0:                                    // initial calibration 
+     tofClust->SetTotMax(256.);         // range in bin number 
+     tofClust->SetTotPreRange(256.);
+     //tofClust->SetTotMin(1.);
+     tofClust->SetTRefDifMax(50.);   // in ns 
+     tofClust->PosYMaxScal(10.);     // in % of length 
+     tofClust->SetMaxTimeDist(0.);      // no cluster building  
+     //tofClust->SetTimePeriod(0.);       // inspect coarse time 
+     break;
+   case 1:                                      // save offsets, update walks, for diamonds 
+     tofClust->SetTotMax(256.);           // range in bin number 
+     tofClust->SetTotPreRange(256.);
+     tofClust->SetTRefDifMax(6.25);       // in ns 
+     //tofClust->SetTimePeriod(6.25);       // inspect coarse time 
+     tofClust->PosYMaxScal(10.);       // in % of length 
+     break;
+   case 11:
+     tofClust->SetTRefDifMax(5.);       // in ns 
+     tofClust->PosYMaxScal(3.0);        // in % of length 
+     break;   
+   case 21:
+     tofClust->SetTRefDifMax(2.5);      // in ns 
+     tofClust->PosYMaxScal(2.0);        // in % of length 
+     break;
+   case 31:
+     tofClust->SetTRefDifMax(2.);    // in ns 
+     tofClust->PosYMaxScal(1.5);        // in % of length 
+     break;
+   case 41:
+     tofClust->SetTRefDifMax(1.);    // in ns 
+     tofClust->PosYMaxScal(0.8);        // in % of length 
+     break;   
+   case 51:
+     tofClust->SetTRefDifMax(0.7);     // in ns 
+     tofClust->PosYMaxScal(0.7);        // in % of length 
+     break;
+   case 61:
+     tofClust->SetTRefDifMax(0.5);      // in ns 
+     tofClust->PosYMaxScal(0.7);        // in % of length 
+     break;   
+   case 71:
+     tofClust->SetTRefDifMax(0.4);     // in ns 
+     tofClust->PosYMaxScal(0.6);        // in % of length 
+     break;
+
+   case 2:                                    // time difference calibration
+     tofClust->SetTRefDifMax(300.);     // in ns 
+     tofClust->PosYMaxScal(1000.);      //in % of length
+     break;
+
+   case 3:                                    // time offsets 
+     tofClust->SetTRefDifMax(200.);     // in ns 
+     tofClust->PosYMaxScal(100.);       //in % of length
+     tofClust->SetMaxTimeDist(0.);      // no cluster building  
+     break;
+   case 12:
+   case 13:
+     tofClust->SetTRefDifMax(100.);   // in ns 
+     tofClust->PosYMaxScal(50.);        //in % of length
+     break;
+   case 22:
+   case 23:
+     tofClust->SetTRefDifMax(50.);    // in ns 
+     tofClust->PosYMaxScal(20.);         //in % of length
+     break;
+   case 32:
+   case 33:
+     tofClust->SetTRefDifMax(25.);    // in ns 
+     tofClust->PosYMaxScal(10.);         //in % of length
+     break;
+   case 42:
+   case 43:
+     tofClust->SetTRefDifMax(12.);   // in ns 
+     tofClust->PosYMaxScal(5.);        //in % of length
+     break;
+   case 52:
+   case 53:
+     tofClust->SetTRefDifMax(5.);   // in ns 
+     tofClust->PosYMaxScal(3.);        //in % of length
+     break;
+   case 62:
+   case 63:
+     tofClust->SetTRefDifMax(3.);   // in ns 
+     tofClust->PosYMaxScal(2.);        //in % of length
+     break;
+   case 72:
+   case 73:
+     tofClust->SetTRefDifMax(2.);    // in ns 
+     tofClust->PosYMaxScal(1.5);        //in % of length
+     break;
+   case 82:
+   case 83:
+     tofClust->SetTRefDifMax(1.);    // in ns 
+     tofClust->PosYMaxScal(1.0);        //in % of length   
+     break;
+   case 92:
+   case 93:
+     tofClust->SetTRefDifMax(0.6);    // in ns 
+     tofClust->PosYMaxScal(1.0);        //in % of length   
+     break;
+
+   case 4:                                      // velocity dependence (DelTOF)
+     tofClust->SetTRefDifMax(6.);    // in ns 
+     tofClust->PosYMaxScal(1.5);        //in % of length
+     break;
+   case 14:
+     tofClust->SetTRefDifMax(5.);   // in ns 
+     tofClust->PosYMaxScal(1.);        //in % of length
+     break;
+   case 24:
+     tofClust->SetTRefDifMax(3.);   // in ns 
+     tofClust->PosYMaxScal(1.0);        //in % of length
+     break;
+   case 34:
+     tofClust->SetTRefDifMax(2.);   // in ns 
+     tofClust->PosYMaxScal(1.0);        //in % of length
+     break;
+   case 44:
+     tofClust->SetTRefDifMax(1.);   // in ns 
+     tofClust->PosYMaxScal(1.0);        //in % of length
+     break;
+   case 54:
+     tofClust->SetTRefDifMax(0.7);     // in ns 
+     tofClust->PosYMaxScal(0.7);        //in % of length
+     break;
+   case 64:
+     tofClust->SetTRefDifMax(0.5);     // in ns 
+     tofClust->PosYMaxScal(0.7);        //in % of length
+     break;
+   default:
+     cout << "<E> Calib mode not implemented! stop execution of script"<<endl;
+     return;
+   } 
+   run->AddTask(tofClust);
+
+   Int_t iBRef=iCalSet%1000;
+   Int_t iSet = (iCalSet - iBRef)/1000;
+   Int_t iRSel=0;
+   Int_t iRSelTyp=0;
+   Int_t iRSelSm=0;
+   Int_t iRSelRpc=0;
+   iRSel=iBRef;     // use diamond
+   if(iSel2==0){
+     // iSel2=iBRef;
+   }else{
+     if(iSel2<0) iSel2=-iSel2;
+   }   
+
+   Int_t iRSelin=iRSel;
+   iRSelRpc=iRSel%10;
+   iRSelTyp = (iRSel-iRSelRpc)/10;
+   iRSelSm=iRSelTyp%10;
+   iRSelTyp = (iRSelTyp-iRSelSm)/10;
+
+   tofClust->SetBeamRefId(iRSelTyp);    // define Beam reference counter 
+   tofClust->SetBeamRefSm(iRSelSm);
+   tofClust->SetBeamRefDet(iRSelRpc);
+   tofClust->SetBeamAddRefMul(-1);
+   tofClust->SetBeamRefMulMax(3);
+
+   Int_t iSel2in=iSel2;
+   Int_t iSel2Rpc= iSel2%10;
+   iSel2=(iSel2-iSel2Rpc)/10;
+   Int_t iSel2Sm=iSel2%10;
+   iSel2=(iSel2-iSel2Sm)/10;
+   if(iSel2 > 0) {
+     tofClust->SetSel2Id(iSel2); 
+     tofClust->SetSel2Sm(iSel2Sm); 
+     tofClust->SetSel2Rpc(iSel2Rpc);
+   }
+
+   Int_t iRef = iSet %1000;
+   Int_t iDut = (iSet - iRef)/1000;
+   Int_t iDutRpc = iDut%10;
+   iDut = (iDut - iDutRpc)/10;
+   Int_t iDutSm = iDut%10;
+   iDut = (iDut - iDutSm)/10;
+
+   tofClust->SetDutId(iDut);
+   tofClust->SetDutSm(iDutSm);
+   tofClust->SetDutRpc(iDutRpc);
+
+   Int_t iRefRpc = iRef%10;
+   iRef = (iRef - iRefRpc)/10;
+   Int_t iRefSm = iRef%10;
+   iRef = (iRef - iRefSm)/10;
+
+   tofClust->SetSelId(iRef);
+   tofClust->SetSelSm(iRefSm);
+   tofClust->SetSelRpc(iRefRpc);
+
+   cout << "Run with iRSel = "<<iRSel<<", iSel2 = "<<iSel2in<<endl;
+   
+   // =========================================================================
+   // ===                       Tracking                                    ===
+   // =========================================================================
+   CbmStsDigitize* stsDigitize = new CbmStsDigitize(); //necessary for kalman !!
+   CbmKF* kalman = new CbmKF();
+
+   CbmTofTrackFinder* tofTrackFinder= new CbmTofTrackFinderNN();
+   tofTrackFinder->SetMaxTofTimeDifference(5.);    // in ns/cm 
+   tofTrackFinder->SetTxLIM(1.6);                 // max slope dx/dz
+   tofTrackFinder->SetTyLIM(1.6);                 // max dev from mean slope dy/dz
+   tofTrackFinder->SetTyMean(0.);                 // mean slope dy/dz
+   tofTrackFinder->SetMaxTofTimeDifference(1.);   // in ns/cm
+ 
+   CbmTofTrackFitter* tofTrackFitter= new CbmTofTrackFitterKF(0,211);
+   TFitter *MyFit = new TFitter(1);                // initialize Minuit
+   tofTrackFinder->SetFitter(tofTrackFitter);
+   CbmTofFindTracks* tofFindTracks  = new CbmTofFindTracks("TOF Track Finder");
+   tofFindTracks->UseFinder(tofTrackFinder);
+   tofFindTracks->UseFitter(tofTrackFitter);
+   tofFindTracks->SetCorMode(iGenCor);           // valid options: 0,1,2,3,4,5,6, 10 - 19
+   tofFindTracks->SetTtTarg(0.035);              // target value for inverse velocity, > 0.033 ns/cm!
+   tofFindTracks->SetCalParFileName(cTrkFile);   // Tracker parameter value file name
+   //   tofFindTracks->SetBeamCounter(5,0,0);         // default beam counter 
+   //tofFindTracks->SetBeamCounter(9,0,1);
+   //tofFindTracks->SetBeamCounter(1,1,1);
+  
+   tofFindTracks->SetT0MAX(dScalFac);            // in ns
+   tofFindTracks->SetSIGT(0.08);                 // default in ns
+   tofFindTracks->SetSIGX(0.3);                  // default in cm
+   tofFindTracks->SetSIGY(0.6);                  // default in cm 
+   tofFindTracks->SetSIGZ(0.1);                  // default in cm 
+   tofFindTracks->SetUseSigCalib(bUseSigCalib);        // ignore resolutions in CalPar file
+   tofTrackFinder->SetSIGLIM(dChi2Lim2*2.);      // matching window in multiples of chi2
+   tofTrackFinder->SetChiMaxAccept(dChi2Lim2);   // max tracklet chi2
+
+   Int_t iMinNofHits=-1;
+   Int_t iNStations=0;
+   Int_t iNReqStations=3;
+
+   switch (iTrackingSetup){
+   case 0:            // bypass mode
+     iMinNofHits=-1;   
+     iNStations=1;        
+     tofFindTracks->SetStation(0, 9, 0, 0);       // Diamond 
+     break;
+   case 1:                                           // for calibration mode of full setup
+     iMinNofHits=3;
+     iNStations=10;
+     iNReqStations=4;
+     //     tofFindTracks->SetStation(0, 9, 0, 0);         
+     tofFindTracks->SetStation(0, 9, 1, 0);
+     tofFindTracks->SetStation(1, 9, 2, 0);         
+     tofFindTracks->SetStation(2, 9, 0, 0);         
+     tofFindTracks->SetStation(3, 9, 1, 1);         
+     tofFindTracks->SetStation(4, 9, 2, 1);
+     tofFindTracks->SetStation(5, 9, 0, 1);     
+     tofFindTracks->SetStation(6, 6, 0, 0);         
+     tofFindTracks->SetStation(7, 6, 0, 1);
+     tofFindTracks->SetStation(8, 6, 1, 0);
+     tofFindTracks->SetStation(9, 6, 1, 1);
+     /*
+     tofFindTracks->SetStation(8, 8, 0, 0);       
+     tofFindTracks->SetStation(9, 8, 0, 1);       
+     tofFindTracks->SetStation(10, 8, 0, 2);       
+     tofFindTracks->SetStation(11, 8, 0, 3);       
+     tofFindTracks->SetStation(12, 8, 0, 4);       
+     tofFindTracks->SetStation(13, 8, 0, 5);       
+     tofFindTracks->SetStation(14, 8, 0, 6);       
+     tofFindTracks->SetStation(15, 8, 0, 7);       
+     */
+     break;
+
+   case 11:                                           // for debugging
+     iMinNofHits=3;
+     iNStations=5;
+     iNReqStations=5;
+     //     tofFindTracks->SetStation(0, 9, 0, 0);         
+     tofFindTracks->SetStation(0, 9, 1, 0);
+     tofFindTracks->SetStation(1, 9, 2, 0);         
+     tofFindTracks->SetStation(2, 9, 2, 1);
+     tofFindTracks->SetStation(3, 9, 0, 1);     
+     tofFindTracks->SetStation(4, 1, 0, 1);
+     break;
+     
+   case 2:
+     iMinNofHits=6;
+     iNStations=7;
+     iNReqStations=7;
+     tofFindTracks->SetStation(0, 9, 1, 0);
+     tofFindTracks->SetStation(1, 9, 2, 0);         
+     tofFindTracks->SetStation(2, 9, 1, 1);         
+     tofFindTracks->SetStation(3, 9, 2, 1);  
+     tofFindTracks->SetStation(4, 9, 0, 0);         
+     tofFindTracks->SetStation(5, 9, 0, 1);
+     tofFindTracks->SetStation(6, 1, 0, 1);
+     //     tofFindTracks->SetStation(6, iDut, iDutSm, iDutRpc);
+     break;
+     
+   case 21:
+     iMinNofHits=5;
+     iNStations=6;
+     iNReqStations=6;
+     tofFindTracks->SetStation(0, 9, 1, 0);
+     tofFindTracks->SetStation(1, 9, 2, 0);         
+     tofFindTracks->SetStation(2, 9, 1, 1);         
+     tofFindTracks->SetStation(3, 9, 2, 1);  
+     //tofFindTracks->SetStation(4, 9, 0, 0);         
+     tofFindTracks->SetStation(4, 9, 0, 1);
+     tofFindTracks->SetStation(5, 1, 0, 1);
+     //     tofFindTracks->SetStation(6, iDut, iDutSm, iDutRpc);
+     break;
+     
+   case 3:
+     iMinNofHits=4;
+     iNStations=5;
+     iNReqStations=5; 
+     tofFindTracks->SetStation(0, 9, 1, 1);         
+     tofFindTracks->SetStation(1, 9, 2, 1);  
+     tofFindTracks->SetStation(2, 9, 1, 0);
+     tofFindTracks->SetStation(3, 9, 2, 0);
+     tofFindTracks->SetStation(4, 9, 0, 1);  
+     break;
+
+   case 4:
+     iMinNofHits=5;
+
+     break;
+
+   case 5:     // for calibration of 2-stack and add-on counters (STAR2, CERN)
+     iMinNofHits=5;
+     iNStations=7;
+     iNReqStations=7;
+     //     tofFindTracks->SetStation(0, 9, 0, 0);         
+     tofFindTracks->SetStation(0, 9, 1, 0);
+     tofFindTracks->SetStation(1, 9, 2, 0);         
+     tofFindTracks->SetStation(2, 9, 1, 1);         
+     tofFindTracks->SetStation(3, 9, 2, 1);
+     tofFindTracks->SetStation(4, 9, 0, 1);         
+     tofFindTracks->SetStation(5, 6, 0, 0);         
+     tofFindTracks->SetStation(6, 6, 0, 1);
+     break;
+     
+   default:
+     cout << "Tracking setup "<<iTrackingSetup<<" not implemented "<<endl;
+     return;
+     ;
+   }
+   tofFindTracks->SetMinNofHits(iMinNofHits);
+   tofFindTracks->SetNStations(iNStations);
+   tofFindTracks->SetNReqStations(iNReqStations);
+   tofFindTracks->PrintSetup();
+   run->AddTask(tofFindTracks);
+
+   // =========================================================================
+   // ===                       Analysis                                    ===
+   // =========================================================================
+
+   CbmTofAnaTestbeam* tofAnaTestbeam = new CbmTofAnaTestbeam("TOF TestBeam Analysis",iVerbose);
+
+   //CbmTofAnaTestbeam defaults  
+   tofAnaTestbeam->SetReqTrg(0);   // 0 - no selection
+   tofAnaTestbeam->SetDXMean(0.);
+   tofAnaTestbeam->SetDYMean(0.);
+   tofAnaTestbeam->SetDTMean(0.);  // in ps
+   tofAnaTestbeam->SetDXWidth(0.4);
+   tofAnaTestbeam->SetDYWidth(0.4);
+   tofAnaTestbeam->SetDTWidth(80.); // in ps
+   tofAnaTestbeam->SetCalParFileName(cAnaFile);
+   tofAnaTestbeam->SetPosY4Sel(0.5); // Y Position selection in fraction of strip length
+   tofAnaTestbeam->SetDTDia(0.);   // Time difference to additional diamond
+   tofAnaTestbeam->SetCorMode(RefSel); // 1 - DTD4, 2 - X4 
+   tofAnaTestbeam->SetMul0Max(30);      // Max Multiplicity in dut 
+   tofAnaTestbeam->SetMul4Max(30);      // Max Multiplicity in Ref - RPC 
+   tofAnaTestbeam->SetMulDMax(10);      // Max Multiplicity in Diamond    
+   tofAnaTestbeam->SetHitDistMin(30.);  // initialization
+   tofAnaTestbeam->SetEnableMatchPosScaling(kFALSE);
+
+   tofAnaTestbeam->SetPosYS2Sel(0.5);   // Y Position selection in fraction of strip length
+   tofAnaTestbeam->SetChS2Sel(0.);      // Center of channel selection window
+   tofAnaTestbeam->SetDChS2Sel(100.);   // Width  of channel selection window
+   tofAnaTestbeam->SetTShift(0.);       // Shift DTD4 to 0
+   tofAnaTestbeam->SetSel2TOff(0.);     // Shift Sel2 time peak to 0 
+   tofAnaTestbeam->SetTOffD4(13.);   // Shift DTD4 to physical value
+
+   if(iSel2 > 0) {
+     tofAnaTestbeam->SetMrpcSel2(iSel2);           // initialization of second selector Mrpc Type 
+     tofAnaTestbeam->SetMrpcSel2Sm(iSel2Sm);       // initialization of second selector Mrpc SmId
+     tofAnaTestbeam->SetMrpcSel2Rpc(iSel2Rpc);     // initialization of second selector Mrpc RpcId
+   }
+
+   tofAnaTestbeam->SetDut(iDut);              // Device under test   
+   tofAnaTestbeam->SetDutSm(iDutSm);          // Device under test   
+   tofAnaTestbeam->SetDutRpc(iDutRpc);        // Device under test   
+   tofAnaTestbeam->SetMrpcRef(iRef);          // Reference RPC     
+   tofAnaTestbeam->SetMrpcRefSm(iRefSm);      // Reference RPC     
+   tofAnaTestbeam->SetMrpcRefRpc(iRefRpc);    // Reference RPC     
+
+   tofAnaTestbeam->SetChi2Lim(10.);           // initialization of Chi2 selection limit  
+
+   if(0)
+   switch (iSet) {
+   case 0:                                 // upper part of setup: P2 - P5
+   case 3:                                 // upper part of setup: P2 - P5
+   case 34:                                // upper part of setup: P2 - P5
+   case 400300:
+     switch (iRSel){
+         case 4:
+	   tofAnaTestbeam->SetTShift(0.);   // Shift DTD4 to 0
+	   tofAnaTestbeam->SetTOffD4(16.);   // Shift DTD4 to physical value
+	   tofAnaTestbeam->SetSel2TOff(0.);     // Shift Sel2 time peak to 0
+	   break;
+	   
+         case 5:
+	   tofAnaTestbeam->SetTShift(-3.);     // Shift DTD4 to 0
+	   tofAnaTestbeam->SetTOffD4(16.);   // Shift DTD4 to physical value
+	   tofAnaTestbeam->SetSel2TOff(0.);     // Shift Sel2 time peak to 0
+	   break;
+
+         case 9:
+	   tofAnaTestbeam->SetChi2Lim(100.);   // initialization of Chi2 selection limit  
+	   tofAnaTestbeam->SetMulDMax(3);      // Max Multiplicity in BeamRef // Diamond    
+	   tofAnaTestbeam->SetTShift(0.1);   // Shift DTD4 to 0
+	   tofAnaTestbeam->SetTOffD4(16.);   // Shift DTD4 to physical value
+	   tofAnaTestbeam->SetSel2TOff(0.5);     // Shift Sel2 time peak to 0
+	   break;
+
+         default:
+	   ;
+     }
+
+     default:
+	   cout<<"<E> detector setup "<<iSet<<" unknown, stop!"<<endl;
+	   return;
+	 ;
+   }  // end of different subsets
+
+   //run->AddTask(tofAnaTestbeam);
+   // =========================================================================
+   /*
+   CbmTofOnlineDisplay* display = new CbmTofOnlineDisplay();
+   display->SetUpdateInterval(1000);
+   run->AddTask(display);   
+   */
+  // -----  Parameter database   --------------------------------------------
+ 
+   FairRuntimeDb* rtdb = run->GetRuntimeDb();
+   Bool_t kParameterMerged = kTRUE;
+   FairParRootFileIo* parIo2 = new FairParRootFileIo(kParameterMerged);
+   parIo2->open(ParFile.Data(), "UPDATE");
+   parIo2->print();
+   rtdb->setFirstInput(parIo2);
+  
+   FairParAsciiFileIo* parIo1 = new FairParAsciiFileIo();
+   parIo1->open(parFileList, "in");
+   parIo1->print();
+   rtdb->setSecondInput(parIo1);
+   rtdb->print();
+   rtdb->printParamContexts();
+
+   //  FairParRootFileIo* parInput1 = new FairParRootFileIo();
+   //  parInput1->open(ParFile.Data());
+   //  rtdb->setFirstInput(parInput1);
+
+   FairEventManager *fMan= new FairEventManager();
+ 
+  CbmEvDisTracks *Tracks =  new CbmEvDisTracks("Tof Tracks",1,kFALSE,kTRUE); //name, verbosity, RnrChildren points, RnrChildren track
+   //  CbmEvDisTracks *Tracks =  new CbmEvDisTracks("Tof Tracks",1);
+  fMan->AddTask(Tracks);
+  CbmPixelHitSetDraw *TofUHits = new CbmPixelHitSetDraw ("TofUHit", kRed, kOpenCross );
+  fMan->AddTask(TofUHits);
+   CbmPointSetArrayDraw *TofHits = new CbmPointSetArrayDraw ("TofHit", 1, 1, 4, kTRUE); //name, colorMode, markerMode, verbosity, RnrChildren
+   //  CbmPixelHitSetDraw *TofHits = new CbmPixelHitSetDraw ("TofHit", kRed, kOpenCircle, 4);// kFullSquare);
+  fMan->AddTask(TofHits); 
+
+
+  TGeoVolume* top = gGeoManager->GetTopVolume();
+  gGeoManager->SetVisOption(1);
+  gGeoManager->SetVisLevel(5);
+  TObjArray* allvolumes = gGeoManager->GetListOfVolumes();
+  //cout<<"GeoVolumes  "  << gGeoManager->GetListOfVolumes()->GetEntries()<<endl;
+  for(Int_t i=0; i<allvolumes->GetEntries(); i++){
+    TGeoVolume* vol     = (TGeoVolume*)allvolumes->At(i);
+    TString name = vol->GetName();
+    //    cout << " GeoVolume "<<i<<" Name: "<< name << endl;
+    vol->SetTransparency(90);
+    /* switch (char *) not allowed any more in root 6 :(
+    switch(name.Data()) {
+    case "counter":
+      vol->SetTransparency(95);
+      break;
+
+    case "tof_glass":
+    case "Gap":
+    case "Cell":
+      vol->SetTransparency(99);
+      break;
+
+    case "pcb":
+      vol->SetTransparency(30);
+      break;
+
+    default:
+      vol->SetTransparency(96);
+    }
+    */
+  }
+  //  gGeoManager->SetVisLevel(3);
+  //  top->SetTransparency(80); 
+  //  top->Draw("ogl");
+
+  //  fMan->Init(1,4,10000);                     
+  fMan->Init(1,5); 
+                    
+  cout << "customize TEveManager gEve "<< gEve << endl;
+  gEve->GetDefaultGLViewer()->SetClearColor(kYellow-10);
+  TGLViewer* v = gEve->GetDefaultGLViewer();
+  TGLAnnotation* ann = new TGLAnnotation(v,cFileId,0.01,0.98);
+  ann->SetTextSize(0.03);// % of window diagonal
+  ann->SetTextColor(4);
+
+  //  gEve->TEveProjectionAxes()->SetDrawOrigin(kTRUE);
+
+  {   // from readCurrentCamera(const char* fname)
+  TGLCamera& c = gEve->GetDefaultGLViewer()->CurrentCamera();
+  const char* fname="Cam.sav";
+  TFile* f = TFile::Open(fname, "READ");
+  if (!f) 
+    return;
+  if (f->GetKey(c.ClassName())) {
+    f->GetKey(c.ClassName())->Read(&c);
+    c.IncTimeStamp();
+    gEve->GetDefaultGLViewer()->RequestDraw();
+  }
+  }
+
+  // -----   Intialise and run   --------------------------------------------
+  // run->Init();
+  //  cout << "Starting run" << endl;
+  //  run->Run(0, nEvents);
+  // ------------------------------------------------------------------------
+  // default display 
+  /*
+  TString Display_Status = "pl_over_Mat04D4best.C";
+  TString Display_Funct = "pl_over_Mat04D4best()";  
+  gROOT->LoadMacro(Display_Status);
+
+  gROOT->LoadMacro("fit_ybox.h");
+  gROOT->LoadMacro("pl_all_CluMul.C");
+  gROOT->LoadMacro("pl_all_CluRate.C");
+  gROOT->LoadMacro("pl_over_cluSel.C");
+  gROOT->LoadMacro("pl_over_clu.C");
+  gROOT->LoadMacro("pl_all_dTSel.C");
+  gROOT->LoadMacro("pl_over_MatD4sel.C");
+  gROOT->LoadMacro("save_hst.C");
+
+  switch(iSet){
+    default:
+  case 0:
+  case 3:
+  case 49:
+  case 79:
+  case 34:
+  case 94:
+  case 37:
+  case 97:
+  case 39:
+  case 99:
+  case 93:
+  case 300400:
+  case 400300:
+  case 910900:
+  case 300900:
+  case 400900:
+  case 901900:
+  case 921920:
+  case 300921:
+  case 920921:
+  case 920300:
+  case 921300:
+
+    gInterpreter->ProcessLine("pl_over_clu(6)");
+    gInterpreter->ProcessLine("pl_over_clu(6,0,1)");
+    gInterpreter->ProcessLine("pl_over_clu(9,0,0)");
+    gInterpreter->ProcessLine("pl_over_clu(9,0,1)");
+    gInterpreter->ProcessLine("pl_over_clu(9,1,0)");
+    gInterpreter->ProcessLine("pl_over_clu(9,1,1)");
+    gInterpreter->ProcessLine("pl_over_clu(9,2,0)");
+    gInterpreter->ProcessLine("pl_over_clu(9,2,1)");
+    gInterpreter->ProcessLine("pl_over_cluSel(0,6,0,0)");
+    gInterpreter->ProcessLine("pl_over_cluSel(0,6,0,1)");
+    gInterpreter->ProcessLine("pl_over_cluSel(0,9,0,0)");
+    gInterpreter->ProcessLine("pl_over_cluSel(0,9,0,1)");
+    gInterpreter->ProcessLine("pl_over_cluSel(0,9,1,0)");
+    gInterpreter->ProcessLine("pl_over_cluSel(0,9,1,1)");
+    gInterpreter->ProcessLine("pl_over_cluSel(0,9,2,0)");
+    gInterpreter->ProcessLine("pl_over_cluSel(0,9,2,1)");
+    gInterpreter->ProcessLine("pl_over_cluSel(1,6,0,0)");
+    gInterpreter->ProcessLine("pl_over_cluSel(1,6,0,1)");
+    gInterpreter->ProcessLine("pl_over_cluSel(1,9,0,0)");
+    gInterpreter->ProcessLine("pl_over_cluSel(1,9,0,1)");
+    gInterpreter->ProcessLine("pl_over_cluSel(1,9,1,0)");
+    gInterpreter->ProcessLine("pl_over_cluSel(1,9,1,1)");
+    gInterpreter->ProcessLine("pl_over_cluSel(1,9,2,0)");
+    gInterpreter->ProcessLine("pl_over_cluSel(1,9,2,1)");
+    gInterpreter->ProcessLine("pl_all_CluMul()");
+    gInterpreter->ProcessLine("pl_all_CluRate()");
+    gInterpreter->ProcessLine("pl_all_dTSel()");
+    TString FSave=Form("save_hst(\"cosdev-status%d_%d_Cal_%s.hst.root\")",iCalSet,iSel2in,cCalId.Data());
+    gInterpreter->ProcessLine(FSave.Data());
+    //gInterpreter->ProcessLine("pl_over_MatD4sel()");
+    break;
+    ;
+  }
+  gInterpreter->ProcessLine(Display_Funct);
+  */
+}
diff --git a/macro/beamtime/hd2020/dis_trks.C b/macro/beamtime/hd2020/dis_trks.C
new file mode 100755
index 0000000000000000000000000000000000000000..ef57d51f3732525374c822a2313269d66317a54c
--- /dev/null
+++ b/macro/beamtime/hd2020/dis_trks.C
@@ -0,0 +1,570 @@
+void dis_trks(Int_t nEvents=10, Int_t iSel=1, Int_t iGenCor=1, TString cFileId="48.50.7.1", TString cSet="000010020", Int_t iSel2=20, Int_t iTrackingSetup=2, Double_t dScalFac=1., Double_t dChi2Lim2=500., Double_t dDeadtime=50,TString cCalId="", Int_t iAnaCor=1,  Int_t iCalSet=30040500) 
+{
+   Int_t iVerbose = 1;
+   if( cCalId == "" ) cCalId=cFileId;
+   // Specify log level (INFO, DEBUG, DEBUG1, ...)
+   //TString logLevel = "FATAL";
+   //TString logLevel = "ERROR";
+   //TString logLevel = "INFO";
+   TString logLevel = "DEBUG"; 
+   //TString logLevel = "DEBUG1";
+   //TString logLevel = "DEBUG2";
+   //TString logLevel = "DEBUG3";
+   TString workDir       = gSystem->Getenv("VMCWORKDIR");
+   //TString workDir          = "../../..";
+   //TString paramDir       = workDir  + "/macro/beamtime/mcbm2020";
+   TString paramDir       = "./";
+   TString ParFile           = paramDir + "/data/" + cFileId.Data() + ".params.root";
+   TString InputFile        = paramDir + "/data/" + cFileId.Data() + ".root";
+   TString InputDigiFile = paramDir + "/data/digidev_" + cFileId.Data() + Form("_%s_%02.0f_Cal",cSet.Data(),dDeadtime) + cCalId + ".out.root";
+   TString OutputFile    = paramDir + "/data/dishits_" + cFileId.Data() + Form("_%s_%06d_%03d",cSet.Data(),iSel,iSel2) + ".out.root";
+   TString cHstFile=paramDir + Form("/hst/%s_%03.0f_%s_%06d_%03d_%03.1f_%03.1f_trk%03d_Cal%s_Dis.hst.root",cFileId.Data(),dDeadtime,cSet.Data(),iSel,iSel2,dScalFac,dChi2Lim2,iTrackingSetup,cCalId.Data());
+   TString cTrkFile=Form("%s_tofFindTracks.hst.root",cCalId.Data());
+   TString cAnaFile=Form("%s_TrkAnaTestBeam.hst.root",cCalId.Data());
+
+   cout << " InputDigiFile = "	<< InputDigiFile   << endl;
+
+   TString shcmd = "rm -v " + ParFile;
+   gSystem->Exec( shcmd.Data()  );
+
+   TList *parFileList = new TList();
+
+   TString TofGeo="v20b_mcbm";  //default
+   cout << "Geometry version "<<TofGeo<<endl;
+
+     // -----   Load the geometry setup   -------------------------------------
+   /*
+  const char* setupName = "mcbm_beam_2020_03";
+  TString setupFile = workDir + "/geometry/setup/setup_" + setupName + ".C";
+  TString setupFunct = "setup_";
+  setupFunct = setupFunct + setupName + "()";
+  std::cout << "-I- mcbm_reco: Loading macro " << setupFile << std::endl;
+  gROOT->LoadMacro(setupFile);
+  gROOT->ProcessLine(setupFunct);
+  CbmSetup* setup = CbmSetup::Instance();
+*/
+  
+   TString geoDir  = gSystem->Getenv("VMCWORKDIR");
+   TString geoFile = geoDir + "/geometry/tof/geofile_tof_" + TofGeo + ".root";
+   TFile* fgeo = new TFile(geoFile);
+   TGeoManager *geoMan = (TGeoManager*) fgeo->Get("FAIRGeom");
+   if (NULL == geoMan){
+     cout << "<E> FAIRGeom not found in geoFile"<<endl;
+     return;
+   }
+  
+   //TObjString *tofDigiFile = new TObjString(workDir + "/parameters/tof/tof_" + TofGeo + ".digi.par"); // TOF digi file
+   //parFileList->Add(tofDigiFile);   
+
+   TObjString *tofDigiBdfFile = new TObjString(workDir  + "/parameters/tof/" + TofGeo +".digibdf.par");
+   parFileList->Add(tofDigiBdfFile);
+
+  // -----   Reconstruction run   -------------------------------------------
+  FairRunAna *run= new FairRunAna();
+  cout << "InputFile:     "<<InputFile.Data()<<endl;
+  cout << "InputDigiFile: "<<InputDigiFile.Data()<<endl;
+
+  //run->SetInputFile(InputFile.Data());
+  //run->AddFriend(InputDigiFile.Data());
+  run->SetInputFile(InputDigiFile.Data());
+  //run->AddFriend(InputFile.Data());
+  run->SetOutputFile(OutputFile);
+
+  FairLogger::GetLogger()->SetLogScreenLevel(logLevel.Data());
+  //  FairLogger::GetLogger()->SetLogVerbosityLevel("MEDIUM");
+  FairLogger::GetLogger()->SetLogVerbosityLevel("VERYHIGH");
+
+  // -----   Local selection variables  -------------------------------------------
+
+   Int_t iRef = iSel %1000;
+   Int_t iDut = (iSel - iRef)/1000;
+   Int_t iDutRpc = iDut%10;
+   iDut = (iDut - iDutRpc)/10;
+   Int_t iDutSm = iDut%10;
+   iDut = (iDut - iDutSm)/10;
+   Int_t iRefRpc = iRef%10;
+   iRef = (iRef - iRefRpc)/10;
+   Int_t iRefSm = iRef%10;
+   iRef = (iRef - iRefSm)/10;
+
+   Int_t iSel2in = iSel2;
+   Int_t iSel2Rpc= iSel2%10;
+   iSel2=(iSel2-iSel2Rpc)/10;
+   Int_t iSel2Sm=iSel2%10;
+   iSel2=(iSel2-iSel2Sm)/10;
+   
+   Int_t calMode=93;
+   Int_t calSel=1;
+   Bool_t bOut=kFALSE;
+   
+   CbmTofEventClusterizer* tofClust = new CbmTofEventClusterizer("TOF Event Clusterizer",iVerbose, bOut);
+   Int_t calSelRead = calSel;
+   if (calSel<0) calSelRead=0;
+   TString cFname=Form("%s_set%09d_%02d_%01dtofClust.hst.root",cFileId.Data(),iCalSet,calMode,calSelRead);
+   if(cCalId!="XXX") 
+     cFname=Form("%s_set%09d_%02d_%01dtofClust.hst.root",cCalId.Data(),iCalSet,calMode,calSelRead);
+   tofClust->SetCalParFileName(cFname);
+   TString cOutFname=Form("tofClust_%s_set%09d.hst.root",cFileId.Data(),iCalSet);
+   tofClust->SetOutHstFileName(cOutFname);
+   
+   // =========================================================================
+   // ===                       Tracking                                    ===
+   // =========================================================================
+   /* 
+   CbmTofEventClusterizer* tofClust = new CbmTofEventClusterizer("TOF Event Clusterizer",iVerbose, bOut);
+   tofClust->SetMemoryTime(1000000.);    // internal storage time of hits in ns
+   */
+
+   CbmTofTrackFinder* tofTrackFinder= new CbmTofTrackFinderNN();
+   tofTrackFinder->SetMaxTofTimeDifference(0.4);    // in ns/cm
+   Int_t TrackerPar=0;
+   switch(TrackerPar) {
+   case 0:                                                          // for full mTof setup
+     tofTrackFinder->SetTxLIM(0.3);                 // max slope dx/dz
+     tofTrackFinder->SetTyLIM(0.3);                 // max dev from mean slope dy/dz
+     tofTrackFinder->SetTxMean(0.);                 // mean slope dy/dz
+     tofTrackFinder->SetTyMean(0.);                 // mean slope dy/dz
+     break;
+   case 1:                                                            // for double stack test counters
+     tofTrackFinder->SetTxMean(0.21);                 // mean slope dy/dz
+     tofTrackFinder->SetTyMean(0.18);                 // mean slope dy/dz
+     tofTrackFinder->SetTxLIM(0.15);                 // max slope dx/dz
+     tofTrackFinder->SetTyLIM(0.18);                 // max dev from mean slope dy/dz
+     break;
+   }
+   
+   CbmTofTrackFitter* tofTrackFitter= new CbmTofTrackFitterKF(0,211);
+   TFitter *MyFit = new TFitter(1);                // initialize Minuit
+   tofTrackFinder->SetFitter(tofTrackFitter);
+   CbmTofFindTracks* tofFindTracks  = new CbmTofFindTracks("TOF Track Finder");
+   tofFindTracks->UseFinder(tofTrackFinder);
+   tofFindTracks->UseFitter(tofTrackFitter);
+   tofFindTracks->SetCalOpt(1);                  // 1 - update offsets
+   tofFindTracks->SetCorMode(iGenCor);           // valid options: 0,1,2,3,4,5,6, 10 - 19
+   tofFindTracks->SetTtTarg(0.057);               // target value for inverse velocity, > 0.033 ns/cm!
+   //tofFindTracks->SetTtTarg(0.035);            // target value for inverse velocity, > 0.033 ns/cm!
+   tofFindTracks->SetCalParFileName(cTrkFile);   // Tracker parameter value file name
+   tofFindTracks->SetBeamCounter(5,0,0);         // default beam counter 
+   tofFindTracks->SetStationMaxHMul(30);          // Max Hit Multiplicity in any used station
+  
+   tofFindTracks->SetT0MAX(dScalFac);            // in ns
+   tofFindTracks->SetSIGT(0.08);                 // default in ns
+   tofFindTracks->SetSIGX(0.3);                  // default in cm
+   tofFindTracks->SetSIGY(0.6);                  // default in cm 
+   tofFindTracks->SetSIGZ(0.05);                 // default in cm 
+   tofFindTracks->SetUseSigCalib(kFALSE);        // ignore resolutions in CalPar file
+   tofTrackFinder->SetSIGLIM(dChi2Lim2*2.);      // matching window in multiples of chi2
+   tofTrackFinder->SetChiMaxAccept(dChi2Lim2);   // max tracklet chi2
+
+
+   Int_t iMinNofHits=-1;
+   Int_t iNStations=0;
+   Int_t iNReqStations=3;
+
+   switch (iTrackingSetup){
+   case 0:            // bypass mode
+     iMinNofHits=-1;   
+     iNStations=1;        
+     tofFindTracks->SetStation(0, 5, 0, 0);       // Diamond 
+     break;
+
+   case 1:                                           // for calibration mode of full setup
+     iMinNofHits=3;
+     iNStations=30;
+     iNReqStations=4;
+     tofFindTracks->SetStation(0, 5, 0, 0);          
+     tofFindTracks->SetStation(1, 0, 2, 2);          
+     tofFindTracks->SetStation(2, 0, 1, 2);
+     tofFindTracks->SetStation(3, 0, 0, 2);         
+     tofFindTracks->SetStation(4, 0, 2, 1);     
+     tofFindTracks->SetStation(5, 0, 1, 1);
+     tofFindTracks->SetStation(6, 0, 0, 1);         
+     tofFindTracks->SetStation(7, 0, 2, 3);          
+     tofFindTracks->SetStation(8, 0, 1, 3);
+     tofFindTracks->SetStation(9, 0, 0, 3);         
+     tofFindTracks->SetStation(10, 0, 2, 0);          
+     tofFindTracks->SetStation(11, 0, 1, 0);
+     tofFindTracks->SetStation(12, 0, 0, 0);         
+     tofFindTracks->SetStation(13, 0, 2, 4);          
+     tofFindTracks->SetStation(14, 0, 1, 4);
+     tofFindTracks->SetStation(15, 0, 0, 4);    
+     tofFindTracks->SetStation(16, 0, 4, 0);          
+     tofFindTracks->SetStation(17, 0, 3, 0);
+     tofFindTracks->SetStation(18, 0, 4, 1);          
+     tofFindTracks->SetStation(19, 0, 3, 1);
+     tofFindTracks->SetStation(20, 0, 4, 2);          
+     tofFindTracks->SetStation(21, 0, 3, 2);
+     tofFindTracks->SetStation(22, 0, 4, 3);          
+     tofFindTracks->SetStation(23, 0, 3, 3);
+     tofFindTracks->SetStation(24, 0, 4, 4);          
+     tofFindTracks->SetStation(25, 0, 3, 4);
+     tofFindTracks->SetStation(26, 9, 0, 0);         
+     tofFindTracks->SetStation(27, 9, 0, 1);         
+     tofFindTracks->SetStation(28, 6, 0, 0);     
+     tofFindTracks->SetStation(29, 6, 0, 1);              
+     break;
+
+   case 2:
+     iMinNofHits=3;
+     iNStations=14;
+     iNReqStations=5;
+     tofFindTracks->SetStation(0, 5, 0, 0);          
+     tofFindTracks->SetStation(1, 0, 4, 1);          
+     tofFindTracks->SetStation(2, 0, 3, 1);
+     tofFindTracks->SetStation(3, 0, 4, 0);          
+     tofFindTracks->SetStation(4, 0, 3, 0);
+     tofFindTracks->SetStation(5, 0, 4, 2);          
+     tofFindTracks->SetStation(6, 0, 3, 2);
+     tofFindTracks->SetStation(7, 0, 4, 3);          
+     tofFindTracks->SetStation(8, 0, 3, 3);
+     tofFindTracks->SetStation(9, 0, 4, 4);          
+     tofFindTracks->SetStation(10, 0, 3, 4);
+     tofFindTracks->SetStation(11, 9, 0, 0);         
+     tofFindTracks->SetStation(12, 9, 0, 1);         
+     tofFindTracks->SetStation(13, 7, 0, 0);         
+     break;
+
+   case 3:
+     iMinNofHits=3;
+     iNStations=16;
+     iNReqStations=4;
+     tofFindTracks->SetStation(0, 5, 0, 0);          
+     tofFindTracks->SetStation(1, 0, 2, 2);          
+     tofFindTracks->SetStation(2, 0, 1, 2);
+     tofFindTracks->SetStation(3, 0, 0, 2);         
+
+     tofFindTracks->SetStation(4, 0, 2, 1);          
+     tofFindTracks->SetStation(5, 0, 1, 1);
+     tofFindTracks->SetStation(6, 0, 0, 1);         
+
+     tofFindTracks->SetStation(7, 0, 2, 3);          
+     tofFindTracks->SetStation(8, 0, 1, 3);
+     tofFindTracks->SetStation(9, 0, 0, 3);         
+
+     tofFindTracks->SetStation(10, 0, 2, 0);          
+     tofFindTracks->SetStation(11, 0, 1, 0);
+     tofFindTracks->SetStation(12, 0, 0, 0);         
+
+     tofFindTracks->SetStation(13, 0, 2, 4);          
+     tofFindTracks->SetStation(14, 0, 1, 4);
+     tofFindTracks->SetStation(15, 0, 0, 4);         
+      
+     /*
+     tofFindTracks->SetStation(16, 0, 3, 2);         
+     tofFindTracks->SetStation(17, 0, 4, 2);  
+     tofFindTracks->SetStation(18, 0, 3, 1);         
+     tofFindTracks->SetStation(19, 0, 4, 1);
+     tofFindTracks->SetStation(20, 0, 3, 3);         
+     tofFindTracks->SetStation(21, 0, 4, 3);
+     tofFindTracks->SetStation(22, 0, 3, 0);         
+     tofFindTracks->SetStation(23, 0, 4, 0);
+     tofFindTracks->SetStation(24, 0, 3, 4);         
+     tofFindTracks->SetStation(25, 0, 4, 4); 
+     */
+     break;
+
+   case 4:
+     iMinNofHits=3;
+     iNStations=4;
+     iNReqStations=4;
+     tofFindTracks->SetStation(0, 5, 0, 0);          
+     tofFindTracks->SetStation(1,  0, 4, 1);           
+     tofFindTracks->SetStation(2,  0, 3, 1);            
+     tofFindTracks->SetStation(3, iDut, iDutSm, iDutRpc);      
+     break;
+
+   case 5:     // for calibration of 2-stack and add-on counters (STAR2, BUC)
+     iMinNofHits=3;
+     iNStations=5;
+     iNReqStations=4;
+     tofFindTracks->SetStation(0, 5, 0, 0);          
+     tofFindTracks->SetStation(1, 0, 4, 1);          
+     tofFindTracks->SetStation(2, 0, 3, 1);
+     tofFindTracks->SetStation(3, 9, 0, 0);         
+     tofFindTracks->SetStation(4, 9, 0, 1);          
+     break;
+
+   case 6:
+     iMinNofHits=3;
+     iNStations=4;
+     iNReqStations=4;
+     tofFindTracks->SetStation(0, 5, 0, 0);          
+     tofFindTracks->SetStation(1, 0, 4, 1);          
+     tofFindTracks->SetStation(2, 0, 3, 1);
+     tofFindTracks->SetStation(3, iDut, iDutSm, iDutRpc);         
+     //    tofFindTracks->SetStation(3, 9, 0, 0);         
+     //    tofFindTracks->SetStation(3, 9, 0, 1);         
+     //    tofFindTracks->SetStation(3, 7, 0, 0);         
+     break;
+
+   case 7:     // for calibration of 2-stack and add-on counters (BUC)
+     iMinNofHits=4;
+     iNStations=5;
+     iNReqStations=5;
+     tofFindTracks->SetStation(0, 5, 0, 0);          
+     tofFindTracks->SetStation(1, 0, 4, 3);          
+     tofFindTracks->SetStation(2, 0, 3, 3);
+     tofFindTracks->SetStation(3, 6, 0, 0);         
+     tofFindTracks->SetStation(4, 6, 0, 1);         
+     break;
+
+   case 8:     // evaluation of add-on counters (BUC)
+     iMinNofHits=3;
+     iNStations=4;
+     iNReqStations=4;
+     tofFindTracks->SetStation(0, 5, 0, 0);          
+     tofFindTracks->SetStation(1, 0, 4, 3);          
+     tofFindTracks->SetStation(2, 0, 3, 3);
+     tofFindTracks->SetStation(3, iDut, iDutSm, iDutRpc);         
+     break;
+
+  case 10:
+     iMinNofHits=3;
+     iNStations=4;
+     iNReqStations=4;
+     tofFindTracks->SetStation(0, 5, 0, 0);          
+     tofFindTracks->SetStation(1, 0, 1, 2);          
+     tofFindTracks->SetStation(2, 0, 0, 2);
+     tofFindTracks->SetStation(3, 0, 2, 2);    
+   
+   default:
+     cout << "Tracking setup "<<iTrackingSetup<<" not implemented "<<endl;
+     return;
+     ;
+   }
+   tofFindTracks->SetMinNofHits(iMinNofHits);
+   tofFindTracks->SetNStations(iNStations);
+   tofFindTracks->SetNReqStations(iNReqStations);
+   tofFindTracks->PrintSetup();
+   run->AddTask(tofFindTracks);
+
+   // =========================================================================
+   // ===                       Analysis                                    ===
+   // =========================================================================
+   CbmTofAnaTestbeam* tofAnaTestbeam = new CbmTofAnaTestbeam("TOF TestBeam Analysis",iVerbose);
+   tofAnaTestbeam->SetCorMode(iAnaCor); // 1 - DTD4, 2 - X4, 3 - Y4, 4 - Texp 
+   tofAnaTestbeam->SetHitDistMin(30.);  // initialization
+   tofAnaTestbeam->SetEnableMatchPosScaling(kTRUE);
+   tofAnaTestbeam->SetSpillDuration(4.);
+   //CbmTofAnaTestbeam defaults  
+   tofAnaTestbeam->SetDXMean(0.);
+   tofAnaTestbeam->SetDYMean(0.);
+   tofAnaTestbeam->SetDTMean(0.);      // in ns
+   tofAnaTestbeam->SetDXWidth(0.5);
+   tofAnaTestbeam->SetDYWidth(1.0);
+   tofAnaTestbeam->SetDTWidth(0.1);    // in ns
+   tofAnaTestbeam->SetCalParFileName(cAnaFile);
+   Double_t dScalFacA=0.9;  // dScalFac is used for tracking 
+   tofAnaTestbeam->SetPosY4Sel(0.5*dScalFacA);   // Y Position selection in fraction of strip length
+   tofAnaTestbeam->SetDTDia(0.);      // Time difference to additional diamond
+   tofAnaTestbeam->SetMul0Max(20);    // Max Multiplicity in dut 
+   tofAnaTestbeam->SetMul4Max(30);     // Max Multiplicity in Ref - RPC 
+   tofAnaTestbeam->SetMulDMax(3);     // Max Multiplicity in Diamond / BeamRef    
+   tofAnaTestbeam->SetTOffD4(14.);    // initialization
+   tofAnaTestbeam->SetDTD4MAX(6.);    // initialization of Max time difference Ref - BRef
+
+   //tofAnaTestbeam->SetTShift(-28000.);// initialization
+   tofAnaTestbeam->SetPosYS2Sel(0.55);    // Y Position selection in fraction of strip length
+   tofAnaTestbeam->SetChS2Sel(0.);             // Center of channel selection window
+   tofAnaTestbeam->SetDChS2Sel(100.);      // Width  of channel selection window
+   tofAnaTestbeam->SetSel2TOff(0.);           // Shift Sel2 time peak to 0 
+   tofAnaTestbeam->SetChi2Lim(5.);             // initialization of Chi2 selection limit  
+   tofAnaTestbeam->SetChi2Lim2(3.);           // initialization of Chi2 selection limit for Mref-Sel2 pair   
+   tofAnaTestbeam->SetDutDX(15.);              // limit inspection of tracklets to selected region
+   tofAnaTestbeam->SetDutDY(15.);              // limit inspection of tracklets to selected region
+   tofAnaTestbeam->SetSIGLIM(3.);                // max matching chi2
+   tofAnaTestbeam->SetSIGT(0.08);                // in ns
+   tofAnaTestbeam->SetSIGX(0.3);                  // in cm
+   tofAnaTestbeam->SetSIGY(0.6);                  // in cm
+
+   Int_t iRSel=500;
+   Int_t iRSelTyp=5;
+   Int_t iRSelSm=0;
+   Int_t iRSelRpc=0;
+   Int_t iRSelin=iRSel;
+   
+   tofAnaTestbeam->SetBeamRefSmType(iRSelTyp); // common reaction reference 
+   tofAnaTestbeam->SetBeamRefSmId(iRSelSm);
+   tofAnaTestbeam->SetBeamRefRpc(iRSelRpc);
+
+   if(iSel2 >= 0) {
+     tofAnaTestbeam->SetMrpcSel2(iSel2);                   // initialization of second selector Mrpc Type 
+     tofAnaTestbeam->SetMrpcSel2Sm(iSel2Sm);       // initialization of second selector Mrpc SmId
+     tofAnaTestbeam->SetMrpcSel2Rpc(iSel2Rpc);             // initialization of second selector Mrpc RpcId
+   }
+
+   tofAnaTestbeam->SetDut(iDut);              // Device under test   
+   tofAnaTestbeam->SetDutSm(iDutSm);          // Device under test   
+   tofAnaTestbeam->SetDutRpc(iDutRpc);        // Device under test   
+   tofAnaTestbeam->SetMrpcRef(iRef);          // Reference RPC     
+   tofAnaTestbeam->SetMrpcRefSm(iRefSm);      // Reference RPC     
+   tofAnaTestbeam->SetMrpcRefRpc(iRefRpc);    // Reference RPC  
+
+   cout<< "dispatch iSel = "<<iSel<<", iSel2in = "<<iSel2in<<", iRSelin = "<<iRSelin<<", iRSel = "<<iRSel<<endl;
+   if(1)
+   {
+   switch (iSel) {
+
+   case 10:   
+     switch (iRSelin){
+       case 500:	 
+	   tofAnaTestbeam->SetTShift(2.5);    // Shift DTD4 to 0
+	   tofAnaTestbeam->SetTOffD4(18.);   // Shift DTD4 to physical value
+	   switch(iSel2in){
+	   case 20:
+	     tofAnaTestbeam->SetSel2TOff(0.);     // Shift Sel2 time peak to 0
+	     break;
+	   default:
+	     ;   
+	   }
+	   break;
+        default:
+       ;
+     }
+     break;
+
+   case 700040:   
+   case 900040:   
+   case 901040:   
+     switch (iRSelin){
+       case 500:	 
+	   tofAnaTestbeam->SetTShift(0.3);    // Shift DTD4 to 0
+	   tofAnaTestbeam->SetTOffD4(18.);   // Shift DTD4 to physical value
+
+	   switch(iSel2in){
+	   case 30:
+	     tofAnaTestbeam->SetSel2TOff(-0.3);     // Shift Sel2 time peak to 0
+	     break;
+	   case 31:
+	     tofAnaTestbeam->SetSel2TOff(-0.41);     // Shift Sel2 time peak to 0
+	     break;
+
+	   default:
+	     ;   
+	   }
+	   break;
+        default:
+       ;
+     }
+     break;
+
+   case 700041:   
+   case 900041:   
+   case 901041:   
+     switch (iRSelin){
+       case 500:	 
+	   tofAnaTestbeam->SetTShift(5.3);    // Shift DTD4 to 0
+	   tofAnaTestbeam->SetTOffD4(18.);   // Shift DTD4 to physical value
+
+	   switch(iSel2in){
+	   case 30:
+	     tofAnaTestbeam->SetSel2TOff(-0.3);     // Shift Sel2 time peak to 0
+	     break;
+	   case 31:
+	     tofAnaTestbeam->SetSel2TOff(-0.41);     // Shift Sel2 time peak to 0
+	     break;
+
+	   default:
+	     ;   
+	   }
+	   break;
+        default:
+       ;
+     }
+     break;
+
+   case 600043:   
+   case 601043:   
+     switch (iRSelin){
+     case 500:	 
+	   tofAnaTestbeam->SetTShift(5.3);    // Shift DTD4 to 0
+	   tofAnaTestbeam->SetTOffD4(11.);   // Shift DTD4 to physical value
+
+	   switch(iSel2in){
+	     case 33:
+	       tofAnaTestbeam->SetSel2TOff(-0.55);     // Shift Sel2 time peak to 0
+	       break;
+
+	     default:
+	       ;   
+	   }
+	   break;
+        default:
+       ;
+     }
+     break;
+
+   default:
+     cout << "Better to define analysis setup! Running with default offset parameter... "<< endl;
+     // return;
+   }  // end of different subsets
+
+   cout << " Initialize TSHIFT to "<<tofAnaTestbeam->GetTShift()<<endl;
+   //run->AddTask(tofAnaTestbeam);
+   }
+
+  // -----  Parameter database   --------------------------------------------
+   FairRuntimeDb* rtdb = run->GetRuntimeDb();
+   Bool_t kParameterMerged = kTRUE;
+   FairParRootFileIo* parIo2 = new FairParRootFileIo(kParameterMerged);
+   parIo2->open(ParFile.Data(), "UPDATE");
+   parIo2->print();
+   rtdb->setFirstInput(parIo2);
+   
+   FairParAsciiFileIo* parIo1 = new FairParAsciiFileIo();
+   parIo1->open(parFileList, "in");
+   parIo1->print();
+   rtdb->setSecondInput(parIo1);
+   rtdb->print();
+   rtdb->printParamContexts();
+
+   //  FairParRootFileIo* parInput1 = new FairParRootFileIo();
+   //  parInput1->open(ParFile.Data());
+   //  rtdb->setFirstInput(parInput1);
+
+   FairEventManager *fMan= new FairEventManager();
+ 
+  CbmEvDisTracks *Tracks =  new CbmEvDisTracks("Tof Tracks",1,kFALSE,kTRUE); //name, verbosity, RnrChildren points, RnrChildren track
+   //  CbmEvDisTracks *Tracks =  new CbmEvDisTracks("Tof Tracks",1);
+  fMan->AddTask(Tracks);
+  CbmPixelHitSetDraw *TofUHits = new CbmPixelHitSetDraw ("TofUHit", kRed, kOpenCross );
+  fMan->AddTask(TofUHits);
+   CbmPointSetArrayDraw *TofHits = new CbmPointSetArrayDraw ("TofHit", 1, 1, 1, kTRUE); //name, colorMode, markerMode, verbosity, RnrChildren
+   //  CbmPixelHitSetDraw *TofHits = new CbmPixelHitSetDraw ("TofHit", kRed, kOpenCircle, 4);// kFullSquare);
+  fMan->AddTask(TofHits); 
+
+  TGeoVolume* top = gGeoManager->GetTopVolume();
+  gGeoManager->SetVisOption(1);
+  gGeoManager->SetVisLevel(5);
+  TObjArray* allvolumes = gGeoManager->GetListOfVolumes();
+  //cout<<"GeoVolumes  "  << gGeoManager->GetListOfVolumes()->GetEntries()<<endl;
+  for(Int_t i=0; i<allvolumes->GetEntries(); i++){
+    TGeoVolume* vol     = (TGeoVolume*)allvolumes->At(i);
+    //TString name = vol->GetName();
+    //cout << " GeoVolume "<<i<<" Name: "<< name << endl;
+    vol->SetLineColor(kRed);
+    vol->SetTransparency(80);
+  }
+  fMan->Init(1,5); 
+                    
+  cout << "customize TEveManager gEve "<< gEve << endl;
+  gEve->GetDefaultGLViewer()->SetClearColor(kYellow-10);
+  TGLViewer* v = gEve->GetDefaultGLViewer();
+  TGLAnnotation* ann = new TGLAnnotation(v,cFileId,0.01,0.98);
+  ann->SetTextSize(0.03);// % of window diagonal
+  ann->SetTextColor(4);
+
+  {   // from readCurrentCamera(const char* fname)
+  TGLCamera& c = gEve->GetDefaultGLViewer()->CurrentCamera();
+  const char* fname="Cam.sav";
+  TFile* f = TFile::Open(fname, "READ");
+  if (!f) 
+    return;
+  if (f->GetKey(c.ClassName())) {
+    f->GetKey(c.ClassName())->Read(&c);
+    c.IncTimeStamp();
+    gEve->GetDefaultGLViewer()->RequestDraw();
+  }
+  }
+
+}
diff --git a/macro/beamtime/hd2020/pl_all_2D.C b/macro/beamtime/hd2020/pl_all_2D.C
new file mode 100644
index 0000000000000000000000000000000000000000..6f8e24889228009b3992595780d627e86881f930
--- /dev/null
+++ b/macro/beamtime/hd2020/pl_all_2D.C
@@ -0,0 +1,93 @@
+void pl_all_2D(Int_t iOpt=0, Int_t iNSt=2)
+{
+  //  TCanvas *can = new TCanvas("can22","can22");
+  //  can->Divide(2,2); 
+  //  TCanvas *can = new TCanvas("can","can",48,55,700,900);
+  TCanvas *can = new TCanvas("can","can",48,56,900,900);
+  can->Divide(4,3,0.01,0.01); 
+  //  can->Divide(2,2,0,0); 
+  Float_t lsize=0.07;
+
+ gPad->SetFillColor(0);
+ gStyle->SetPalette(1);
+ gStyle->SetLabelSize(lsize);
+
+ //gStyle->SetOptStat(kTRUE);
+ //gROOT->cd();
+ //gROOT->SetDirLevel(2);
+
+ TH2 *h;
+ TH2 *h2;
+ const Int_t   iType[6]   ={9,6,7,1,6,8};
+ const Int_t  iSmNum[6]={3,2,1,1,1,1};
+ const Int_t iRpcNum[6]={2,2,1,3,2,8};
+ TString cOpt; 
+
+ switch(iOpt){
+ case 0: 
+   cOpt="Size";
+   break;
+ case 1:
+   cOpt="Pos";
+   break;
+ case 2:
+   cOpt="TOff";
+   break;
+ case 3:
+   cOpt="Tot";
+   break;
+ case 4:
+   cOpt="AvWalk";
+   break;
+ case 5:
+   cOpt="AvLnWalk";
+   break;
+ case 6:
+   cOpt="Mul";
+   break;
+ case 7:
+   cOpt="Trms";
+   break;
+ case 8:
+   cOpt="DelPos";
+   break;
+ case 9:
+   cOpt="DelTOff";
+   break;
+ case 10:
+   cOpt="DelMatPos";
+   break;
+ case 11:
+   cOpt="DelMatTOff";
+   break;
+ default:
+   ;
+ }
+
+ Int_t iCanv=0;
+ // if (h!=NULL) h->Delete();
+
+ for(Int_t iSt=0; iSt<iNSt; iSt++){
+   // cout << "plot station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+  for(Int_t iSm=0; iSm<iSmNum[iSt];iSm++){
+    //cout << "plot module at station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+    for(Int_t iRp=0; iRp<iRpcNum[iSt];iRp++){
+      //cout << "plot rpc at station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+      can->cd(iCanv+1); iCanv++;
+      gROOT->cd();
+      TString hname=Form("cl_SmT%01d_sm%03d_rpc%03d_%s",iType[iSt],iSm,iRp,cOpt.Data());
+      h=(TH2 *)gROOT->FindObjectAny(hname);
+      if (h!=NULL) {
+	if (iOpt==6){h->Draw();}
+	else{
+	gPad->SetLogz();
+	h->Draw("colz");
+	h->ProfileX()->Draw("same");
+	}
+      }else{cout<<"Histogram "<<hname<<" not existing. "<<endl;}
+    }
+  }
+ }
+ can->SaveAs(Form("pl_all_%s.pdf",cOpt.Data()));
+
+} 
diff --git a/macro/beamtime/hd2020/pl_all_3D.C b/macro/beamtime/hd2020/pl_all_3D.C
new file mode 100644
index 0000000000000000000000000000000000000000..9582590d67d502633894b06877dc29a1528b6cb6
--- /dev/null
+++ b/macro/beamtime/hd2020/pl_all_3D.C
@@ -0,0 +1,65 @@
+void pl_all_3D(Int_t iOpt=0, Int_t iSel=0, Int_t iNSt=2)
+{
+  //  TCanvas *can = new TCanvas("can22","can22");
+  //  can->Divide(2,2); 
+  //  TCanvas *can = new TCanvas("can","can",48,55,700,900);
+  TCanvas *can = new TCanvas("can","can",48,56,900,900);
+  can->Divide(4,3,0.01,0.01); 
+  //  can->Divide(2,2,0,0); 
+  Float_t lsize=0.07;
+
+ gPad->SetFillColor(0);
+ gStyle->SetPalette(1);
+ gStyle->SetLabelSize(lsize);
+
+ //gStyle->SetOptStat(kTRUE);
+ //gROOT->cd();
+ //gROOT->SetDirLevel(2);
+
+ TH3 *h;
+ TH2 *h2;
+ const Int_t   iType[6]={9,6,7,5,6,8};
+ const Int_t  iSmNum[6]={3,2,1,1,1,1};
+ const Int_t iRpcNum[6]={2,2,1,1,2,8};
+ TString cOpt; 
+
+ Int_t iCanv=0;
+ // if (h!=NULL) h->Delete();
+
+ for(Int_t iSt=0; iSt<iNSt; iSt++){
+  for(Int_t iSm=0; iSm<iSmNum[iSt];iSm++){
+    for(Int_t iRp=0; iRp<iRpcNum[iSt];iRp++){
+      can->cd(iCanv+1); iCanv++;
+      gROOT->cd();
+      TString hname=Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_Walk2",iType[iSt],iSm,iRp,iSel);
+      h=(TH3 *)gROOT->FindObjectAny(hname);
+      if (h!=NULL) {
+		switch(iOpt){
+		  case 0: 
+			cOpt="yx";
+			h->Project3D(cOpt)->Draw("colz");
+			break;
+
+		  case 1:
+			cOpt="yx";
+			h->Project3D(cOpt)->Draw("colz");
+		    gPad->SetLogz();
+			break;
+
+		  case 2: 
+			cOpt="pfyx";
+			h2 = (TH2 *)h->Project3DProfile("yx");
+			h2->SetMinimum(-0.2);
+			h2->SetMaximum(0.2);
+			h2->Draw("colz");
+			break;
+		  default:
+			;
+		}
+      }else{cout<<"Histogram "<<hname<<" not existing. "<<endl;}
+    }
+  }
+ }
+ can->SaveAs(Form("pl_all_%s.pdf",cOpt.Data()));
+
+} 
diff --git a/macro/beamtime/hd2020/pl_all_CluMul.C b/macro/beamtime/hd2020/pl_all_CluMul.C
new file mode 100644
index 0000000000000000000000000000000000000000..99d1efbbb4612611f793624602ede7415e37269b
--- /dev/null
+++ b/macro/beamtime/hd2020/pl_all_CluMul.C
@@ -0,0 +1,49 @@
+void pl_all_CluMul(Int_t iNSt=2, Double_t MulMax=100)
+{
+  //  TCanvas *can = new TCanvas("can22","can22");
+  //  can->Divide(2,2); 
+  //  TCanvas *can = new TCanvas("can","can",48,55,700,900);
+  TCanvas *can = new TCanvas("can","can",48,56,900,900);
+  can->Divide(4,3,0.01,0.01); 
+  //  can->Divide(2,2,0,0); 
+  Float_t lsize=0.07;
+
+ gPad->SetFillColor(0);
+ gStyle->SetPalette(1);
+ gStyle->SetLabelSize(lsize);
+
+ //gStyle->SetOptStat(kTRUE);
+ //gROOT->cd();
+ //gROOT->SetDirLevel(2);
+
+ TH1 *h;
+ TH2 *h2;
+ const Int_t   iType[6]={9,6,7,5,6,8};
+ const Int_t  iSmNum[6]={3,2,1,1,1,1};
+ const Int_t iRpcNum[6]={2,2,1,1,2,8};
+ 
+ Int_t iCanv=0;
+ // if (h!=NULL) h->Delete();
+
+ for(Int_t iSt=0; iSt<iNSt; iSt++){
+   //   cout << "plot station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+  for(Int_t iSm=0; iSm<iSmNum[iSt];iSm++){
+    // cout << "plot module at station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+   for(Int_t iRp=0; iRp<iRpcNum[iSt];iRp++){
+     //cout << "plot rpc at station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+    can->cd(iCanv+1); iCanv++;
+    gROOT->cd();
+    TString hname=Form("cl_SmT%01d_sm%03d_rpc%03d_Mul",iType[iSt],iSm,iRp);
+    h=(TH1 *)gROOT->FindObjectAny(hname);
+    if (h!=NULL) {
+     h->GetXaxis()->SetRange(0.,MulMax);
+     h->Draw("");
+     gPad->SetLogy();
+    }else{cout<<"Histogram "<<hname<<" not existing. "<<endl;}
+    if(iRp==10) break;
+   }
+  }
+ }
+ can->SaveAs(Form("pl_all_CluMul.pdf"));
+
+} 
diff --git a/macro/beamtime/hd2020/pl_all_CluPosEvol.C b/macro/beamtime/hd2020/pl_all_CluPosEvol.C
new file mode 100644
index 0000000000000000000000000000000000000000..2ab35430b2b5ea1859dc8743c25b5ebde8cd1888
--- /dev/null
+++ b/macro/beamtime/hd2020/pl_all_CluPosEvol.C
@@ -0,0 +1,59 @@
+void pl_all_CluPosEvol(Int_t iNSt=2, Int_t iTmax=0)
+{
+  //  TCanvas *can = new TCanvas("can22","can22");
+  //  can->Divide(2,2); 
+  //  TCanvas *can = new TCanvas("can","can",48,55,700,900);
+  TCanvas *can = new TCanvas("can","can",48,56,900,900);
+  can->Divide(3,4,0.01,0.01); 
+  //  can->Divide(4,4,0.01,0.01); 
+  //  can->Divide(2,2,0,0); 
+  Float_t lsize=0.09;
+
+ gPad->SetFillColor(0);
+ gStyle->SetPalette(1);
+ gStyle->SetLabelSize(lsize);
+
+ //gStyle->SetStatY(0.5);
+ //gStyle->SetStatX(0.5);
+ gStyle->SetStatW(0.5);
+ gStyle->SetStatH(0.3); 
+
+ gStyle->SetOptStat(kFALSE);
+ //gROOT->cd();
+ //gROOT->SetDirLevel(2);
+
+ TH1 *h;
+ TH2 *h2;
+ const Int_t   iType[7]={9,6,1,2,9,8,1};
+ const Int_t  iSmNum[7]={3,2,1,2,3,3,3};
+ const Int_t iRpcNum[7]={2,2,3,1,2,1,3};
+ 
+ Int_t iCanv=0;
+ // if (h!=NULL) h->Delete();
+
+ for(Int_t iSt=0; iSt<iNSt; iSt++){
+   // cout << "plot station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+  for(Int_t iSm=0; iSm<iSmNum[iSt];iSm++){
+    //cout << "plot module at station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+   for(Int_t iRp=0; iRp<iRpcNum[iSt];iRp++){
+     //cout << "plot rpc at station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+    gROOT->cd();
+    TString hname=Form("cl_SmT%01d_sm%03d_rpc%03d_PosEvol",iType[iSt],iSm,iRp);
+    h=(TH1 *)gROOT->FindObjectAny(hname);
+    if (h!=NULL) {
+      can->cd(iCanv+1); iCanv++;
+      Double_t dLMargin=0.35;
+      Double_t dTitOffset=1.6;
+      gPad->SetLeftMargin(dLMargin);
+      h->UseCurrentStyle();
+      h->GetYaxis()->SetTitleOffset(dTitOffset);
+      if(iTmax>0) h->GetXaxis()->SetRange(0.,iTmax);
+      h->Draw("");
+       //gPad->SetLogy();
+    }else{cout<<"Histogram "<<hname<<" not existing. "<<endl;}
+   }
+  }
+ }
+ can->SaveAs(Form("pl_all_CluPosEvol.pdf"));
+
+} 
diff --git a/macro/beamtime/hd2020/pl_all_CluRate.C b/macro/beamtime/hd2020/pl_all_CluRate.C
new file mode 100644
index 0000000000000000000000000000000000000000..d179e2792b11aa234267ef83513a5f70e94144f2
--- /dev/null
+++ b/macro/beamtime/hd2020/pl_all_CluRate.C
@@ -0,0 +1,98 @@
+void pl_all_CluRate(Int_t iNSt=2, Int_t iOpt=0, Double_t Tstart=0., Double_t Tend=800., Int_t iMode=0)
+{
+  //  TCanvas *can = new TCanvas("can22","can22");
+  //  can->Divide(2,2); 
+  //  TCanvas *can = new TCanvas("can","can",48,55,700,900);
+  TCanvas *can = new TCanvas("can","can",48,56,900,900);
+  switch(iMode) {
+  case 0:
+    switch(iNSt){
+    case 6:
+      can->Divide(5,7,0.01,0.01);
+      break;
+    case 5:
+      can->Divide(5,6,0.01,0.01);
+      break;
+    default:
+      can->Divide(5,6,0.01,0.01);
+      break;
+    }
+    break;
+  case 1:
+      can->Divide(1,2,0.01,0.01);
+    ;
+    break;
+  }
+  //  can->Divide(2,2,0,0); 
+  Float_t lsize=0.06;
+
+ gPad->SetFillColor(0);
+ gStyle->SetPalette(1);
+ gStyle->SetLabelSize(lsize);
+
+ //gStyle->SetOptStat(kTRUE);
+ //gROOT->cd();
+ //gROOT->SetDirLevel(2);
+ gStyle->SetPadLeftMargin(0.4);
+ gStyle->SetTitleOffset(1.0, "x");
+ gStyle->SetTitleOffset(1.2, "y");
+ gStyle->SetTitleFontSize(0.03);
+ 
+
+ TH1 *h;
+ TH2 *h2;
+ const Int_t       iType[6]={9,6,7,5,6,8};
+ const Int_t  iSmNum[6]={3,2,1,1,1,1};
+ const Int_t iRpcNum[6]={2,2,1,1,2,8};
+ 
+ Int_t iCanv=0;
+ Int_t iCol=0;
+ // if (h!=NULL) h->Delete();
+
+ for(Int_t iSt=0; iSt<iNSt; iSt++){
+   //   cout << "plot station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+  for(Int_t iSm=0; iSm<iSmNum[iSt];iSm++){
+    //cout << "plot module at station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+   for(Int_t iRp=0; iRp<iRpcNum[iSt];iRp++){
+     //cout << "plot rpc at station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+    gROOT->cd();
+    TString hname="";
+    switch(iOpt){
+    case 0:
+      hname=Form("cl_SmT%01d_sm%03d_rpc%03d_rate",iType[iSt],iSm,iRp);
+      break;
+    case 1:
+      hname=Form("cl_SmT%01d_sm%03d_rpc%03d_rate10s",iType[iSt],iSm,iRp);
+      break;
+    }
+    h=(TH1 *)gROOT->FindObjectAny(hname);
+    if (h!=NULL) {
+     h->GetXaxis()->SetRangeUser(Tstart,Tend);
+     switch(iMode) {
+     case 0:
+       can->cd(iCanv+1); iCanv++;
+       h->Draw("");
+     //h->UseCurrentStyle();
+       gPad->SetLogy();
+       break;
+     case 1:
+       can->cd(iSt+1);
+       if ( iSm>0 ) continue;
+       h->UseCurrentStyle();  // set current defaults 
+       h->SetMarkerStyle(20+iSm);
+       iCol=iRp+1;
+       if(iCol == 5) iCol++;
+       h->SetLineColor(iCol);
+       h->SetLineStyle(1);
+       h->SetMarkerColor(iCol);
+       if(iSm==0 && iRp==0) h->Draw("LPE");
+       else h->Draw("LPEsame");
+       break;
+     }
+    }else{cout<<"Histogram "<<hname<<" not existing. "<<endl;}
+   }
+  }
+ }
+ can->SaveAs(Form("pl_all_CluRate%d_%d.pdf",iNSt,iOpt));
+
+} 
diff --git a/macro/beamtime/hd2020/pl_all_CluTimeEvol.C b/macro/beamtime/hd2020/pl_all_CluTimeEvol.C
new file mode 100644
index 0000000000000000000000000000000000000000..115aa97258c8eb4231816349a319a966fab90724
--- /dev/null
+++ b/macro/beamtime/hd2020/pl_all_CluTimeEvol.C
@@ -0,0 +1,58 @@
+void pl_all_CluTimeEvol(Int_t iNSt=3, Int_t iTmax=0)
+{
+  //  TCanvas *can = new TCanvas("can22","can22");
+  //  can->Divide(2,2); 
+  //  TCanvas *can = new TCanvas("can","can",48,55,700,900);
+  TCanvas *can = new TCanvas("can","can",48,56,900,900);
+  can->Divide(3,4,0.01,0.01); 
+  //  can->Divide(4,4,0.01,0.01); 
+  //  can->Divide(2,2,0,0); 
+  Float_t lsize=0.09;
+
+ gPad->SetFillColor(0);
+ gStyle->SetPalette(1);
+ gStyle->SetLabelSize(lsize);
+
+ //gStyle->SetStatY(0.5);
+ //gStyle->SetStatX(0.5);
+ gStyle->SetStatW(0.5);
+ gStyle->SetStatH(0.3); 
+
+ gStyle->SetOptStat(kFALSE);
+ //gROOT->cd();
+ //gROOT->SetDirLevel(2);
+
+ TH1 *h;
+ TH2 *h2;
+ const Int_t   iType[7]={9,6,1,2,9,8,1};
+ const Int_t  iSmNum[7]={3,1,1,2,3,3,3};
+ const Int_t iRpcNum[7]={2,2,3,1,2,1,3};
+ Int_t iCanv=0;
+ // if (h!=NULL) h->Delete();
+
+ for(Int_t iSt=0; iSt<iNSt; iSt++){
+   // cout << "plot station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+  for(Int_t iSm=0; iSm<iSmNum[iSt];iSm++){
+    //cout << "plot module at station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+   for(Int_t iRp=0; iRp<iRpcNum[iSt];iRp++){
+     //cout << "plot rpc at station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+    gROOT->cd();
+    TString hname=Form("cl_SmT%01d_sm%03d_rpc%03d_TimeEvol",iType[iSt],iSm,iRp);
+    h=(TH1 *)gROOT->FindObjectAny(hname);
+    if (h!=NULL) {
+      can->cd(iCanv+1); iCanv++;
+      Double_t dLMargin=0.35;
+      Double_t dTitOffset=1.6;
+      gPad->SetLeftMargin(dLMargin);
+      h->UseCurrentStyle();
+      h->GetYaxis()->SetTitleOffset(dTitOffset);
+      if(iTmax>0) h->GetXaxis()->SetRange(0.,iTmax);
+      h->Draw("");
+       //gPad->SetLogy();
+    }else{cout<<"Histogram "<<hname<<" not existing. "<<endl;}
+   }
+  }
+ }
+ can->SaveAs(Form("pl_all_CluTimeEvol.pdf"));
+
+} 
diff --git a/macro/beamtime/hd2020/pl_all_DTLastHits.C b/macro/beamtime/hd2020/pl_all_DTLastHits.C
new file mode 100644
index 0000000000000000000000000000000000000000..a994c550a5cba6eba8dba49c4f069fb132d19d26
--- /dev/null
+++ b/macro/beamtime/hd2020/pl_all_DTLastHits.C
@@ -0,0 +1,50 @@
+void pl_all_DTLastHits(Int_t iNSt=6, Double_t Tstart=1., Double_t Tend=1000.)
+{
+  //  TCanvas *can = new TCanvas("can22","can22");
+  //  can->Divide(2,2); 
+  //  TCanvas *can = new TCanvas("can","can",48,55,700,900);
+  TCanvas *can = new TCanvas("can","can",48,56,900,900);
+  can->Divide(5,7,0.01,0.01); 
+  //  can->Divide(2,2,0,0); 
+  Float_t lsize=0.07;
+
+ gPad->SetFillColor(0);
+ gStyle->SetPalette(1);
+ gStyle->SetLabelSize(lsize);
+
+ //gStyle->SetOptStat(kTRUE);
+ //gROOT->cd();
+ //gROOT->SetDirLevel(2);
+
+ TH1 *h;
+ TH2 *h2;
+ const Int_t       iType[6]={0,9,7,5,6,8};
+ const Int_t  iSmNum[6]={5,1,1,1,1,0};
+ const Int_t iRpcNum[6]={5,2,1,1,2,8};
+ 
+ 
+ Int_t iCanv=0;
+ // if (h!=NULL) h->Delete();
+
+ for(Int_t iSt=0; iSt<iNSt; iSt++){
+   //   cout << "plot station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+  for(Int_t iSm=0; iSm<iSmNum[iSt];iSm++){
+    //cout << "plot module at station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+   for(Int_t iRp=0; iRp<iRpcNum[iSt];iRp++){
+     //cout << "plot rpc at station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+    can->cd(iCanv+1); iCanv++;
+    gROOT->cd();
+    TString hname=Form("cl_SmT%01d_sm%03d_rpc%03d_DTLastHits",iType[iSt],iSm,iRp);
+    h=(TH1 *)gROOT->FindObjectAny(hname);
+    if (h!=NULL) {
+     h->GetXaxis()->SetRange(Tstart,Tend);
+     h->Draw("");
+     gPad->SetLogy();
+    }else{cout<<"Histogram "<<hname<<" not existing. "<<endl;}
+    if(iRp==10) break;
+   }
+  }
+ }
+  can->SaveAs(Form("pl_all_CluRate.pdf"));
+
+} 
diff --git a/macro/beamtime/hd2020/pl_all_DigiCor.C b/macro/beamtime/hd2020/pl_all_DigiCor.C
new file mode 100644
index 0000000000000000000000000000000000000000..8b220f6d3a47dceab6d700c0d47b97cfbdc2964c
--- /dev/null
+++ b/macro/beamtime/hd2020/pl_all_DigiCor.C
@@ -0,0 +1,44 @@
+void pl_all_DigiCor(Int_t iNDet=2)
+{
+  //  TCanvas *can = new TCanvas("can22","can22");
+  //  can->Divide(2,2); 
+    TCanvas *can = new TCanvas("can","can",48,55,900,900);
+    //TCanvas *can = new TCanvas("can","can",48,56,900,700);
+    //can->Divide(4,4,0.01,0.01); 
+    //  can->Divide(2,3,0.01,0.01); 
+    can->Divide(4,3,0.01,0.01); 
+  Float_t lsize=0.07;
+
+ gPad->SetFillColor(0);
+ gStyle->SetPalette(1);
+ gStyle->SetLabelSize(lsize);
+
+ //gStyle->SetOptStat(kTRUE);
+ //gROOT->cd();
+ //gROOT->SetDirLevel(2);
+
+ TH1 *h;
+ TH2 *h2;
+ const Int_t   iType[6]={9,6,7,8,6,8};
+ const Int_t  iSmNum[6]={3,2,1,1,1,1};
+ const Int_t iRpcNum[6]={2,2,1,8,2,8};
+
+ Int_t iCanv=0;
+ // if (h!=NULL) h->Delete();
+
+ for(Int_t iCh=0; iCh<iNDet; iCh++){
+  for(Int_t iSm=0; iSm<iSmNum[iCh];iSm++){
+    for(Int_t iRpc=0; iRpc<iRpcNum[iCh];iRpc++){
+      can->cd(iCanv+1); iCanv++;
+      gROOT->cd();
+      TString hname=Form("cl_SmT%01d_sm%03d_rpc%03d_DigiCor",iType[iCh],iSm,iRpc);
+      h2=(TH2 *)gROOT->FindObjectAny(hname);
+      if (h2!=NULL) {
+	h2->Draw("colz");
+	//     gPad->SetLogy();
+      }else{cout<<"Histogram "<<hname<<" not existing. "<<endl;}
+    }
+  }
+ } 
+ can->SaveAs(Form("pl_all_DigiCor.pdf"));
+} 
diff --git a/macro/beamtime/hd2020/pl_all_DigiDTLD.C b/macro/beamtime/hd2020/pl_all_DigiDTLD.C
new file mode 100644
index 0000000000000000000000000000000000000000..51472fd76d39b3e4e93f5ac818b1d1facbe46dad
--- /dev/null
+++ b/macro/beamtime/hd2020/pl_all_DigiDTLD.C
@@ -0,0 +1,129 @@
+void pl_all_DigiDTLD(Int_t iNDet=6, Double_t dDTthr=2., Int_t iOpt=0)
+{
+  //  TCanvas *can = new TCanvas("can22","can22");
+  //  can->Divide(2,2); 
+    TCanvas *can = new TCanvas("can","can",48,55,900,900);
+    //TCanvas *can = new TCanvas("can","can",48,56,900,700);
+    //can->Divide(4,4,0.01,0.01); 
+    //  can->Divide(2,3,0.01,0.01); 
+    can->Divide(5,7,0.01,0.01); 
+  Float_t lsize=0.07;
+
+ gPad->SetFillColor(0);
+ gStyle->SetPalette(1);
+ gStyle->SetLabelSize(lsize);
+
+ //gStyle->SetOptStat(kTRUE);
+ //gROOT->cd();
+ //gROOT->SetDirLevel(2);
+
+ TH1 *h;
+ TH2 *h2;
+ TH1 *hTime;
+ TString hnameT;
+
+ const Int_t       iType[6]={0,9,7,5,6,8};
+ const Int_t  iSmNum[6]={5,1,1,1,1,1};
+ const Int_t iRpcNum[6]={5,2,1,1,2,8};
+
+ Double_t dTime = 0.;
+ Int_t iCanv=0;
+
+ Int_t jSmType=5;
+ Int_t jSm=0;
+ Int_t jRp=0;
+ 
+// if (h!=NULL) h->Delete();
+
+ for(Int_t iCh=0; iCh<iNDet; iCh++){
+  for(Int_t iSm=0; iSm<iSmNum[iCh];iSm++){
+    for(Int_t iRpc=0; iRpc<iRpcNum[iCh];iRpc++){
+      can->cd(iCanv+1); iCanv++;
+      gROOT->cd();
+      TString hname=Form("cl_SmT%01d_sm%03d_rpc%03d_DigiDTLD",iType[iCh],iSm,iRpc);
+      switch(iOpt) {
+      case 0:
+	;
+	break;
+
+      case 1:
+	hname=Form("%s_fdead",hname.Data());
+	break;
+
+      case 2:
+	hname=Form("%s_py",hname.Data());
+	break;
+
+      case 3:
+	hname=Form("%s_py_dead",hname.Data());
+	break;
+
+      default:
+	;
+      }
+
+      h2=(TH2 *)gROOT->FindObjectAny(hname);
+      TH1D* hx;
+      TH1D* hy;
+      TH1D* hy_dead;
+      Double_t dDeadtimeSum=0.;
+      if (h2!=NULL) {
+	switch(iOpt) {
+	 case 0: 	  
+	   h2->Draw("colz");
+	   gPad->SetLogz();
+
+	   // Determine time duration an data taking
+	   hnameT=Form("cl_SmT%01d_sm%03d_rpc%03d_rate",jSmType,jSm,jRp);
+	   hTime=(TH1 *)gROOT->FindObjectAny(hnameT);
+	   for(dTime=0; dTime<hTime->GetNbinsX(); dTime++) 
+	     if(hTime->GetBinContent(dTime+1) == 0) break;
+	   cout << "Normalize for a run duration of " << dTime << " s" << endl;
+
+	   // create result histograms
+	   hx=h2->ProjectionX(Form("%s_fdead",hname.Data()));
+	   hx->GetYaxis()->SetTitle("deadtime fraction");
+	   hx->Reset();
+	   hy=h2->ProjectionY(Form("%s_py",hname.Data()));
+	   hy_dead=h2->ProjectionY(Form("%s_py_dead",hname.Data()));
+	   hy_dead->GetYaxis()->SetTitle("deadtime fraction");
+	   hy_dead->Reset();
+	   for(Int_t iT=hy->GetNbinsX(); iT>0; iT--){
+	     dDeadtimeSum+= hy->GetBinContent(iT)*hy->GetXaxis()->GetBinLowEdge(iT);
+	     Double_t dDeadFrac= dDeadtimeSum / h2->GetNbinsX()/ dTime;
+	      hy_dead->SetBinContent(iT,dDeadFrac);
+	   }
+	  
+	   for (Int_t iCh=0; iCh<h2->GetNbinsX(); iCh++) {
+	     TH1D* hCh=h2->ProjectionY(Form("%s_%d_py",hname.Data(),iCh),iCh+1,iCh+1);
+	     Double_t dAll=hCh->GetEntries();
+	     Double_t dTAllMean=hCh->GetMean();
+	     if(dAll > 0) {
+	       Double_t BL= hCh->GetXaxis()->FindBin(dDTthr);
+	       Double_t dLate=hCh->Integral(BL,hCh->GetNbinsX(),"");
+	       hCh->GetXaxis()->SetRange(BL,hCh->GetNbinsX());
+	       Double_t dTLateMean=hCh->GetMean();
+	       //Double_t dLateRatio=dLate*dTLateMean/dAll/dTAllMean;
+	       Double_t dLateRatio=dLate*dTLateMean/dTime;
+	       cout << Form("Long DT fraction for %s, ch %d: %6.3f, dTAll %6.3f, dTLate %6.3f",hname.Data(),iCh,dLateRatio, dTAllMean, dTLateMean) << endl;
+	       hx->SetBinContent(iCh+1,dLateRatio);
+	     }
+	   }
+	   break;
+	   
+	 case 1:
+	   h2->Draw();
+	   break;
+
+	 case 2:
+	 case 3:
+	   h2->Draw();
+	   gPad->SetLogy();
+	   break;
+ 	} 
+       } else{cout<<"Histogram "<<hname<<" not existing. "<<endl;}
+    }
+  }
+ } 
+ can->SaveAs(Form("pl_all_DigiDTLD.pdf"));
+} 
diff --git a/macro/beamtime/hd2020/pl_all_DigiMul.C b/macro/beamtime/hd2020/pl_all_DigiMul.C
new file mode 100644
index 0000000000000000000000000000000000000000..841a6f39bcca93bfdfadd2852393572379890bd9
--- /dev/null
+++ b/macro/beamtime/hd2020/pl_all_DigiMul.C
@@ -0,0 +1,43 @@
+void pl_all_DigiMul(Int_t iNDet=2)
+{
+  //  TCanvas *can = new TCanvas("can22","can22");
+  //  can->Divide(2,2); 
+    TCanvas *can = new TCanvas("can","can",48,55,900,900);
+    //TCanvas *can = new TCanvas("can","can",48,56,900,700);
+    //can->Divide(4,4,0.01,0.01); 
+    //  can->Divide(2,3,0.01,0.01); 
+    can->Divide(3,3,0.01,0.01); 
+  Float_t lsize=0.07;
+
+ gPad->SetFillColor(0);
+ gStyle->SetPalette(1);
+ gStyle->SetLabelSize(lsize);
+
+ //gStyle->SetOptStat(kTRUE);
+ //gROOT->cd();
+ //gROOT->SetDirLevel(2);
+
+ TH1 *h;
+ TH2 *h2;
+ Int_t iType[6]  ={9,6,6,5,9,8};
+ Int_t iNumSm[6] ={3,2,1,1,3,2};
+ Int_t iNumRpc[6]={2,2,2,1,2,1};
+ Int_t iCanv=0;
+ // if (h!=NULL) h->Delete();
+
+ for(Int_t iCh=0; iCh<iNDet; iCh++){
+  for(Int_t iSm=0; iSm<iNumSm[iCh];iSm++){
+    for(Int_t iRpc=0; iRpc<iNumRpc[iCh];iRpc++){
+      can->cd(iCanv+1); iCanv++;
+      gROOT->cd();
+      TString hname=Form("cl_SmT%01d_sm%03d_rpc%03d_DigiMul",iType[iCh],iSm,iRpc);
+      h2=(TH2 *)gROOT->FindObjectAny(hname);
+      if (h2!=NULL) {
+	h2->Draw("colz");
+	gPad->SetLogz();
+      }else{cout<<"Histogram "<<hname<<" not existing. "<<endl;}
+    }
+  }
+ } 
+ can->SaveAs(Form("pl_all_DigiMul.pdf"));
+} 
diff --git a/macro/beamtime/hd2020/pl_all_DigiStatus.C b/macro/beamtime/hd2020/pl_all_DigiStatus.C
new file mode 100644
index 0000000000000000000000000000000000000000..972538c41b9290317e6fc3dd2741922177011965
--- /dev/null
+++ b/macro/beamtime/hd2020/pl_all_DigiStatus.C
@@ -0,0 +1,43 @@
+void pl_all_DigiStatus(Int_t iNDet=2)
+{
+  //  TCanvas *can = new TCanvas("can22","can22");
+  //  can->Divide(2,2); 
+    TCanvas *can = new TCanvas("can","can",48,55,900,900);
+    //TCanvas *can = new TCanvas("can","can",48,56,900,700);
+    //can->Divide(4,4,0.01,0.01); 
+    //  can->Divide(2,3,0.01,0.01); 
+    can->Divide(4,3,0.01,0.01); 
+  Float_t lsize=0.07;
+
+ gPad->SetFillColor(0);
+ gStyle->SetPalette(1);
+ gStyle->SetLabelSize(lsize);
+
+ //gStyle->SetOptStat(kTRUE);
+ //gROOT->cd();
+ //gROOT->SetDirLevel(2);
+
+ TH1 *h;
+ TH2 *h2;
+ Int_t iType[6]  ={9,6,6,5,9,8};
+ Int_t iNumSm[6] ={3,2,1,1,3,2};
+ Int_t iNumRpc[6]={2,2,2,1,2,1};
+ Int_t iCanv=0;
+ // if (h!=NULL) h->Delete();
+
+ for(Int_t iCh=0; iCh<iNDet; iCh++){
+  for(Int_t iSm=0; iSm<iNumSm[iCh];iSm++){
+    for(Int_t iRpc=0; iRpc<iNumRpc[iCh];iRpc++){
+      can->cd(iCanv+1); iCanv++;
+      gROOT->cd();
+      TString hname=Form("cl_SmT%01d_sm%03d_rpc%03d_DigiStatus",iType[iCh],iSm,iRpc);
+      h2=(TH2 *)gROOT->FindObjectAny(hname);
+      if (h2!=NULL) {
+	h2->Draw("colz");
+	gPad->SetLogz();
+      }else{cout<<"Histogram "<<hname<<" not existing. "<<endl;}
+    }
+  }
+ } 
+ can->SaveAs(Form("pl_all_DigiStatus.pdf"));
+} 
diff --git a/macro/beamtime/hd2020/pl_all_DigiTot.C b/macro/beamtime/hd2020/pl_all_DigiTot.C
new file mode 100644
index 0000000000000000000000000000000000000000..e0e7966a3596f3bea2fce71016b06d7022db1fea
--- /dev/null
+++ b/macro/beamtime/hd2020/pl_all_DigiTot.C
@@ -0,0 +1,44 @@
+void pl_all_DigiTot(Int_t iNDet=2)
+{
+  //  TCanvas *can = new TCanvas("can22","can22");
+  //  can->Divide(2,2); 
+    TCanvas *can = new TCanvas("can","can",48,55,900,900);
+    //TCanvas *can = new TCanvas("can","can",48,56,900,700);
+    //can->Divide(4,4,0.01,0.01); 
+    //  can->Divide(2,3,0.01,0.01); 
+    can->Divide(4,3,0.01,0.01); 
+  Float_t lsize=0.07;
+
+ gPad->SetFillColor(0);
+ gStyle->SetPalette(1);
+ gStyle->SetLabelSize(lsize);
+
+ //gStyle->SetOptStat(kTRUE);
+ //gROOT->cd();
+ //gROOT->SetDirLevel(2);
+
+ TH1 *h;
+ TH2 *h2;
+ const Int_t   iType[6]={9,6,7,8,6,8};
+ const Int_t  iSmNum[6]={3,2,1,1,1,1};
+ const Int_t iRpcNum[6]={2,2,1,8,2,8};
+
+ Int_t iCanv=0;
+ // if (h!=NULL) h->Delete();
+
+ for(Int_t iCh=0; iCh<iNDet; iCh++){
+  for(Int_t iSm=0; iSm<iSmNum[iCh];iSm++){
+    for(Int_t iRpc=0; iRpc<iRpcNum[iCh];iRpc++){
+      can->cd(iCanv+1); iCanv++;
+      gROOT->cd();
+      TString hname=Form("hDigiTot_SmT%01d_sm%03d_rpc%03d",iType[iCh],iSm,iRpc);
+      h2=(TH2 *)gROOT->FindObjectAny(hname);
+      if (h2!=NULL) {
+	h2->Draw("colz");
+	gPad->SetLogz();
+      }else{cout<<"Histogram "<<hname<<" not existing. "<<endl;}
+    }
+  }
+ } 
+ can->SaveAs(Form("pl_all_DigiTot.pdf"));
+} 
diff --git a/macro/beamtime/hd2020/pl_all_Sel2D.C b/macro/beamtime/hd2020/pl_all_Sel2D.C
new file mode 100644
index 0000000000000000000000000000000000000000..0209c6b85ba949f399d0bb152d456500bbafe822
--- /dev/null
+++ b/macro/beamtime/hd2020/pl_all_Sel2D.C
@@ -0,0 +1,122 @@
+void pl_all_Sel2D(Int_t iOpt=0, Int_t iSel=0, Int_t iOpt2=0, Int_t iNSt=2)
+{
+  //  TCanvas *can = new TCanvas("can22","can22");
+  //  can->Divide(2,2); 
+  //  TCanvas *can = new TCanvas("can","can",48,55,700,900);
+  TCanvas *can = new TCanvas("can","can",48,56,900,900);
+  can->Divide(4,3,0.01,0.01); 
+  //  can->Divide(2,2,0,0); 
+  Float_t lsize=0.07;
+
+ gPad->SetFillColor(0);
+ gStyle->SetPalette(1);
+ gStyle->SetLabelSize(lsize);
+
+ //gStyle->SetOptStat(kTRUE);
+ //gROOT->cd();
+ //gROOT->SetDirLevel(2);
+
+ TH1* hp;
+ TH2 *h;
+ TH2 *h2;
+ const Int_t   iType[6]   ={9,6,1,7,6,8};
+ const Int_t  iSmNum[6]={3,2,1,1,1,1};
+ const Int_t iRpcNum[6]={2,2,3,1,2,8};
+ TString cOpt; 
+
+ switch(iOpt){
+ case 0: 
+   cOpt="Size";
+   break;
+ case 1:
+   cOpt="Pos";
+   break;
+ case 2:
+   cOpt="TOff";
+   break;
+ case 3:
+   cOpt="Tot";
+   break;
+ case 4:
+   cOpt="AvWalk";
+   break;
+ case 5:
+   cOpt="DelTof";
+   break;
+ case 6:
+   cOpt="dXdY";
+   break;
+ default:
+   ;
+ }
+
+ Int_t iCanv=0;
+ // if (h!=NULL) h->Delete();
+ TString  FitHName[100];
+ Double_t FitIntegral[100];
+ Double_t FitMean[100];
+ Double_t FitWidth[100];
+ Int_t NFit=0;
+
+ for(Int_t iSt=0; iSt<iNSt; iSt++){
+   // cout << "plot station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+  for(Int_t iSm=0; iSm<iSmNum[iSt];iSm++){
+    //cout << "plot module at station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+    for(Int_t iRp=0; iRp<iRpcNum[iSt];iRp++){
+      //cout << "plot rpc at station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+      can->cd(iCanv+1); iCanv++;
+      gROOT->cd();
+      TString hname=Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_%s",iType[iSt],iSm,iRp,iSel,cOpt.Data());
+      h=(TH2 *)gROOT->FindObjectAny(hname);
+      if (h!=NULL) {
+	h->Draw("colz");
+	gPad->SetLogz();
+	h->ProfileX()->Draw("same");
+	if(iOpt2>0)
+	switch(iOpt){
+	case 6:
+	  {
+	  switch(iOpt2){
+	  case 1: // x-projection
+	    hp=h->ProjectionX();
+	    hp->Draw();
+	    break;
+	  case 2: // y-projection
+	    hp=h->ProjectionY();
+	    hp->Draw();
+	    break;
+	  }
+	  Double_t dFMean=hp->GetMean();
+	  Double_t dFLim=2.5*hp->GetRMS();
+	  if(hp->Integral()>10){
+	    TFitResultPtr fRes=hp->Fit("gaus","S","HEsame",dFMean-dFLim,dFMean+dFLim);
+	    FitHName[NFit]=hp->GetName();
+	    FitIntegral[NFit]=hp->Integral();
+	    FitMean[NFit]= fRes->Parameter(1);
+	    FitWidth[NFit]= fRes->Parameter(2);
+	    NFit++;
+	  ;
+	  }
+	  }
+	  break;  // case 6 end 
+
+	default:
+	  ;
+	}
+
+      }else{cout<<"Histogram "<<hname<<" not existing. "<<endl;}
+      if(iRp==10) break;
+    }
+  }
+ }
+ if(iOpt2>0){
+   for(Int_t iFit=0; iFit<NFit; iFit++){
+	    cout << "FitRes "<< FitHName[iFit] << ", Stat: "<< FitIntegral[iFit]
+		 << ", Mean " << FitMean[iFit]
+		 << ", Width " << FitWidth[iFit] << endl;
+   }
+ }
+
+ can->SaveAs(Form("pl_all_Sel%d_%s.pdf",iSel,cOpt.Data()));
+
+} 
diff --git a/macro/beamtime/hd2020/pl_all_Track2D.C b/macro/beamtime/hd2020/pl_all_Track2D.C
new file mode 100755
index 0000000000000000000000000000000000000000..405a8c1d22a47d5d897308fe913ccdbdf00356c8
--- /dev/null
+++ b/macro/beamtime/hd2020/pl_all_Track2D.C
@@ -0,0 +1,136 @@
+void pl_all_Track2D(Int_t iOpt=1, Int_t iNSt=2)
+{
+  //  TCanvas *can = new TCanvas("can22","can22");
+  //  can->Divide(2,2); 
+  //  TCanvas *can = new TCanvas("can","can",48,55,700,900);
+  TCanvas *can = new TCanvas("can","can",48,56,900,900);
+  can->Divide(4,2,0.01,0.01); 
+  //  can->Divide(2,2,0,0); 
+  Float_t lsize=0.07;
+
+ gPad->SetFillColor(0);
+ gStyle->SetPalette(1);
+ gStyle->SetLabelSize(lsize);
+
+ //gStyle->SetOptStat(kTRUE);
+ //gROOT->cd();
+ //gROOT->SetDirLevel(2);
+
+ TH2 *h;
+ TH2 *h2;
+ const Int_t  iType[6]={9,6,6,5,7,8};
+ const Int_t  iSmNum[6]={3,1,1,1,1,1};
+ const Int_t iRpcNum[6]={2,2,2,1,1,8};
+ TString cOpt; 
+
+ switch(iOpt){
+ case 0: 
+   cOpt="Size";
+   break;
+ case 1:
+   cOpt="Pos";
+   break;
+ case 2:
+   cOpt="TOff";
+   break;
+ case 3:
+   cOpt="Tot";
+   break;
+ case 4:
+   cOpt="Walk";
+   break;
+ case 5:
+   cOpt="Walk";
+   break;
+ case 6:
+   cOpt="Mul";
+   break;
+ case 7:
+   cOpt="Trms";
+   break;
+ case 8:
+   cOpt="DelPos";
+   break;
+ case 9:
+   cOpt="DelTOff";
+   break;
+ case 10:
+   cOpt="DelMatPos";
+   break;
+ case 11:
+   cOpt="DelMatTOff";
+   break;
+ default:
+   ;
+ }
+
+ Int_t iDet=0;
+ Double_t dAvMean=0.;
+ Double_t dAvRMS=0.;
+ Int_t iCanv=0;
+
+ for(Int_t iSt=0; iSt<iNSt; iSt++){
+   // cout << "plot station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+  for(Int_t iSm=0; iSm<iSmNum[iSt];iSm++){
+    //cout << "plot module at station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+    for(Int_t iRp=0; iRp<iRpcNum[iSt];iRp++){
+      //cout << "plot rpc at station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+      can->cd(iCanv+1); iCanv++;
+      gROOT->cd();
+      TString hname="";
+      Int_t iCol=1;
+      switch(iOpt) {
+		case 4:
+		  for (Int_t iSide=0; iSide<2; iSide++)
+		    for (Int_t iCh=0; iCh<32; iCh++) {
+			  hname=Form("cal_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%d_%s",iType[iSt],iSm,iRp,iCh,iSide,cOpt.Data());
+              h=(TH2 *)gROOT->FindObjectAny(hname);
+              if (h!=NULL) {
+                TProfile *hProf=h->ProfileX(Form("%s_pfx%d%d",hname.Data(),iCh,iSide));
+                hProf->SetLineColor(iCol);
+                hProf->SetLineStyle(1);
+                hProf->SetMarkerColor(iCol);
+                hProf->SetMarkerStyle(24+iSide);
+                iCol++;
+                if(iCh==0) iCol=1;
+                if(iCh==0 && iSide==0) {
+				  hProf->SetMaximum(0.4);
+				  hProf->SetMinimum(-0.4);
+				  hProf->GetXaxis()->SetRangeUser(0.,10.);
+				  hProf->Draw("LP");
+			    } else {
+				  hProf->Draw("LPsame");
+			    }
+			  }
+              
+		    }
+		  break;
+		default:
+            hname=Form("cal_SmT%01d_sm%03d_rpc%03d_%s",iType[iSt],iSm,iRp,cOpt.Data());
+            h=(TH2 *)gROOT->FindObjectAny(hname);
+            if (h!=NULL) {
+				if (iOpt == 2 || iOpt==2 ) {
+					gPad->SetLogz();
+				}
+				h->Draw("colz");
+				h->ProfileX()->Draw("same");
+				iDet++;
+				dAvMean += h->ProfileX()->GetMean(2);
+				dAvRMS  += h->ProfileX()->GetRMS(2);
+				cout << "TrackQA " << cOpt.Data() <<" for TSR " << iType[iSt] << iSm << iRp << ": Off "<< h->ProfileX()->GetMean(2) 
+				<< ", RMS "<< h->ProfileX()->GetRMS(2) << endl;
+			}
+		}	  
+	}
+  }
+ }
+ dAvMean /= (Double_t)iDet;
+ dAvRMS  /= (Double_t)iDet;
+ cout << "TrackQA " << cOpt.Data() <<": AvOff "<< dAvMean << ", AvRMS "<< dAvRMS << endl;
+ dAvMean=TMath::Abs(dAvMean);
+ gROOT->ProcessLine(Form(".! echo %d > %sAvOff.res", (Int_t)(dAvMean*1.E4),cOpt.Data()));
+ gROOT->ProcessLine(Form(".! echo %d > %sAvRMS.res", (Int_t)(dAvRMS*1.E4),cOpt.Data()));
+
+ can->SaveAs(Form("pl_all_Track_%s.pdf",cOpt.Data()));
+
+} 
diff --git a/macro/beamtime/hd2020/pl_cmp_CluRate.C b/macro/beamtime/hd2020/pl_cmp_CluRate.C
new file mode 100755
index 0000000000000000000000000000000000000000..572eb587553bb3c80707e973a01e95985832dff8
--- /dev/null
+++ b/macro/beamtime/hd2020/pl_cmp_CluRate.C
@@ -0,0 +1,123 @@
+void pl_cmp_CluRate(Int_t iNSt=3, Long_t iSet=900032500, Int_t iOpt=0, Double_t Tstart=0., Double_t Tend=10., Int_t iMode=1)
+{
+  //  TCanvas *can = new TCanvas("can22","can22");
+  //  can->Divide(2,2); 
+  //  TCanvas *can = new TCanvas("can","can",48,55,700,900);
+  TCanvas *can = new TCanvas("can","can",48,56,900,900);
+  //  can->Divide(2,2,0,0); 
+  can->Divide(1,3,0.01,0.01);
+  Float_t lsize=0.09;
+
+ gPad->SetFillColor(0);
+ gStyle->SetPalette(1);
+ gStyle->SetLabelSize(lsize);
+
+ gStyle->SetOptStat(kFALSE);
+ //gROOT->cd();
+ //gROOT->SetDirLevel(2);
+ gStyle->SetPadLeftMargin(0.4);
+ gStyle->SetTitleOffset(1.0, "x");
+ gStyle->SetTitleOffset(0.8, "y");
+ gStyle->SetTitleFontSize(0.08);
+ 
+
+ TH1 *h[iNSt];
+ TH1 *hrat1[iNSt];
+ TH1 *hrat2[iNSt];
+
+ TH2 *h2;
+ Int_t   iType[6]={0,9,5,6,7,8};
+ const Int_t  iSmNum[6]={5,1,1,1,1,1};
+ const Int_t iRpcNum[6]={5,2,1,2,1,8};
+ Int_t iId_full[iNSt];
+
+ Int_t iCanv=0;
+ Int_t iCol=0;
+ // if (h!=NULL) h->Delete();
+ TLegend *leg = new TLegend(0.78,0.6,0.85,0.9);
+ leg->SetTextSize(0.05);
+
+ Long_t iSet_tmp=iSet;
+ 
+ can->cd(1);
+ for(Int_t iSt=0; iSt<iNSt; iSt++){
+    gROOT->cd();
+    Int_t iId=iSet_tmp%1000;
+    iId_full[iSt]=iId;
+    iSet_tmp=(iSet_tmp-iId)/1000;
+    Int_t iRp=iId%10;
+    iId -= iRp;    iId /=10;
+    Int_t iSm=iId%10;
+    iId -= iSm;    iId /=10;
+    iType[iSt]=iId;
+    
+    TString hname="";
+    switch(iOpt){
+    case 0:
+      hname=Form("cl_SmT%01d_sm%03d_rpc%03d_rate",iType[iSt],iSm,iRp);
+      break;
+    case 1:
+      hname=Form("cl_SmT%01d_sm%03d_rpc%03d_rate10s",iType[iSt],iSm,iRp);
+      break;
+    }
+    h[iSt]=(TH1 *)gROOT->FindObjectAny(hname);
+    if (h[iSt]!=NULL) {
+     h[iSt]->GetXaxis()->SetRangeUser(Tstart,Tend);
+     cout << "Draw " << iSt << ": " << hname << endl;
+     leg->AddEntry(h[iSt],Form("%03d",iId_full[iSt]),"p");
+    
+     switch(iMode) {
+     case 0:
+       h[iSt]->Draw("");
+     //h->UseCurrentStyle();
+       gPad->SetLogy();
+       break;
+     case 1:
+       h[iSt]->UseCurrentStyle();  // set current defaults 
+       h[iSt]->SetMarkerStyle(20+iSm);
+       h[iSt]->SetTitle("");
+       iCol=iSt+1;
+       if(iCol == 5) iCol++;
+       h[iSt]->SetLineColor(iCol);
+       h[iSt]->SetLineStyle(1);
+       h[iSt]->SetMarkerColor(iCol);
+       if(iSt==0) h[iSt]->Draw("LPE");
+       else h[iSt]->Draw("LPEsame");
+       break;
+     }
+    }else{cout<<"Histogram "<<hname<<" not existing. "<<endl;}
+   }
+  leg->Draw();
+  
+  can->cd(2);
+  Double_t RatMax=1.5;
+   gPad->SetLogy();
+   gPad->SetGridx();
+   gPad->SetGridy();
+  for(Int_t iSt=1; iSt<iNSt; iSt++){
+    hrat1[iSt]=(TH1 *)h[iSt]->Clone();
+    hrat1[iSt]->Divide(h[iSt],h[0],1.,1.);
+    hrat1[iSt]->SetYTitle(Form("Ratio to %03d",iId_full[0]));
+    hrat1[iSt]->SetMaximum(RatMax);
+    hrat1[iSt]->SetMinimum(1.E-2);
+    if(iSt==1) hrat1[iSt]->Draw("LPE");
+    else hrat1[iSt]->Draw("LPEsame");
+  }
+
+  can->cd(3);
+   gPad->SetLogy();
+   gPad->SetGridx();
+   gPad->SetGridy();
+   for(Int_t iSt=2; iSt<iNSt; iSt++){
+    hrat2[iSt]=(TH1 *)h[iSt]->Clone();
+    hrat2[iSt]->Divide(h[iSt],h[1],1.,1.);
+    hrat2[iSt]->SetYTitle(Form("Ratio to %03d",iId_full[1]));
+    hrat2[iSt]->SetMinimum(1.E-2);
+    gPad->SetLogy();
+    if(iSt==2) hrat2[iSt]->Draw("LPE");
+    else hrat2[iSt]->Draw("LPEsame");
+  }
+  
+ can->SaveAs(Form("pl_cmp_CluRate%ld_%d.pdf",iSet,iOpt));
+
+} 
diff --git a/macro/beamtime/hd2020/rootlogon.C b/macro/beamtime/hd2020/rootlogon.C
new file mode 100644
index 0000000000000000000000000000000000000000..638638603193b656eb482bd2fa68ea40a35a893c
--- /dev/null
+++ b/macro/beamtime/hd2020/rootlogon.C
@@ -0,0 +1,65 @@
+{
+  //}
+  //void rootlogon_nh()
+  //{
+  pTime=new TDatime();
+  cout << " Executing rootlogon.C (nh) at " <<pTime->GetDate() << ", "<<pTime->GetTime() <<endl;
+  gStyle->SetOptStat(111);
+  gStyle->SetLineWidth(2.);
+  gStyle->SetFrameLineWidth(2.);
+  gStyle->SetTitleOffset(1.01, "x");
+  gStyle->SetTitleOffset(0.9, "y");
+  gStyle->SetTitleSize(0.08, "x"); // axis labels 
+  gStyle->SetTitleSize(0.08, "y");
+  gStyle->SetHistLineWidth(2.);
+  gStyle->SetPadRightMargin(0.15);
+  gStyle->SetPadLeftMargin(0.2);
+  gStyle->SetPadBottomMargin(0.2);
+  gStyle->SetTitleFontSize(0.07); // histogram title
+  gStyle->SetLabelSize(0.07, "x"); //numbers
+  gStyle->SetLabelSize(0.07, "y"); //numbers
+  gStyle->SetLabelSize(0.05, "z"); //numbers
+  gStyle->SetLabelOffset(0.02, "x"); //numbers
+  gStyle->SetLabelOffset(0.02, "y"); //numbers
+  gStyle->SetLabelOffset(0.02, "z"); //numbers
+  gStyle->SetTextSize(0.3);
+  gStyle->SetNdivisions(505, "x");
+  gStyle->SetNdivisions(505, "y");
+  gStyle->SetNdivisions(505, "z");
+  gStyle->SetTickLength(0.06, "x");
+  gStyle->SetTickLength(0.03, "y");
+  gStyle->SetTickLength(0.06, "z");
+  gStyle->SetPalette(1,0);
+  //  gStyle->SetOptDate(23);
+  //  gStyle->SetDateX(0.5);
+  // gStyle->SetDateY(0.5);
+
+  gStyle->SetLineScalePS(1.0);
+
+  gSystem->AddIncludePath(Form("-I%s/include",gSystem->Getenv("FAIRROOTPATH")));
+
+  gSystem->AddIncludePath(Form("-I%s/roclight",gSystem->Getenv("VMCWORKDIR")));
+  gSystem->AddIncludePath(Form("-I%s/data",gSystem->Getenv("VMCWORKDIR")));
+  gSystem->AddIncludePath(Form("-I%s/data/tof",gSystem->Getenv("VMCWORKDIR")));
+  gSystem->AddIncludePath(Form("-I%s/base ",gSystem->Getenv("VMCWORKDIR")));
+  gSystem->AddIncludePath(Form("-I%s/beamtime/tof",gSystem->Getenv("VMCWORKDIR")));
+  gSystem->AddIncludePath(Form("-I%s/beamtime/tof/hadaq",gSystem->Getenv("VMCWORKDIR")));
+  gSystem->AddIncludePath(Form("-I%s/beamtime/tof/tdc",gSystem->Getenv("VMCWORKDIR")));
+  gSystem->AddIncludePath(Form("-I%s/beamtime/tof/tdc/v1290",gSystem->Getenv("VMCWORKDIR")));
+  gSystem->AddIncludePath(Form("-I%s/beamtime/tof/tdc/vftx",gSystem->Getenv("VMCWORKDIR")));
+  gSystem->AddIncludePath(Form("-I%s/beamtime/tof/tdc/trb",gSystem->Getenv("VMCWORKDIR")));
+  gSystem->AddIncludePath(Form("-I%s/beamtime/tof/tdc/get4",gSystem->Getenv("VMCWORKDIR")));
+  gSystem->AddIncludePath(Form("-I%s/beamtime/tof/scalers",gSystem->Getenv("VMCWORKDIR")));
+  gSystem->AddIncludePath(Form("-I%s/beamtime/tof/scalers/triglog",gSystem->Getenv("VMCWORKDIR")));
+  gSystem->AddIncludePath(Form("-I%s/beamtime/tof/scalers/scalormu",gSystem->Getenv("VMCWORKDIR")));
+  gSystem->AddIncludePath(Form("-I%s/beamtime/tof/scalers/scal2014",gSystem->Getenv("VMCWORKDIR")));
+  gSystem->AddIncludePath(Form("-I%s/beamtime/tof/scalers/orgen",gSystem->Getenv("VMCWORKDIR")));
+  gSystem->AddIncludePath(Form("-I%s/beamtime/tof/unpMoni",gSystem->Getenv("VMCWORKDIR")));
+  gSystem->AddIncludePath(Form("-I%s/beamtime/tof/calib",gSystem->Getenv("VMCWORKDIR")));
+  gSystem->AddIncludePath(Form("-I%s/beamtime/tof/calib/tdc",gSystem->Getenv("VMCWORKDIR")));
+  gSystem->AddIncludePath(Form("-I%s/beamtime/tof/calib/scaler",gSystem->Getenv("VMCWORKDIR")));
+  gSystem->AddIncludePath(Form("-I%s/beamtime/tof/mapping",gSystem->Getenv("VMCWORKDIR")));
+  gSystem->AddIncludePath(Form("-I%s/beamtime/tof/output",gSystem->Getenv("VMCWORKDIR")));
+  gSystem->AddIncludePath(Form("-I%s/beamtime/tof/display",gSystem->Getenv("VMCWORKDIR")));
+
+}