diff --git a/doc/musrSim.tex b/doc/musrSim.tex
index a21af0b..9b34337 100644
--- a/doc/musrSim.tex
+++ b/doc/musrSim.tex
@@ -845,7 +845,7 @@ The list of variables that can be stored in the Root tree:
 \item{\bf muIniMomX, muIniMomY, muIniMomZ}  (Double\_t) -- initial momentum of the muon when it was generated (in MeV/c).
 \item{\bf muIniPolX, muIniPolY, muIniPolZ}  (Double\_t) -- initial polarisation of the muon when it was generated.
 \item{\bf muDecayDetID} (Int\_t) -- ID number of the detector in which the muon stopped and decayed.
-\item{\bf muDecayTime} (Double\_t) -- the time at which the muon stopped and decayed (in $\mu$s).
+\item{\bf muDecayTime} (Double\_t) -- the time at which the muon decayed (in $\mu$s).
 \item{\bf muDecayPosX, muDecayPosY, muDecayPosZ}  (Double\_t) -- the position where the muon stopped and decayed (in mm).
 \item{\bf muDecayPolX, muDecayPolY, muDecayPolZ}  (Double\_t) -- polarisation of the muon when it stopped and decayed.
 \item{\bf muTargetTime}  (Double\_t) -- time at which the muon entered the volume whose name starts by ``target'' -- usually the sample (in $\mu$s).
diff --git a/include/musrDetectorConstruction.hh b/include/musrDetectorConstruction.hh
index e235e30..0c3cac9 100644
--- a/include/musrDetectorConstruction.hh
+++ b/include/musrDetectorConstruction.hh
@@ -29,6 +29,7 @@
 #include "G4ThreeVector.hh"
 #include "G4RotationMatrix.hh"
 #include "G4FieldManager.hh"
+#include "G4OpticalSurface.hh"
 #include <map>
 
 class G4Tubs;
@@ -47,7 +48,7 @@ class musrDetectorConstruction : public G4VUserDetectorConstruction
 {
 public:
   
-  musrDetectorConstruction();
+  musrDetectorConstruction(G4String steeringFileName);
   ~musrDetectorConstruction();
 
 public:
@@ -61,6 +62,7 @@ public:
   void ReportProblemInStearingFile(char* myString);
   G4Material* CharToMaterial(char myString[100]);
   G4LogicalVolume* FindLogicalVolume(G4String LogicalVolumeName);
+  G4VPhysicalVolume* FindPhysicalVolume(G4String PhysicalVolumeName);
   void SetColourOfLogicalVolume(G4LogicalVolume* pLogVol,char* colour);
 
 private:
@@ -73,6 +75,10 @@ private:
 
   std::map<std::string,G4RotationMatrix*> pointerToRotationMatrix;
   std::map<std::string,G4FieldManager*> pointerToField;
+
+  std::map<std::string,G4MaterialPropertiesTable*> materialPropertiesTableMap;
+  std::map<std::string,G4MaterialPropertiesTable*>::iterator  itMPT;
+
   
 private:
   void DefineMaterials(); 
diff --git a/include/musrRootOutput.hh b/include/musrRootOutput.hh
index aabe36c..af7a7cc 100644
--- a/include/musrRootOutput.hh
+++ b/include/musrRootOutput.hh
@@ -79,6 +79,9 @@ class musrRootOutput  {
 			     G4double ekVertex, G4double xVertex, G4double yVertex, G4double zVertex, 
 			     G4int idVolVertex, G4int idProcVertex, G4int idTrackVertex, G4int particleID) ;
 
+    void SetOPSAinfo    (G4int nDetectors, G4int ID, G4int nPhot, G4double timeFirst, G4double timeA, 
+			 G4double timeB, G4double timeC, G4double timeD, G4double timeE, G4double timeLast);
+
     void SetSaveDetectorInfo (G4int ID, G4int particleID, G4double ke, G4double x, G4double y, G4double z, G4double time, 
 			      G4double px, G4double py, G4double pz, G4double polx, G4double poly, G4double polz) ;
 
@@ -208,7 +211,15 @@ class musrRootOutput  {
     static G4bool store_fieldIntegralBz1;
     static G4bool store_fieldIntegralBz2;
     static G4bool store_fieldIntegralBz3;
-
+    static G4bool store_odet_ID;
+    static G4bool store_odet_nPhot;
+    static G4bool store_odet_timeFirst;
+    static G4bool store_odet_timeA;
+    static G4bool store_odet_timeB;
+    static G4bool store_odet_timeC;
+    static G4bool store_odet_timeD;
+    static G4bool store_odet_timeE;
+    static G4bool store_odet_timeLast;
 
     static G4int oldEventNumberInG4EqEMFieldWithSpinFunction;
 
@@ -299,6 +310,21 @@ class musrRootOutput  {
     G4int     det_VvvTrackID[det_nMax];
     G4int     det_VvvParticleID[det_nMax];
 
+  public:
+    static const Int_t odet_nMax=det_nMax;    
+  
+  private:
+    G4int odet_n;
+    G4int odet_ID[odet_nMax];
+    G4int odet_nPhot[odet_nMax];
+    G4double odet_timeFirst[odet_nMax];
+    G4double odet_timeA[odet_nMax];
+    G4double odet_timeB[odet_nMax];
+    G4double odet_timeC[odet_nMax];
+    G4double odet_timeD[odet_nMax];
+    G4double odet_timeE[odet_nMax];
+    G4double odet_timeLast[odet_nMax];
+
   public:
     static const Int_t save_nMax=1000;
 
diff --git a/include/musrScintSD.hh b/include/musrScintSD.hh
index 591e405..fcecd38 100644
--- a/include/musrScintSD.hh
+++ b/include/musrScintSD.hh
@@ -31,26 +31,87 @@
 class G4Step;
 class G4HCofThisEvent;
 
+class signalInfo {
+  public:
+    signalInfo(G4int id, G4int nP, G4double tFirst, 
+  	       G4double tA, G4double tB, G4double tC, G4double tD, G4double tE, G4double tLast) {
+               detID=id; nPhot=nP; timeFirst=tFirst;
+               timeA=tA; timeB=tB; timeC=tC; timeD=tD; timeE=tE; timeLast=tLast;}
+    ~signalInfo() {}
+    void transferDataToRoot(musrRootOutput* myRootOut, G4int nn) {
+                myRootOut->SetOPSAinfo(nn,detID,nPhot,timeFirst,timeA,timeB,timeC,timeD,timeE,timeLast);
+    }
+
+  private:
+    G4int detID;
+    G4int nPhot;
+    G4int nPhot_abs;
+    G4int nPhot_refl;
+    G4int nPhot_other;
+    G4double timeFirst;
+    G4double timeA;
+    G4double timeB;
+    G4double timeC;
+    G4double timeD;
+    G4double timeE;
+    G4double timeLast;
+};
+
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 class musrScintSD : public G4VSensitiveDetector
 {
   public:
-      musrScintSD(G4String);
-     ~musrScintSD();
+    static musrScintSD* GetInstance();
+    musrScintSD(G4String);
+    ~musrScintSD();
 
-      void Initialize(G4HCofThisEvent*);
-      G4bool ProcessHits(G4Step*, G4TouchableHistory*);
-      void EndOfEvent(G4HCofThisEvent*);
+    void Initialize(G4HCofThisEvent*);
+    G4bool ProcessHits(G4Step*, G4TouchableHistory*);
+    void EndOfEvent(G4HCofThisEvent*);
+
+    // Optical Photon Signal Analysis (OPSA)
+    void Set_OPSA_minNrOfDetectedPhotons(G4int val) {OPSA_minNrOfDetectedPhotons=val;}
+    void Set_OPSA_SignalSeparationTime(G4double val) {OPSA_signalSeparationTime=val;}
+    void Set_OPSA_frac(G4double a, G4double b, G4double c, G4double d, G4double e) 
+             {OPSA_fracA=a; OPSA_fracB=b; OPSA_fracC=c; OPSA_fracD=d; OPSA_fracE=e;}
+    void AddEventIDToMultimapOfEventIDsForOPSAhistos (G4int ev_ID, G4int detector_ID)
+          {bool_multimapOfEventIDsForOPSAhistosEXISTS=true; multimapOfEventIDsForOPSAhistos.insert(std::pair<G4int,G4int>(ev_ID,detector_ID));}
+    void SetOPSAhistoBinning(Int_t nBins, Double_t min, Double_t max) {OPSAhistoNbin=nBins; OPSAhistoMin=min; OPSAhistoMax=max;}
+    void ProcessOpticalPhoton(G4Step*);
+    void EndOfEvent_OptiacalPhotons();
 
   private:
+      static musrScintSD* pointer;
       musrScintHitsCollection* scintCollection;
       G4bool   myStoreOnlyEventsWithHits;
       G4int    myStoreOnlyEventsWithHitInDetID;
       G4double mySignalSeparationTime;
       G4bool   myStoreOnlyTheFirstTimeHit;
       G4bool   boolIsVvvInfoRequested;
+      musrRootOutput* myRootOutput;
 
+      // for optical photon counting
+      typedef std::multimap<G4double,Int_t> optHitDetectorMapType;
+      typedef std::map<Int_t,optHitDetectorMapType*> optHitMapType;
+      optHitMapType optHitMap;
+      // for optical photon signal analysis (OPSA)
+      G4int    OPSA_minNrOfDetectedPhotons;
+      G4double OPSA_signalSeparationTime;
+      G4double OPSA_fracA;
+      G4double OPSA_fracB;
+      G4double OPSA_fracC;
+      G4double OPSA_fracD;
+      G4double OPSA_fracE;
+      typedef std::multimap<G4int,signalInfo*>  OPSA_signal_MapType;
+      OPSA_signal_MapType OPSA_signal_Map;
+  
+      G4bool bool_multimapOfEventIDsForOPSAhistosEXISTS;
+      typedef std::multimap<G4int,G4int> multimapOfEventIDsForOPSAhistos_Type;
+      multimapOfEventIDsForOPSAhistos_Type multimapOfEventIDsForOPSAhistos;
+      TH1D* OPSAhisto;
+      Int_t OPSAhistoNbin;
+      Double_t OPSAhistoMin, OPSAhistoMax;
 };
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/include/musrSteppingAction.hh b/include/musrSteppingAction.hh
index 4e79e01..9837d23 100644
--- a/include/musrSteppingAction.hh
+++ b/include/musrSteppingAction.hh
@@ -23,7 +23,6 @@
 
 #ifndef musrSteppingAction_h
 #define musrSteppingAction_h 1
-
 #include "G4UserSteppingAction.hh"
 #include "G4ProcessManager.hh"
 #include "globals.hh"
@@ -32,6 +31,7 @@
 
 class G4Timer;
 
+
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 class musrSteppingAction : public G4UserSteppingAction
@@ -55,10 +55,10 @@ class musrSteppingAction : public G4UserSteppingAction
     void   SetCalculationOfFieldIntegralRequested(G4bool decision) {boolCalculateFieldIntegral = decision;}
 
   private:
-    musrRootOutput* myRootOutput;
     G4Timer* timer;
     time_t realTimeWhenThisEventStarted;
     static musrSteppingAction* pointer;
+    musrRootOutput* myRootOutput;
 
     G4bool muAlreadyWasInTargetInThisEvent;
     G4bool muAlreadyWasInM0InThisEvent;
diff --git a/musrSim.cc b/musrSim.cc
index d884660..38fc621 100644
--- a/musrSim.cc
+++ b/musrSim.cc
@@ -3,6 +3,7 @@
 #include "musrPrimaryGeneratorAction.hh"
 #include "musrRunAction.hh"
 #include "musrEventAction.hh"
+//#include "StackingAction.hh"
 #include "musrSteppingAction.hh"
 #include "musrSteppingVerbose.hh"
 
@@ -11,6 +12,10 @@
 #include "G4UIterminal.hh"
 #include "G4UItcsh.hh"
 
+//#include <TApplication.h>
+//#include <TSystem.h>
+
+
 // The following two lines are needed to cope with the problem of 
 // "Error in <TPluginManager::FindHandler>: Cannot find plugin handler for TVirtualStreamerInfo! 
 //               Does $ROOTSYS/etc/plugins/TVirtualStreamerInfo exist?"
@@ -53,6 +58,8 @@ int main(int argc,char** argv) {
   // Create class "musrErrorMessage"
   musrErrorMessage* myErrorMessage = new musrErrorMessage();
 
+  //  TApplication* myapp=new TApplication("myapp",0,0);
+  
   // Create Root class for storing the output of the Geant simulation
   musrRootOutput* myRootOutput = new musrRootOutput();
 
@@ -68,13 +75,14 @@ int main(int argc,char** argv) {
 
 
   // UserInitialization classes (mandatory)
-  musrDetectorConstruction* musrdetector = new musrDetectorConstruction;
+  //  musrDetectorConstruction* musrdetector = new musrDetectorConstruction;
   if (argc>1) {
     G4int myRunNr=atoi(argv[1]);       // Get the run number from the name of the 
                                        // parameter file, if it starts with a number.
     if (myRunNr>0)  {runManager->SetRunIDCounter(myRunNr);}
-    musrdetector->SetInputParameterFileName(argv[1]);
+    //    musrdetector->SetInputParameterFileName(argv[1]);
   }
+  musrDetectorConstruction* musrdetector = new musrDetectorConstruction(steeringFileName);
   runManager->SetUserInitialization(musrdetector);
   runManager->SetUserInitialization(new musrPhysicsList);
   
@@ -89,6 +97,7 @@ int main(int argc,char** argv) {
   runManager->SetUserAction(new musrPrimaryGeneratorAction(musrdetector));
   runManager->SetUserAction(new musrRunAction);  
   runManager->SetUserAction(new musrEventAction);
+  //  runManager->SetUserAction(new StackingAction);
   runManager->SetUserAction(new musrSteppingAction);
 
   //Initialize G4 kernel
@@ -133,6 +142,8 @@ int main(int argc,char** argv) {
     }
   }
 
+  //  myapp->Run(kTRUE);
+
 #ifdef G4VIS_USE
   delete visManager;
 #endif
diff --git a/musrSimAna/musrAnalysis.cxx b/musrSimAna/musrAnalysis.cxx
index 76dae6a..9970a10 100644
--- a/musrSimAna/musrAnalysis.cxx
+++ b/musrSimAna/musrAnalysis.cxx
@@ -75,12 +75,12 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
     exit(1);
   }
   std::cout << "Configuration file \"" << charV1190FileName << "\" was opened."<<std::endl;
-  char  line[1001];
+  char  line[10001];
 
   // Write out the configuration file into the output file:
   std::cout << "\n\n....oooOO0OOooo........oooOO0OOooo......Configuration file used for this run....oooOO0OOooo........oooOO0OOooo......"<<std::endl; 
   while (!feof(fSteeringFile)) {
-    fgets(line,1000,fSteeringFile);
+    fgets(line,10000,fSteeringFile);
     //    if ((line[0]!='#')&&(line[0]!='\n')&&(line[0]!='\r'))  // TS: Do not print comments
     std::cout << line;
   }
@@ -120,7 +120,7 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
   //    }
   //  }
   while (!feof(fSteeringFile)) {
-    fgets(line,1000,fSteeringFile);
+    fgets(line,10000,fSteeringFile);
     if ((line[0]!='#')&&(line[0]!='\n')&&(line[0]!='\r')&&(line[0]!='$')&&(line[0]!='!')) {
       char tmpString0[200]="Unset";
       char tmpString1[200]="Unset";
@@ -194,6 +194,18 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
 	sscanf(&line[strlen("MUSRMODE=")],"%s",&tmpString1);
 	musrMode = tmpString1[0];
       }
+      else if (strcmp(tmpString0,"CLONECHANNEL")==0) {
+	int ichannel_orig_tmp, ichannel_new_tmp;
+	sscanf(&line[0],"%*s %d %d",&ichannel_orig_tmp,&ichannel_new_tmp);
+	bool_clonedChannelsMultimap_NotEmpty = true;
+	clonedChannelsMultimap.insert(std::pair<int,int>(ichannel_orig_tmp,ichannel_new_tmp));
+      }
+      else if (strcmp(tmpString0,"DEBUGEVENT")==0) {
+	int ieventToDebug_tmp, iLevelToDebug_tmp;
+	sscanf(&line[0],"%*s %d %d",&ieventToDebug_tmp,&iLevelToDebug_tmp);
+	bool_debugingRequired=true;
+	debugEventMap.insert(std::pair<int,int>(ieventToDebug_tmp,iLevelToDebug_tmp));
+      }
       else if (strncmp(tmpString0,"musrTH",strlen("musrTH"))==0) {
 	// Definition of the histograms - either musrTH1D or musrTH2D
 	int beginningOfHistoTitle=0, endOfHistoTitle =0;
@@ -355,6 +367,18 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
 	else if (strcmp(conditionNameTMP,"pileupEventCandidate")==0)    conditionMap[iConditionTMP]=&pileupEventCandidate;
 	else if (strcmp(conditionNameTMP,"pileupEvent")==0)             conditionMap[iConditionTMP]=&pileupEvent;
 	else if (strcmp(conditionNameTMP,"goodEvent_det_AND_muonDecayedInSample_gen")==0) conditionMap[iConditionTMP]=&goodEvent_det_AND_muonDecayedInSample_gen;
+	else if (strcmp(conditionNameTMP,"goodEvent_F_det")==0)         conditionMap[iConditionTMP]=&goodEvent_F_det;
+	else if (strcmp(conditionNameTMP,"goodEvent_B_det")==0)         conditionMap[iConditionTMP]=&goodEvent_B_det;
+	else if (strcmp(conditionNameTMP,"goodEvent_U_det")==0)         conditionMap[iConditionTMP]=&goodEvent_U_det;
+	else if (strcmp(conditionNameTMP,"goodEvent_D_det")==0)         conditionMap[iConditionTMP]=&goodEvent_D_det;
+	else if (strcmp(conditionNameTMP,"goodEvent_L_det")==0)         conditionMap[iConditionTMP]=&goodEvent_L_det;
+	else if (strcmp(conditionNameTMP,"goodEvent_R_det")==0)         conditionMap[iConditionTMP]=&goodEvent_R_det;
+	else if (strcmp(conditionNameTMP,"goodEvent_F_det_AND_pileupEvent")==0) conditionMap[iConditionTMP]=&goodEvent_F_det_AND_pileupEvent;
+	else if (strcmp(conditionNameTMP,"goodEvent_B_det_AND_pileupEvent")==0) conditionMap[iConditionTMP]=&goodEvent_B_det_AND_pileupEvent;
+	else if (strcmp(conditionNameTMP,"goodEvent_U_det_AND_pileupEvent")==0) conditionMap[iConditionTMP]=&goodEvent_U_det_AND_pileupEvent;
+	else if (strcmp(conditionNameTMP,"goodEvent_D_det_AND_pileupEvent")==0) conditionMap[iConditionTMP]=&goodEvent_D_det_AND_pileupEvent;
+	else if (strcmp(conditionNameTMP,"goodEvent_L_det_AND_pileupEvent")==0) conditionMap[iConditionTMP]=&goodEvent_L_det_AND_pileupEvent;
+	else if (strcmp(conditionNameTMP,"goodEvent_R_det_AND_pileupEvent")==0) conditionMap[iConditionTMP]=&goodEvent_R_det_AND_pileupEvent;
 	else {
 	  std::cout<<" !!!  ERROR: Condition of the name \""<<conditionNameTMP<<"\" not predefined ==> Add it in the musrAnalysis.cxx   S T O P !!!"<<std::endl;
 	  exit(1);
@@ -404,6 +428,47 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
 	  nscan = sscanf(pch,"%d %g",&N1,&PHASE_SHIFT);
 	} while (nscan==2);
       }
+      else if (strcmp(tmpString0,"counterGrouping")==0) {
+	int nscan; int N1; char NAME[100];
+	char *pch = line + strlen("counterGrouping");
+	char counterGroupName[100];
+	nscan = sscanf(pch,"%s",counterGroupName);
+	char* pch1 = strstr(pch,counterGroupName)+strlen(counterGroupName);
+	pch = pch1;
+	do {
+	  nscan = sscanf(pch,"%d",&N1);
+	  if      (strcmp(counterGroupName,"F")==0) F_posCounterList.push_back(N1);
+	  else if (strcmp(counterGroupName,"B")==0) B_posCounterList.push_back(N1);
+	  else if (strcmp(counterGroupName,"U")==0) U_posCounterList.push_back(N1);
+	  else if (strcmp(counterGroupName,"D")==0) D_posCounterList.push_back(N1);
+	  else if (strcmp(counterGroupName,"L")==0) L_posCounterList.push_back(N1);
+	  else if (strcmp(counterGroupName,"R")==0) R_posCounterList.push_back(N1);
+	  else {
+	    std::cout<<"\n\n UNKNOWN COUNTER GROUP REQUIRED !!!    =====>    S T O P    F O R C E D"<<std::endl;
+	    std::cout<<line<<std::endl;
+	    exit(1);
+	  }
+	  std::cout<<"counterGroupName="<<counterGroupName<<"   N1="<<N1<<std::endl;
+	  nscan = sscanf(pch,"%s",NAME);
+	  char* pch2 = strstr(pch ,NAME)+strlen(NAME);
+	  pch=pch2;
+	  nscan = sscanf(pch,"%d",&N1);
+	} while (nscan==1);
+      }
+      else if  (strcmp(tmpString0,"artificiallyChangeMuDecayTime")==0) {
+	float min, max, mmmin, mmmax;
+	sscanf(&line[0],"%*s %g %g %g %g %g %g %g",&min,&max,&mmmin,&mmmax);
+	bool_muDecayTimeTransformation = true;
+	muDecayTime_t_min              = min;
+        muDecayTime_t_max              = max;
+	muDecayTime_Transformation_min = mmmin;
+        muDecayTime_Transformation_max = mmmax;
+	if ((muDecayTime_t_max <= muDecayTime_t_min) || (muDecayTime_Transformation_max <= muDecayTime_Transformation_min)) {
+	  std::cout<<" ERROR!  musrAnalysis: error when setting the \"artificiallyChangeMuDecayTime\" parameters!  Min > Max !!!"<<std::endl;
+	  std::cout<<"    ==> S T O P "<<std::endl;
+	  exit(1);
+	}
+      }
       else if  (strcmp(tmpString0,"fit")==0) {
 	char histoName[100]; char functionName[100]; float xMin; float xMax; float p0, p1, p2, p3, p4, p5, p6;
 	sscanf(&line[0],"%*s %s %s %g %g %g %g %g %g %g",histoName,functionName,&xMin,&xMax,&p0,&p1,&p2,&p3,&p4,&p5,&p6);
@@ -453,7 +518,8 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
 	  funct -> SetParameter(1,p1);
 	}
 	else if (strcmp(functionName,"rotFrameTime20")==0) {
-	  funct = new TF1("rotFrameTime20","[2]*exp(-x/2.19703)*cos(x*[0]+[1]) ");
+	  //	  funct = new TF1("rotFrameTime20","[2]*exp(-x/2.19703)*cos(x*[0]+[1]) ");
+	  funct = new TF1("rotFrameTime20","[2]*cos(x*[0]+[1]) ");
 	  funct -> SetParameter(0,p0);
 	  funct -> SetParameter(1,p1);
 	  funct -> SetParameter(2,p2);
@@ -574,22 +640,8 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
   pileupWindowBinMax = Long64_t(pileupWindowMax/tdcresolution);
   dataWindowBinMin   = Long64_t(  dataWindowMin/tdcresolution);
   dataWindowBinMax   = Long64_t(  dataWindowMax/tdcresolution);
+
   for (counterMapType::const_iterator it = allCounterMap.begin(); it!=allCounterMap.end(); ++it) {
-    char DetectorType = (it->second)->GetCounterType();
-    if      (DetectorType=='M') {
-      (it->second)->SetMaxCoincidenceTimeWindow(pileupWindowBinMin);
-      (it->second)->SetCoincidenceTimeWindow(pileupWindowBinMin,pileupWindowBinMax);
-      (it->second)->SetCoincidenceTimeWindowOfAllCoincidenceDetectors(-mcoincwin,-mcoincwin,mcoincwin);
-    }
-    else if (DetectorType=='P') {
-      (it->second)->SetMaxCoincidenceTimeWindow(dataWindowBinMin);
-      (it->second)->SetCoincidenceTimeWindow(dataWindowBinMin,dataWindowBinMax);
-      (it->second)->SetCoincidenceTimeWindowOfAllCoincidenceDetectors(-pcoincwin,-pcoincwin,pcoincwin);
-    }
-    else if (DetectorType=='V') {
-      (it->second)->SetMaxCoincidenceTimeWindow(-vcoincwin);
-      (it->second)->SetCoincidenceTimeWindow(-vcoincwin,vcoincwin);
-    }
     int iChanNr = (it->second)->GetCounterNr();
     for (std::multimap<int,int>::iterator itCoinc = tmpCoincidenceMultimap.begin(); itCoinc != tmpCoincidenceMultimap.end();  ++itCoinc)  {
       if ((*itCoinc).first == iChanNr) {
@@ -606,6 +658,26 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
     }
   }
 
+  for (counterMapType::const_iterator it = allCounterMap.begin(); it!=allCounterMap.end(); ++it) {
+    char DetectorType = (it->second)->GetCounterType();
+    if      (DetectorType=='M') {
+      (it->second)->SetMaxCoincidenceTimeWindow(pileupWindowBinMin);
+      (it->second)->SetCoincidenceTimeWindow_M(pileupWindowBinMin,pileupWindowBinMax);
+      (it->second)->SetCoincidenceTimeWindowOfAllCoincidenceDetectors('M',-mcoincwin+pileupWindowBinMin,-mcoincwin,mcoincwin);
+      (it->second)->SetCoincidenceTimeWindowOfAllVetoDetectors(-mcoincwin+pileupWindowBinMin,-vcoincwin,vcoincwin);
+    }
+    else if (DetectorType=='P') {
+      (it->second)->SetMaxCoincidenceTimeWindow(dataWindowBinMin);
+      (it->second)->SetCoincidenceTimeWindow_P(dataWindowBinMin,dataWindowBinMax);
+      (it->second)->SetCoincidenceTimeWindowOfAllCoincidenceDetectors('P',-pcoincwin+dataWindowBinMin,-pcoincwin,pcoincwin);
+      (it->second)->SetCoincidenceTimeWindowOfAllVetoDetectors(-pcoincwin+dataWindowBinMin,-vcoincwin,vcoincwin);
+    }
+    //    else if (DetectorType=='V') {
+    //      (it->second)->SetMaxCoincidenceTimeWindow(-vcoincwin);
+    //      (it->second)->SetCoincidenceTimeWindow(-vcoincwin,vcoincwin);
+    //    }
+  }
+
     //    for (int j=0; j<maxChannels; j++) {
     //    for (multimap<int,int>::iterator itCoinc = tmpCoincidenceMultimap.begin(); itCoinc != tmpCoincidenceMultimap.end();  ++itCoinc)  {
     //      std::cout << "   houby  [" << (*itCoinc).first << ", " << (*itCoinc).second << "]" << std::endl;
@@ -656,7 +728,6 @@ void musrAnalysis::CreateHistograms() {
   //  omega = 850.62*fieldValue;  // value used in Geant ?
   omega = 851.610*fieldValue;       // value from the fits of the data
 
-
   hInfo->Fill(1, fieldValue);
   hInfo->Fill(6, runID);
   hInfo->Fill(7, hGeantParameters->GetBinContent(7));
@@ -731,13 +802,19 @@ void musrAnalysis::AnalyseEvent(Long64_t iiiEntry) {
 
   // Loop over several next event and preprocess them (i.e. fill
   // them into the lists/maps of the class musrCounter).
+  if (bool_debugingRequired) {
+    if (debugEventMap[eventID]>0) {std::cout<<"DEBUGEVENT "<<eventID<<"_________________(before \"PreprocessEvent\"_________"<<std::endl;}
+  }
   while (((iiiEntry>lastPreprocessedEntry)||(((nextUnfilledEventTime-currentTime)<safeTimeWindow))&&(!boolInfinitelyLowMuonRate)) && (lastPreprocessedEntry+1<nentries)) {
     Double_t deltaT = PreprocessEvent(lastPreprocessedEntry+1);
     nextUnfilledEventTime+=deltaT;
   };
-  
   fChain->GetEntry(iiiEntry);  InitialiseEvent();
-  
+  if (bool_debugingRequired) {
+    if (debugEventMap[eventID]>2)  PrintHitsInAllCounters();
+    if (debugEventMap[eventID]>1) {std::cout<<"DEBUGEVENT "<<eventID<<"_________________(after \"PreprocessEvent\"_________"<<std::endl;}
+  }
+
   // Loop over all interesting "moments", which are:
   //    1) any hit in the muon counter
   //    2) any new event (even if nothing was detected in the muon counter)
@@ -751,7 +828,14 @@ void musrAnalysis::AnalyseEvent(Long64_t iiiEntry) {
   //  for (counterMapType::const_iterator it = allCounterMap.begin(); it!=allCounterMap.end(); ++it) {
   //    (*it).second->myPrintThisCounter(eventID);
   //  }
+  if (bool_debugingRequired) {
+    if (debugEventMap[eventID]>1) {std::cout<<"DEBUGEVENT "<<eventID<<"_________________(before \"FillHistograms\"_________"<<std::endl;}
+  }
   FillHistograms(iiiEntry);
+  if (bool_debugingRequired) {
+    if (debugEventMap[eventID]>1) {std::cout<<"DEBUGEVENT "<<eventID<<"_________________(after \"FillHistograms\"_________"<<std::endl;}
+  }
+
 
 
   currentTime+=timeToNextEvent*muonRateFactor;
@@ -790,6 +874,14 @@ void musrAnalysis::RewindAllTimeInfo(Long64_t timeBinsToRewind) {
 
 //================================================================
 
+void musrAnalysis::PrintHitsInAllCounters() {
+  std::cout<<"___________________\n";
+  for (counterMapType::const_iterator it = allCounterMap.begin(); it!=allCounterMap.end(); ++it) {
+    (*it).second->myPrintThisCounter(eventID,0);
+  }
+}
+//================================================================
+
 void musrAnalysis::FillHistograms(Int_t iiiEntry) {
   //  std::cout<<"musrAnalysis::FillHistograms() event="<<eventID<<" , bool="<<generatedInfo<<","<<detectorInfo<<std::endl;
 
@@ -823,7 +915,10 @@ void musrAnalysis::FillHistograms(Int_t iiiEntry) {
     // Check whether there was good hit in the Positron counter
     //    Long64_t dataBinMin = (mCounterHitExistsForThisEventID) ? timeBin0+dataWindowBinMin : timeBinOfThePreviouslyProcessedHit-100000000;
     //    Long64_t dataBinMax = (mCounterHitExistsForThisEventID) ? timeBin0+dataWindowBinMax : timeBinOfThePreviouslyProcessedHit+100000000;
+    pileup_eventID             = -1001;
     pileup_muDecayDetID_double = -1001;
+    pileup_muDecayPosZ         = -1000000000;
+    pileup_muDecayPosR         = -1000000000;
     if (mCounterHitExistsForThisEventID) {
       numberOfGoodMuons++;
       Long64_t dataBinMin = timeBin0+dataWindowBinMin;
@@ -833,7 +928,10 @@ void musrAnalysis::FillHistograms(Int_t iiiEntry) {
       if (pCounterHitExistsForThisEventID && (kEntry>0)&&(posEntry>0)&&(kEntry!=posEntry)) {
 	// This must be a pileup event (positron counter hit comes from the different event than the muon counter hit)
 	fChain->GetEntry(posEntry);
+	pileup_eventID             = eventID;
 	pileup_muDecayDetID_double = muDecayDetID;
+	pileup_muDecayPosZ         = muDecayPosZ;
+	pileup_muDecayPosR         = sqrt(muDecayPosX*muDecayPosX+muDecayPosY*muDecayPosY);
 	//	if (pileup_muDecayDetID_double==-1000) {
 	//	  std::cout<<"DEBUG:  pileup_muDecayDetID_double==-1000, posEntry="<<posEntry<<", eventID="<<eventID<<", idetP_edep="<<idetP_edep<<", idetP="<<idetP<<",  idetP_ID="<<idetP_ID<<std::endl;
 	//	}
@@ -875,7 +973,40 @@ void musrAnalysis::FillHistograms(Int_t iiiEntry) {
     pileupEventCandidate = ((kEntry>0)&&(posEntry>0)&&(kEntry!=posEntry)) ? true:false;
     pileupEvent       = pileupEventCandidate&&goodEvent_det;
     goodEvent_det_AND_muonDecayedInSample_gen = goodEvent_det && muonDecayedInSample_gen;
-    //    if (pileupEvent) std::cout<<"pileupEvent (eventID="<<eventID<<std::endl;
+
+     //    posCounterList_Iterator = find(F_posCounterList.begin(), F_posCounterList.end(), idetP_ID);    
+     //    goodEvent_F_det   = posCounterList_Iterator != F_posCounterList.end()
+    
+    goodEvent_F_det = goodEvent_det && ( (find(F_posCounterList.begin(), F_posCounterList.end(), idetP_ID)) != F_posCounterList.end() );
+    goodEvent_B_det = goodEvent_det && ( (find(B_posCounterList.begin(), B_posCounterList.end(), idetP_ID)) != B_posCounterList.end() );
+    goodEvent_U_det = goodEvent_det && ( (find(U_posCounterList.begin(), U_posCounterList.end(), idetP_ID)) != U_posCounterList.end() );
+    goodEvent_D_det = goodEvent_det && ( (find(D_posCounterList.begin(), D_posCounterList.end(), idetP_ID)) != D_posCounterList.end() );
+    goodEvent_L_det = goodEvent_det && ( (find(L_posCounterList.begin(), L_posCounterList.end(), idetP_ID)) != L_posCounterList.end() );
+    goodEvent_R_det = goodEvent_det && ( (find(R_posCounterList.begin(), R_posCounterList.end(), idetP_ID)) != R_posCounterList.end() );
+    //    std::cout<<"goodEvent_F_det="<<goodEvent_F_det<<std::endl;
+    if (pileupEvent&&goodEvent_F_det) {
+      std::cout<<" DEBUG: Pileup Event: eventID = "<<eventID<<"     pileup_eventID = "<<pileup_eventID<<"     det_time10 = "<<det_time10<<std::endl;
+      //      debugEventMap.insert(std::pair<int,int>(eventID,10));
+    }
+    
+    goodEvent_F_det_AND_pileupEvent = goodEvent_F_det && pileupEvent;
+    goodEvent_B_det_AND_pileupEvent = goodEvent_B_det && pileupEvent;
+    goodEvent_U_det_AND_pileupEvent = goodEvent_U_det && pileupEvent;
+    goodEvent_D_det_AND_pileupEvent = goodEvent_D_det && pileupEvent;
+    goodEvent_L_det_AND_pileupEvent = goodEvent_L_det && pileupEvent;
+    goodEvent_R_det_AND_pileupEvent = goodEvent_R_det && pileupEvent;
+
+    //    if (bool_debugingRequired && muonTriggered_det) {
+    //      std::cout<<"DEBUG: goodEvent_det: eventID="<<eventID<<std::endl;
+    //      if (goodEvent_det) std::cout<<"    ___DETECTED___"<<std::endl;
+    //      MyPrintTree();
+    //      MyPrintConditions();
+    //    }
+    
+    //    if (pileupEvent) {
+    //      std::cout<<"\n NEW: pileupEvent: eventID="<<eventID<<",  kEntry="<<kEntry<<",  posEntry="<<posEntry<< std::endl;
+    //      std::cout<<"det_time10 = "<<det_time10<<std::endl;
+    //    }
 
     // Fill pileup-variables, but only if positron comes from different muon than the trigger signal
     //      if (mCounterHitExistsForThisEventID&&pCounterHitExistsForThisEventID)
@@ -938,25 +1069,110 @@ void musrAnalysis::InitialiseEvent() {
 Double_t musrAnalysis::PreprocessEvent(Long64_t iEn) {
   //  std::cout<<"musrAnalysis::PreprocessEvent()"<<std::endl;
   fChain->GetEntry(iEn); InitialiseEvent();
+
+  // Clone some channels into different one, if requested by user
+  // (This is usefull when e.g. user splits a signal from a veto
+  // and uses it in two different ways - e.g. once for vetoing 
+  // muons, and second (with a different threshold) for validating
+  // a positron candidate.  This is initiated by the 
+  // keyword "CLONECHANNEL" in the *.v1190 file
+  
+
+  if (bool_clonedChannelsMultimap_NotEmpty) {
+    //   std::cout<<"det_n="<<det_n<<std::endl;
+    Int_t det_n_OLD=det_n;
+    for (Int_t i=0; i<det_n_OLD; i++) {
+      //      std::cout<<"   det_ID["<<i<<"]="<<det_ID[i]<<"  edep="<<det_edep[i]<<std::endl;
+      clonedChannelsMultimapType::const_iterator it = clonedChannelsMultimap.find(det_ID[i]);
+      //      std::cout<<"              clonedChannelsMultimap[i]="<<clonedChannelsMultimap[i]<<std::endl;
+      if (it!=clonedChannelsMultimap.end()) {
+	
+	int chNumTMP = it->first;
+	//	std::cout<<"DEBUG:  chNumTMP="<<chNumTMP<<"  ";
+	std::pair<clonedChannelsMultimapType::const_iterator,clonedChannelsMultimapType::const_iterator> ret = clonedChannelsMultimap.equal_range(chNumTMP);
+	for (clonedChannelsMultimapType::const_iterator ittt=ret.first; ittt!=ret.second; ++ittt) {
+	  //	  std::cout << "  ittt->second=" << ittt->second;
+	  int chNumNewTMP = ittt->second;
+	  det_ID[det_n]            = chNumNewTMP;
+          det_edep[det_n]          = det_edep[i];
+          det_edep_el[det_n]       = det_edep_el[i];
+	  det_edep_pos[det_n]      = det_edep_pos[i];
+          det_edep_gam[det_n]      = det_edep_gam[i];
+          det_edep_mup[det_n]      = det_edep_mup[i];
+	  det_nsteps[det_n]        = det_nsteps[i];
+          det_length[det_n]        = det_length[i];
+          det_time_start[det_n]    = det_time_start[i];
+	  det_time_end[det_n]      = det_time_end[i];
+          det_x[det_n]             = det_x[i];
+          det_y[det_n]             = det_y[i];
+          det_z[det_n]             = det_z[i];
+	  det_kine[det_n]          = det_kine[i];
+          det_VrtxKine[det_n]      = det_VrtxKine[i];
+          det_VrtxX[det_n]         = det_VrtxX[i];
+	  det_VrtxY[det_n]         = det_VrtxY[i];
+          det_VrtxZ[det_n]         = det_VrtxZ[i];
+          det_VrtxVolID[det_n]     = det_VrtxVolID[i];
+	  det_VrtxProcID[det_n]    = det_VrtxProcID[i];
+	  det_VrtxTrackID[det_n]   = det_VrtxTrackID[i];
+	  det_VrtxParticleID[det_n]= det_VrtxParticleID[i];
+	  det_VvvKine[det_n]       = det_VvvKine[i];
+          det_VvvX[det_n]          = det_VvvX[i];
+	  det_VvvY[det_n]          = det_VvvY[i];
+	  det_VvvZ[det_n]          = det_VvvZ[i];
+	  det_VvvVolID[det_n]      = det_VvvVolID[i];
+	  det_VvvProcID[det_n]     = det_VvvProcID[i];
+	  det_VvvTrackID[det_n]    = det_VvvTrackID[i];
+	  det_VvvParticleID[det_n] = det_VvvParticleID[i];
+	  det_n++;
+	  
+	}
+      }
+    }
+  }
+
   //  std::cout<<"\n musrAnalysis::PreprocessEvent() Filling event "<<eventID<<std::endl;
   //  MyPrintTree();
-
   Double_t globTime = nextUnfilledEventTime;
   for (Int_t i=0; i<det_n; i++) {
     // //    Int_t detID=det_ID[i];
     std::map<int,musrCounter*>::iterator it;
     it = allCounterMap.find(det_ID[i]);
     if (it==allCounterMap.end()) {
-      // uncomment      std::cout<<"Active detector with det_ID="<<det_ID[i]<<" not found !!!"<<std::endl;
+      // std::cout<<"Active detector with det_ID="<<det_ID[i]<<" not found !!!"<<std::endl;
     }
     else {
       //      Double_t omega = 851.372*fieldNomVal[0];
       Double_t dPhaseShift = (omega==0) ? 0:phaseShiftMap[det_ID[i]]/omega;
       //cDEL      if (det_ID[i]<20) std::cout<<"phaseShift = "<<phaseShiftMap[det_ID[i]]<<"    dPhaseShift="<<dPhaseShift<<"    tdcresolution="<<tdcresolution<<std::endl;
-      Long64_t timeBin  = Long64_t( (globTime+det_time_start[i]            )/tdcresolution );
-      Long64_t timeBin2 = Long64_t( (globTime+det_time_start[i]+dPhaseShift)/tdcresolution );
-      //cDEL      if (det_ID[i]<20)       std::cout<<" timeBin = "<<timeBin<<"   ("<<globTime+det_time_start[i]<<")"<<std::endl;
-      //cDEL      if (det_ID[i]<20)       std::cout<<" timeBin2 = "<<timeBin2<<"   ("<<globTime+det_time_start[i]+dPhaseShift<<")"<<std::endl;
+
+      Double_t t1,t2;
+      if (bool_muDecayTimeTransformation) {
+	//  THIS OPTION WAS SUPOSED TO ALLOW THE USER TO SOMEHOW (IN A TRICKY AND NOT 100% RELIABLE WAY)
+	//  SHIFT THE POSITRON COUNTS TO A RESTRICTED TIME WINDOW, AS IF THE MUON DECAYED IN SOME RESTRICTED
+	//  TIME INTERVAL.  THIS, HOWEVER, DOES NOT WORK, IN THE WAY IT IS IMPLEMENTED HERE, BECAUSE THE
+	//  MUON SPIN ROTATION AT REST IS IGNORED HERE - IT WOULD NEED SOME FURTHER WORK TO DO, AND 
+	//  BECOMES DIFFICULT TO INTERPRET LATER ==> WORK ON THIS WAS STOPPED.
+	Double_t diff = muDecayTime_t_max - muDecayTime_t_min;
+	Double_t ttt1 = det_time_start[i];
+	if ((ttt1 > muDecayTime_Transformation_min) && (ttt1 < muDecayTime_Transformation_max)) {
+	  if (muDecayTime < muDecayTime_t_min) {
+	    ttt1 = ( Long64_t((muDecayTime_t_min-muDecayTime)/diff)+1) * diff + det_time_start[i];
+	  }
+	  else if (muDecayTime > muDecayTime_t_max) {
+	    ttt1 = ( Long64_t((muDecayTime_t_min-muDecayTime)/diff)) * diff + det_time_start[i];
+	  }
+	}
+	t1 = (globTime+ttt1            )/tdcresolution;
+	t2 = (globTime+ttt1+dPhaseShift)/tdcresolution;
+	//	std::cout<<"DEBUG:  det_time_start[i]="<<det_time_start[i]<<"    ttt1="<<ttt1<<"   t1="<<t1<<"   t2="<<t2<<std::endl;
+      }
+      else {
+	t1 = (globTime+det_time_start[i]            )/tdcresolution;
+	t2 = (globTime+det_time_start[i]+dPhaseShift)/tdcresolution;
+      }
+      Long64_t timeBin  = Long64_t(t1);
+      Long64_t timeBin2 = Long64_t(t2);
+
       (*it).second->FillHitInCounter(det_edep[i],timeBin,timeBin2,iEn,eventID,i,det_ID[i]);
     }
   }
@@ -969,24 +1185,37 @@ Double_t musrAnalysis::PreprocessEvent(Long64_t iEn) {
 
 //================================================================
 Bool_t musrAnalysis::PositronCounterHit(Int_t evID, Long64_t dataBinMin, Long64_t dataBinMax, Long64_t& tBin1, Long64_t& tBin2, Int_t& kEntry, Int_t& idetP, Int_t& idetP_ID, Double_t& idetP_edep, Bool_t& doubleHit) {
+
+  if (bool_debugingRequired) {
+    if (debugEventMap[eventID]>2) {std::cout<<"DEBUGEVENT:"<<eventID<<"\"PositronCounterHit\":  pCounterMap.size()="<<pCounterMap.size()<<std::endl;}
+  }
   if (pCounterMap.empty()) return false;
 
-  Bool_t positronHitFound = false;
+  Bool_t positronHitFound  = false;
+  Bool_t goodPositronFound = false;
   if (musrMode=='D') {
     // Loop over all positron counters
     for (counterMapType::const_iterator it = pCounterMap.begin(); it!=pCounterMap.end(); ++it) {
       //      std::cout<<" ===POSITRON==="<< pCounterMap.size()<<std::endl;
       //      Bool_t thereWasHit = (it->second)->GetNextGoodHitInThisEvent(evID,dataBinMin,dataBinMax,'P',tBin1,kEntry,idetP,idetP_ID,idetP_edep,doubleHit);
-      Bool_t thereWasHit = (it->second)->GetNextGoodPositron(evID,dataBinMin,dataBinMax,tBin1,tBin2,kEntry,idetP,idetP_ID,idetP_edep,doubleHit);
+      //      Bool_t thereWasHit = (it->second)->GetNextGoodPositron(evID,dataBinMin,dataBinMax,tBin1,tBin2,kEntry,idetP,idetP_ID,idetP_edep,doubleHit);
+      Int_t positronQuality = (it->second)->GetNextGoodPositron(evID,dataBinMin,dataBinMax,tBin1,tBin2,kEntry,idetP,idetP_ID,idetP_edep,doubleHit);
       //      std::cout<<"000000000000  tBin1="<<tBin1<<"  tBin2="<<tBin2<<std::endl;
       if (doubleHit) {return false;}  // There were two hits in the same positron counter
-      if (thereWasHit) {
-	if (positronHitFound) {doubleHit = true; return false;}   // There were two hits in two different positron counters
+      if (positronQuality>0) {
+	if (positronHitFound) {
+	  if (bool_debugingRequired) {
+	    if (debugEventMap[eventID]>2) {std::cout<<"DEBUGEVENT:"<<eventID<<"\"PositronCounterHit\":  Coincidence with other positron candidate in other counter."<<std::endl;}
+	  }
+	  doubleHit = true; 
+	  return false;
+	}   // There were two hits in two different positron counters
 	positronHitFound = true;
+	if (positronQuality>1) goodPositronFound = true;
       }
     }
   }
-  return positronHitFound;
+  return goodPositronFound;
 }
 
 //================================================================
diff --git a/musrSimAna/musrAnalysis.hh b/musrSimAna/musrAnalysis.hh
index 9fcf924..ef76ecd 100644
--- a/musrSimAna/musrAnalysis.hh
+++ b/musrSimAna/musrAnalysis.hh
@@ -22,6 +22,13 @@
 
 class musrTH;
 
+//#include "musrSimGlobal.hh"
+
+typedef std::map<int,int> debugEventMapType;
+extern debugEventMapType debugEventMap;
+extern Bool_t bool_debugingRequired;
+
+
 class musrAnalysis {
 public :
    TTree          *fChain;   //!pointer to the analyzed TTree or TChain
@@ -199,7 +206,10 @@ public :
    Double_t        gen_time10;
    Double_t        det_time10_MINUS_gen_time10;
    Double_t        det_time20;
+   Double_t        pileup_eventID;
    Double_t        pileup_muDecayDetID_double;
+   Double_t        pileup_muDecayPosZ;
+   Double_t        pileup_muDecayPosR;
 
 
    typedef std::map<int, Bool_t*> conditionMapType;
@@ -216,6 +226,18 @@ public :
    Bool_t          pileupEventCandidate;
    Bool_t          pileupEvent;
    Bool_t          goodEvent_det_AND_muonDecayedInSample_gen;
+   Bool_t          goodEvent_F_det;
+   Bool_t          goodEvent_B_det;
+   Bool_t          goodEvent_U_det;
+   Bool_t          goodEvent_D_det;
+   Bool_t          goodEvent_L_det;
+   Bool_t          goodEvent_R_det;
+   Bool_t          goodEvent_F_det_AND_pileupEvent;
+   Bool_t          goodEvent_B_det_AND_pileupEvent;
+   Bool_t          goodEvent_U_det_AND_pileupEvent;
+   Bool_t          goodEvent_D_det_AND_pileupEvent;
+   Bool_t          goodEvent_L_det_AND_pileupEvent;
+   Bool_t          goodEvent_R_det_AND_pileupEvent;
 
    musrAnalysis(TTree *tree=0);
    virtual ~musrAnalysis();
@@ -234,6 +256,7 @@ public :
    virtual void     RemoveOldHitsFromCounters(Long64_t timeBinLimit);
   //   virtual void     RewindAllTimeInfo(Double_t timeToRewind);
    virtual void     RewindAllTimeInfo(Long64_t timeBinsToRewind);
+   virtual void     PrintHitsInAllCounters();
    virtual void     InitialiseEvent();
    virtual Double_t PreprocessEvent(Long64_t iEn);
    virtual Bool_t   PositronCounterHit(Int_t evID, Long64_t dataBinMin, Long64_t dataBinMax, Long64_t& tBin1, Long64_t& tBin2, Int_t& kEntry, Int_t& idetP, Int_t& idetP_ID, Double_t& idetP_edep, Bool_t& doubleHit);
@@ -247,6 +270,10 @@ public :
    TH1D* hGeantParameters;
    TH1D* hInfo;
 
+  //   typedef std::map<int,int> debugEventMapType;
+  //   debugEventMapType debugEventMap;
+  //   Bool_t bool_debugingRequired;
+
    static const Double_t pi=3.14159265358979324;
 
  private:
@@ -290,6 +317,8 @@ public :
    counterMapType allCounterMap;
    Int_t testIVar1;
    Double_t omega;
+   Bool_t   bool_muDecayTimeTransformation;
+   Double_t muDecayTime_Transformation_min, muDecayTime_Transformation_max, muDecayTime_t_min, muDecayTime_t_max;
 
    static const Double_t microsecond=1.;
    static const Double_t nanosecond=0.001;
@@ -327,6 +356,19 @@ private:
   typedef std::multimap<Int_t,Int_t> humanDecayMultimapType;
   humanDecayMapType humanDecayMap;
   humanDecayMultimapType humanDecayMultimap;
+
+  typedef std::multimap<int,int> clonedChannelsMultimapType;
+  clonedChannelsMultimapType clonedChannelsMultimap;
+  Bool_t bool_clonedChannelsMultimap_NotEmpty;
+
+  // List of group of detectors: F,B,U,D,L,R:
+  std::list <Int_t> F_posCounterList;
+  std::list <Int_t> B_posCounterList;
+  std::list <Int_t> U_posCounterList;
+  std::list <Int_t> D_posCounterList;
+  std::list <Int_t> L_posCounterList;
+  std::list <Int_t> R_posCounterList;
+  //  std::list <Int_t>::iterator posCounterList_Iterator;
 };
 
 #endif
@@ -399,7 +441,10 @@ musrAnalysis::musrAnalysis(TTree *tree)
   variableMap["det_time10"]=&det_time10;
   variableMap["gen_time10"]=&gen_time10;
   variableMap["det_time10_MINUS_gen_time10"]=&det_time10_MINUS_gen_time10;
+  variableMap["pileup_eventID"]=&pileup_eventID;
   variableMap["pileup_muDecayDetID"]=&pileup_muDecayDetID_double;
+  variableMap["pileup_muDecayPosZ"]=&pileup_muDecayPosZ;
+  variableMap["pileup_muDecayPosR"]=&pileup_muDecayPosR;
   variableMap["det_time20"]=&det_time20;
 
   testIVar1=0;
@@ -407,6 +452,9 @@ musrAnalysis::musrAnalysis(TTree *tree)
   motherOfHumanDecayHistograms=NULL;
   humanDecayPileupHistograms=NULL;
   motherOfHumanDecayPileupHistograms=NULL;
+  bool_muDecayTimeTransformation = false;
+  bool_clonedChannelsMultimap_NotEmpty = false;
+  bool_debugingRequired = false;
 
 // if parameter tree is not specified (or zero), connect the file
 // used to generate this class and read the Tree.
diff --git a/musrSimAna/musrCounter.cxx b/musrSimAna/musrCounter.cxx
index 61eea88..a323661 100644
--- a/musrSimAna/musrCounter.cxx
+++ b/musrSimAna/musrCounter.cxx
@@ -2,6 +2,10 @@
 #include "musrCounter.hh"
 #include "TCanvas.h"
 
+typedef std::map<int,int> debugEventMapType;
+debugEventMapType debugEventMap;
+Bool_t bool_debugingRequired;
+
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 musrCounter::musrCounter(int CHANNEL_NR, char CHANNEL_NAME[200], char CHANNEL_TYPE, float E_THRESH, int TIME_SHIFT) {
@@ -17,8 +21,12 @@ musrCounter::musrCounter(int CHANNEL_NR, char CHANNEL_NAME[200], char CHANNEL_TY
   TDC_t1=0;
   TDC_t2=0;
   TDC_histoNrAdd=0;
-  coincidenceTimeWindowMin=0;
-  coincidenceTimeWindowMax=0;
+  antiCoincidenceTimeWindowMin=0;
+  antiCoincidenceTimeWindowMax=0;
+  coincidenceTimeWindowMin_M=0;
+  coincidenceTimeWindowMax_M=0;
+  coincidenceTimeWindowMin_P=0;
+  coincidenceTimeWindowMax_P=0;
   maxCoincidenceTimeWindow=0;
   strcpy(TDC_histoNameAdd,"Unset");
   doubleHitN=0;
@@ -74,7 +82,7 @@ void musrCounter::RemoveHitsInCounter(Long64_t timeBinLimit) {
   //  myPrintThisCounter();
   //  if (counterNr==1) {std::cout<<"ooooo1  timeBinLimit="<<timeBinLimit<<std::endl;  myPrintThisCounter();}
   for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) {
-    //    std::cout<<"musrCounter::RemoveHitsInCounter: counterNr="<<counterNr<<"  hitMap.size()="<<hitMap.size()<<"  maxCoincidenceTimeWindow="<<maxCoincidenceTimeWindow<<" bins of tdcresolution."<<std::endl;
+    //    std::cout<<"  musrCounter::RemoveHitsInCounter:  counterNr="<<counterNr<<"  timeBinLimit="<<timeBinLimit<<"  maxCoincidenceTimeWindow="<<maxCoincidenceTimeWindow<<"    counterTimeShift="<<counterTimeShift<<std::endl;
     if ((it->first)>(timeBinLimit+maxCoincidenceTimeWindow-counterTimeShift)) return;  //note that maxCoincidenceTimeWindow is usually negative number
     else {
       //      std::cout<<"       Deleting hit from counter "<<counterNr<<",  time bin = "<<(it->first)<<std::endl;
@@ -104,19 +112,31 @@ void musrCounter::RewindHitsInCounter(Long64_t timeBinsToRewind) {
 }
 
 //================================================================
-Bool_t musrCounter::IsInCoincidence(Long64_t timeBin, Bool_t ignoreHitsAtBinZero, Long64_t timeBinMinimum, Long64_t timeBinMaximum){
+Bool_t musrCounter::IsInCoincidence(Long64_t timeBin, char motherCounter, Bool_t ignoreHitsAtBinZero, Long64_t timeBinMinimum, Long64_t timeBinMaximum){
   // timeBin ... time bin, at which the coincidence is searched
   // counterTimeShiftOfRequestingCounter  ... time shift (in bin units) of the counter, for which the coincidence is searched
   // ignoreHitsAtBinZero  ... if "true", hits at timeBin will be ignored (needed for searching of coincidence of M counter 
   //                                     with other M counters  or P counters with other P counters)
   //                             "false" should be used for coincidence detectors and vetos.
   if (hitMap.empty()) return false;
-  //  Long64_t timeBinToTest = timeBin + counterTimeShift - counterTimeShiftOfRequestingCounter;
-  Long64_t timeBinToTest = timeBin;
+  //  Long64_t timeBinToTest = timeBin;
+  Long64_t timeBinMin;
+  Long64_t timeBinMax;
+
   // If timeBinMinimum and timeBinMaximum are not specified, use internal time window of the detector (koincidence or veto detectors).
   // Otherwise use timeBinMinimum and timeBinMaximum (e.g.coincidence of a positron counter with other positron counters).
-  Long64_t timeBinMin = (timeBinMinimum==-123456789) ? timeBinToTest + coincidenceTimeWindowMin : timeBinToTest + timeBinMinimum;
-  Long64_t timeBinMax = (timeBinMaximum==-123456789) ? timeBinToTest + coincidenceTimeWindowMax : timeBinToTest + timeBinMaximum;
+  if (timeBinMinimum!=-123456789)  timeBinMin = timeBin + timeBinMinimum;               // time window requested through "timeBinMinimum"
+  else if (counterType == 'V')     timeBinMin = timeBin + antiCoincidenceTimeWindowMin; // this is veto detector
+  else if (motherCounter=='M')     timeBinMin = timeBin + coincidenceTimeWindowMin_M;   // this is coinc. detector (or M counter) connected to M
+  else if (motherCounter=='P')     timeBinMin = timeBin + coincidenceTimeWindowMin_P;   // this is coinc. detector connected to P
+
+  if (timeBinMaximum!=-123456789)  timeBinMax = timeBin + timeBinMaximum;               // time window requested through "timeBinMinimum"
+  else if (counterType == 'V')     timeBinMax = timeBin + antiCoincidenceTimeWindowMax; // this is veto detector
+  else if (motherCounter=='M')     timeBinMax = timeBin + coincidenceTimeWindowMax_M;   // this is coinc. detector (or M counter) connected to M
+  else if (motherCounter=='P')     timeBinMax = timeBin + coincidenceTimeWindowMax_P;   // this is coinc. detector connected to P
+
+  //  Long64_t timeBinMin = (timeBinMinimum==-123456789) ? timeBin + coincidenceTimeWindowMin : timeBin + timeBinMinimum;
+  //  Long64_t timeBinMax = (timeBinMaximum==-123456789) ? timeBin + coincidenceTimeWindowMax : timeBin + timeBinMaximum;
   for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) {
     Long64_t timeBinOfCount_tmp = it->first;
     if ((timeBinOfCount_tmp >= timeBinMin) && (timeBinOfCount_tmp <= timeBinMax)) {
@@ -157,17 +177,14 @@ Bool_t musrCounter::GetNextGoodMuon(Int_t evtID, Long64_t timeBinMin, Long64_t&
     //    std::cout<<"*** "<<evtID<<"   eventNumber="<<eventNumber<<"   canditas="<<numberOfMuonCandidates<<std::endl;
     // Hit candidate was found.  Now check its coincidences and vetos
     for (counterMapType::const_iterator itCounter = koincidenceCounterMap.begin(); itCounter!=koincidenceCounterMap.end(); ++itCounter) {
-      //      if (!( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,counterTimeShift) )) goto endOfThisHit;   // no coincidence found ==> skip hit 
-      if (!( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp) )) goto endOfThisHit;   // no coincidence found ==> skip hit 
+      if (!( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,'M') )) goto endOfThisHit;   // no coincidence found ==> skip hit 
     }
     for (counterMapType::const_iterator itCounter = vetoCounterMap.begin(); itCounter!=vetoCounterMap.end(); ++itCounter) {
-      //      if ( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,counterTimeShift) ) goto endOfThisHit;   //  coincidence with veto found ==> skip hit 
-      if ( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp) ) goto endOfThisHit;   //  coincidence with veto found ==> skip hit 
+      if ( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,'M') ) goto endOfThisHit;   //  coincidence with veto found ==> skip hit 
     }
     // Check coincidences with other hits in the M counter
     // it is expected that there is only one M counter
-    //      if (this->IsInCoincidence(timeBinOfCount_tmp,counterTimeShift,true) ) {   //  coincidence with another M-counter hit ==> skip hit 
-    if (this->IsInCoincidence(timeBinOfCount_tmp,true) ) {   //  coincidence with another M-counter hit ==> skip hit 
+    if (this->IsInCoincidence(timeBinOfCount_tmp,'M',true) ) {   //  coincidence with another M-counter hit ==> skip hit 
       doubleHitFound = true;
       //     std::cout<<"UJGUR double hit found"<<std::endl;
       goto endOfThisHit;   
@@ -187,35 +204,58 @@ Bool_t musrCounter::GetNextGoodMuon(Int_t evtID, Long64_t timeBinMin, Long64_t&
 }
 
 //================================================================
-Bool_t musrCounter::GetNextGoodPositron(Int_t evtID, Long64_t timeBinMin, Long64_t timeBinMax, Long64_t& timeBinOfNextGoodHit, Long64_t& timeBinOfNextGoodHit_phaseShifted, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep, Bool_t& doubleHitFound) {
+Int_t musrCounter::GetNextGoodPositron(Int_t evtID, Long64_t timeBinMin, Long64_t timeBinMax, Long64_t& timeBinOfNextGoodHit, Long64_t& timeBinOfNextGoodHit_phaseShifted, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep, Bool_t& doubleHitFound) {
   // INPUT PARAMETERS:  evtID, timeBinMin
   // OUTPUT PARAMETERS: timeBinOfNextGoodHit
-  //
+  //                    positronQuality = 0 ... no positron candidate found
+  //                                    = 1 ... positron candidate found, but vetoed or not in a required coincidence with other detector
+  //                                    = 2 ... good positron found
   // Loop over the hits in the counter
-  if (hitMap.empty()) return false;
+  Int_t positronQuality=0;
+  if (bool_debugingRequired) {if (debugEventMap[evtID]>4) myPrintThisCounter(evtID);}
+  if (hitMap.empty()) return 0;
   if (counterType!='P') {std::cout<<"\n!!! FATAL ERROR !!! musrCounter::GetNextGoodPositron: not the positron counter! ==> S T O P !!!\n"; exit(1);}
   doubleHitFound = false;
   for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) {
     Int_t    eventNumber = (it->second)->eventIDnumber;
     Long64_t timeBinOfCount_tmp = it->first;
-    if ((timeBinOfCount_tmp <= timeBinMin) || (timeBinOfCount_tmp > timeBinMax)) continue;   // This hit is out of the data interval ==> skip it
+    if ((timeBinOfCount_tmp <= timeBinMin) || (timeBinOfCount_tmp > timeBinMax)) {
+      if (bool_debugingRequired) {
+	if (debugEventMap[evtID]>3) {std::cout<<"DEBUGEVENT:"<<evtID<<"\"GetNextGoodPositron\":  Hit out of data interval"<<std::endl;}
+      }
+      continue;   // This hit is out of the data interval ==> skip it
+    }
 
     // Hit candidate was found.  Now check its coincidences and vetos
+    positronQuality=1;
     for (counterMapType::const_iterator itCounter = koincidenceCounterMap.begin(); itCounter!=koincidenceCounterMap.end(); ++itCounter) {
-      //      if (!( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,counterTimeShift) )) goto endOfThisHit;   // no coincidence found ==> skip hit 
-      if (!( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp) )) goto endOfThisHit;   // no coincidence found ==> skip hit 
+      if (bool_debugingRequired) {
+	if (debugEventMap[evtID]>4) { (itCounter->second)->myPrintThisCounter(evtID); }
+      }
+      if (!( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,'P') )) {
+	if (bool_debugingRequired) {
+	  if (debugEventMap[evtID]>3) {std::cout<<"DEBUGEVENT:"<<evtID<<"\"GetNextGoodPositron\":  Coincidence required but not found"<<std::endl;}
+	}
+	goto endOfThisHit;   // no coincidence found ==> skip hit 
+      }
     }
     for (counterMapType::const_iterator itCounter = vetoCounterMap.begin(); itCounter!=vetoCounterMap.end(); ++itCounter) {
-      //      if ( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,counterTimeShift) ) goto endOfThisHit;   //  coincidence with veto found ==> skip hit 
-      if ( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp) ) goto endOfThisHit;   //  coincidence with veto found ==> skip hit 
+      if ( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,'P') ) {
+	if (bool_debugingRequired) {
+	  if (debugEventMap[evtID]>3) {std::cout<<"DEBUGEVENT:"<<evtID<<"\"GetNextGoodPositron\":  Coincidence vith veto detector found"<<std::endl;}
+	}
+	goto endOfThisHit;   //  coincidence with veto found ==> skip hit
+      } 
     }
     // Check coincidences with other P counters    
     // Coincidences with other P counters must be checked in musrAnalysis.cxx
-    //      if (this->IsInCoincidence(timeBinOfCount_tmp,counterTimeShift,true) ) {  //  coincidence with another P-counter hit ==> skip hit 
-    // hovno
+    //
     // ADD HERE THE CHECK THAT THE POSITRON IS FOUND WITHIN THE DATA INTERVAL, i.e. within <timeBinMin,timeBinMax>.
     // AND ALSO CHECK THAT THIS IS TRUE FOR THE COINCIDENCE WITH OTHER POSITRON COUNTERS!!!
-    if (this->IsInCoincidence(timeBinOfCount_tmp,true,timeBinMin,timeBinMax) ) {  //  coincidence with another P-counter hit ==> skip hit 
+    if (this->IsInCoincidence(timeBinOfCount_tmp,'P',true,timeBinMin,timeBinMax) ) {  //  coincidence with another P-counter hit ==> skip hit 
+      	if (bool_debugingRequired) {
+	  if (debugEventMap[evtID]>3) {std::cout<<"DEBUGEVENT:"<<evtID<<"\"GetNextGoodPositron\":  Coincidence with other positron candidate in the same counter."<<std::endl;}
+	}
       doubleHitFound = true;
       //	std::cout<<"tttttttttttttttttttt doubleHitN="<<doubleHitN++<<std::endl;
       goto endOfThisHit;   
@@ -227,50 +267,57 @@ Bool_t musrCounter::GetNextGoodPositron(Int_t evtID, Long64_t timeBinMin, Long64
     idetEdep = (it->second)->det_edep;
     timeBinOfNextGoodHit = timeBinOfCount_tmp;
     timeBinOfNextGoodHit_phaseShifted = (it->second) -> timeBin2;
+    if (bool_debugingRequired) {
+      if (debugEventMap[evtID]>3) {std::cout<<"DEBUGEVENT:"<<evtID<<"\"GetNextGoodPositron\":  Good positron candidate found in this counter."<<std::endl;}
+    }
     //cDEL    std::cout<<"timeBinOfNextGoodHit ="<<timeBinOfNextGoodHit<<"        timeBinOfNextGoodHit_phaseShifted = "<<timeBinOfNextGoodHit_phaseShifted<<std::endl;
-    return true;
+    return 2;
     
     endOfThisHit:
-    continue;
+    ;
+    //    continue;
   }
-  return false;
+  return positronQuality;
 }
 
 //================================================================
-void musrCounter::SetCoincidenceTimeWindowOfAllCoincidenceDetectors(Long64_t maxCoinc, Long64_t min, Long64_t max) {
+void musrCounter::SetCoincidenceTimeWindowOfAllCoincidenceDetectors(char motherCounter, Long64_t maxCoinc, Long64_t min, Long64_t max) {
+  //  std::cout<<"QQQQQQQQQQQQQQQQQQQQQ koincidenceCounterMap.size()="<<koincidenceCounterMap.size()<<std::endl;
   for (counterMapType::const_iterator it = koincidenceCounterMap.begin(); it!=koincidenceCounterMap.end(); ++it) {
-    if ( ( ((it->second)->GetMaxCoincidenceTimeWindow()) !=0) &&
-	 ( ((it->second)->GetMaxCoincidenceTimeWindow()) !=maxCoinc) )  {
-      std::cout<<" !!!!  ERROR  SetCoincidenceTimeWindowOfAllCoincidenceDetectors : coincidenceTimeWindow set multiple times! ==> S T O P !!!"<<std::endl;
+    Long64_t maxCoinc_AlreadySet = ((it->second)->GetMaxCoincidenceTimeWindow());
+    if (maxCoinc < maxCoinc_AlreadySet)    (it->second)->SetMaxCoincidenceTimeWindow(maxCoinc);
+ 
+    if      (motherCounter=='M') (it->second)->SetCoincidenceTimeWindow_M(min,max);
+    else if (motherCounter=='P') (it->second)->SetCoincidenceTimeWindow_P(min,max);
+    else {
+      std::cout<<"musrCounter::SetCoincidenceTimeWindowOfAllCoincidenceDetectors ERROR: Strange motherCounter "
+	       <<motherCounter<<"\n   ==> S T O P "<<std::endl; 
       exit(1);
     }
-    (it->second)->SetMaxCoincidenceTimeWindow(maxCoinc);
-    (it->second)->SetCoincidenceTimeWindow(min,max);
   }
 }
 
 //================================================================
 void musrCounter::SetCoincidenceTimeWindowOfAllVetoDetectors(Long64_t maxCoinc, Long64_t min, Long64_t max) {
   for (counterMapType::const_iterator it = vetoCounterMap.begin(); it!=vetoCounterMap.end(); ++it) {
-    if ( ( ((it->second)->GetMaxCoincidenceTimeWindow()) !=0) &&
-	 ( ((it->second)->GetMaxCoincidenceTimeWindow()) !=maxCoinc) )  {
-      std::cout<<" !!!!  ERROR  SetCoincidenceTimeWindowOfAllVetoDetectors : coincidenceTimeWindow set multiple times! ==> S T O P !!!"<<std::endl;
-      exit(1);
-    }
-    (it->second)->SetMaxCoincidenceTimeWindow(maxCoinc);
-    (it->second)->SetCoincidenceTimeWindow(min,max);
+    Long64_t maxCoinc_AlreadySet = ((it->second)->GetMaxCoincidenceTimeWindow());
+    if (maxCoinc < maxCoinc_AlreadySet)    (it->second)->SetMaxCoincidenceTimeWindow(maxCoinc);
+
+    (it->second)->SetAntiCoincidenceTimeWindow(min,max);
   }
 }
 
 //================================================================
-void musrCounter::myPrintThisCounter(Int_t evtID) {
+void musrCounter::myPrintThisCounter(Int_t evtID, Int_t detail) {
   Bool_t eventMixing=false;
-  std::cout<<"musrCounter::myPrintThisCounter:  counterNr="<<counterNr<<":   ";
+  if ((hitMap.begin()==hitMap.end()) && (detail<=1) ) return;
+  if (detail>1) std::cout<<"musrCounter::myPrintThisCounter:  counterNr = "<<counterNr<<":   ";
+  else          std::cout<<"   counter = "<<counterNr<<":   ";
   for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) {
     std::cout<<"\t"<<it->first<<"  "<<(it->second)->eventIDnumber<<",";
     if (evtID != (it->second)->eventIDnumber) {eventMixing=true;}
   }
-  if (eventMixing) {std::cout<<"   Potential event mmmmmmmmixing";}
+  if (eventMixing) {std::cout<<"   Potential event mixing";}
   std::cout<<std::endl;
 }
 
diff --git a/musrSimAna/musrCounter.hh b/musrSimAna/musrCounter.hh
index 54cec8c..89e4a84 100644
--- a/musrSimAna/musrCounter.hh
+++ b/musrSimAna/musrCounter.hh
@@ -18,6 +18,11 @@ public:
   Int_t det_id;
   Double_t det_edep;
   Long64_t timeBin2;
+
+  //  extern double GlobalKamil;
+  //  typedef std::map<int,int> debugEventMapType;
+  //  extern debugEventMapType debugEventMap;
+  //  extern Bool_t bool_debugingRequired;
 };
 
 class musrCounter {
@@ -39,19 +44,21 @@ class musrCounter {
     }
     void SetMaxCoincidenceTimeWindow(Long64_t val) {maxCoincidenceTimeWindow=val;}
     Long64_t GetMaxCoincidenceTimeWindow() {return maxCoincidenceTimeWindow;}
-    void SetCoincidenceTimeWindowOfAllCoincidenceDetectors(Long64_t maxCoinc, Long64_t min, Long64_t max);
+    void SetCoincidenceTimeWindowOfAllCoincidenceDetectors(char motherCounter, Long64_t maxCoinc, Long64_t min, Long64_t max);
     void SetCoincidenceTimeWindowOfAllVetoDetectors(Long64_t maxCoinc, Long64_t min, Long64_t max);
-    void SetCoincidenceTimeWindow(Long64_t min, Long64_t max) {coincidenceTimeWindowMin=min; coincidenceTimeWindowMax=max;}
+    void SetCoincidenceTimeWindow_M(Long64_t min, Long64_t max) {coincidenceTimeWindowMin_M=min; coincidenceTimeWindowMax_M=max;}
+    void SetCoincidenceTimeWindow_P(Long64_t min, Long64_t max) {coincidenceTimeWindowMin_P=min; coincidenceTimeWindowMax_P=max;}
+    void SetAntiCoincidenceTimeWindow(Long64_t min, Long64_t max) {antiCoincidenceTimeWindowMin=min; antiCoincidenceTimeWindowMax=max;}
     void SetTDChistogram(char hName[200],int t0,int t1,int t2,int hNr,char hNameAdd[200]);
     void FillTDChistogram(Double_t variable, Double_t vaha);
     void DrawTDChistogram();
     void FillHitInCounter(Double_t edep, Long64_t timeBin, Long64_t timeBin2, Int_t kEntry, Int_t eveID, Int_t iDet, Int_t detectorID);
     void RemoveHitsInCounter(Long64_t timeBinLimit);
     void RewindHitsInCounter(Long64_t timeBinsToRewind);
-    Bool_t IsInCoincidence(Long64_t timeBin, Bool_t ignoreHitsAtBinZero=false, Long64_t timeBinMinimum=-123456789, Long64_t timeBinMaximum=-123456789);
+    Bool_t IsInCoincidence(Long64_t timeBin, char motherCounter, Bool_t ignoreHitsAtBinZero=false, Long64_t timeBinMinimum=-123456789, Long64_t timeBinMaximum=-123456789);
     Bool_t GetNextGoodMuon(Int_t evtID, Long64_t timeBinMin, Long64_t& timeBinOfNextHit, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep, Bool_t& doubleHitFound);
-    Bool_t GetNextGoodPositron(Int_t evtID, Long64_t timeBinMin, Long64_t timeBinMax, Long64_t& timeBinOfNextGoodHit, Long64_t& timeBinOfNextGoodHit_phaseShifted, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep, Bool_t& doubleHitFound);
-    void myPrintThisCounter(Int_t evtID);
+    Int_t GetNextGoodPositron(Int_t evtID, Long64_t timeBinMin, Long64_t timeBinMax, Long64_t& timeBinOfNextGoodHit, Long64_t& timeBinOfNextGoodHit_phaseShifted, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep, Bool_t& doubleHitFound);
+    void myPrintThisCounter(Int_t evtID,  Int_t detail=2);
     Long64_t GetNumberOfMuonCandidates(){return numberOfMuonCandidates;}
 
 
@@ -69,8 +76,8 @@ class musrCounter {
   int   TDC_t0, TDC_t1, TDC_t2, TDC_histoNrAdd;
   char  TDC_histoNameAdd[200];
   TH1D* histTDC;
-  Long64_t coincidenceTimeWindowMin;
-  Long64_t coincidenceTimeWindowMax;
+  Long64_t antiCoincidenceTimeWindowMin, coincidenceTimeWindowMin_M, coincidenceTimeWindowMin_P;
+  Long64_t antiCoincidenceTimeWindowMax, coincidenceTimeWindowMax_M, coincidenceTimeWindowMax_P;
   Long64_t maxCoincidenceTimeWindow;
   //  typedef std::map<Long64_t,Int_t> hitMap_TYPE;  // Long64_t = timeBin, Int_t=eventID
   typedef std::map<Long64_t,hitInfo*> hitMap_TYPE;  // Long64_t = timeBin, hitInfo = eventID and det_i
diff --git a/musrSimAna/musrTH.cxx b/musrSimAna/musrTH.cxx
index 095dbfb..01f95ac 100644
--- a/musrSimAna/musrTH.cxx
+++ b/musrSimAna/musrTH.cxx
@@ -47,7 +47,8 @@ void musrTH::FillTH1D(Double_t vaha, Bool_t* cond){
   for (Int_t i=0; i<musrAnalysis::nrConditions; i++) {
     if (bool_rotating_reference_frame) {
       //      Double_t var = *variableToBeFilled_X;
-      if (cond[i])  histArray1D[i]->Fill(*variableToBeFilled_X,vaha*cos(2*musrAnalysis::pi*rot_ref_frequency*(*variableToBeFilled_X)+rot_ref_phase));
+      Double_t waha = vaha*exp((*variableToBeFilled_X)/2.19703)*cos(2*musrAnalysis::pi*rot_ref_frequency*(*variableToBeFilled_X)+rot_ref_phase);
+      if (cond[i])  histArray1D[i]->Fill(*variableToBeFilled_X,waha);
     }
     else {
       if (cond[i])  histArray1D[i]->Fill(*variableToBeFilled_X,vaha);
@@ -223,7 +224,7 @@ void  musrTH::FitHistogramsIfRequired(Double_t omega) {
       
   for (Int_t i=0; i<musrAnalysis::nrConditions; i++) {
     //    std::cout<<"fitted histogram pointer="<<histArray1D[i]<<std::endl;
-    histArray1D[i]->Fit(funct,"","",funct_xMin,funct_xMax);
+    histArray1D[i]->Fit(funct,"WW","",funct_xMin,funct_xMax);
     //    histArray1D[i]->Fit(funct,"LL","",funct_xMin,funct_xMax);
   }
   //  if (strcmp(funct->GetName(),"simpleExpoPLUSconst")==0) {
diff --git a/src/musrDetectorConstruction.cc b/src/musrDetectorConstruction.cc
index d866237..49d671d 100644
--- a/src/musrDetectorConstruction.cc
+++ b/src/musrDetectorConstruction.cc
@@ -64,6 +64,9 @@
 #include "G4RegionStore.hh"
 #include "G4ProductionCuts.hh"
 
+//#include "G4OpticalSurface.hh"
+#include "G4LogicalBorderSurface.hh"
+
 #include "musrRootOutput.hh"   //cks for storing some info in the Root output file
 #include "musrParameters.hh"
 #include "musrErrorMessage.hh"
@@ -71,9 +74,10 @@
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-musrDetectorConstruction::musrDetectorConstruction()
-:parameterFileName("Unset"), checkOverlap(true), aScintSD(0)
+musrDetectorConstruction::musrDetectorConstruction(G4String steeringFileName)
+:checkOverlap(true), aScintSD(0)
 {  
+  parameterFileName = steeringFileName;
   DefineMaterials();
   detectorMessenger = new musrDetectorMessenger(this);
 }
@@ -102,7 +106,9 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
   //  G4double roundingErr=0.01*mm;        // 0.01mm precision is OK for displaying subtracted volumes, while 0.001mm is not
   //----------------------
 
-  musrRootOutput* myRootOutput = musrRootOutput::GetRootInstance();
+  musrRootOutput*     myRootOutput     = musrRootOutput::GetRootInstance();
+  musrSteppingAction* mySteppingAction = musrSteppingAction::GetInstance();
+
 
   G4VPhysicalVolume* pointerToWorldVolume=NULL;
   //  Read detector configuration parameters from the steering file:
@@ -483,7 +489,6 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	if (strstr(name,"save")!=NULL) {
 	  if (volumeID!=0) {    // do not store hits in special "save" volume if ID of this volume is 0 
 	                        // (due to difficulties to distinguish between ID=0 and no save volume when using std::map)
-	    musrSteppingAction* mySteppingAction = musrSteppingAction::GetInstance();
 	    mySteppingAction->SetLogicalVolumeAsSpecialSaveVolume(logicName,volumeID);
 	    musrRootOutput::GetRootInstance()->SetSpecialSaveVolumeDefined();
 	  }
@@ -562,6 +567,164 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	}
       }
 
+      // cks:  Optical Boundary  (needed only when simulating optical photons)
+      else if (strcmp(tmpString1,"opticalSurface")==0) {
+	if (musrParameters::boolG4OpticalPhotons) {
+	  char optSurfaceName[100];
+	  char physVolName1[100];
+	  char physVolName2[100];
+	  char type[100];
+	  char finish[100];
+	  char model[100];
+	  char materialPropertiesTableName[100]="Undefined";
+	  sscanf(&line[0],"%*s %*s %s %s %s %s %s %s %s",optSurfaceName,physVolName1,physVolName2,type,finish,model,materialPropertiesTableName);
+	  G4VPhysicalVolume* pPhysVol1 = FindPhysicalVolume(physVolName1);
+	  G4VPhysicalVolume* pPhysVol2 = FindPhysicalVolume(physVolName2);
+	  if ((pPhysVol1==NULL)||(pPhysVol2==NULL)) {
+	    G4cout << "ERROR! musrDetectorConstruction::Construct():  Physical Volume not found!"<<G4endl;
+	    G4cout << "      ==>  S T O P    F O R C E D"<<G4endl;
+	    ReportGeometryProblem(line);
+	    exit(1);
+	  }
+	  
+	  G4OpticalSurface* optSurfTMP = new G4OpticalSurface(optSurfaceName);
+	  new G4LogicalBorderSurface(optSurfaceName,pPhysVol1,pPhysVol2,optSurfTMP);
+	  
+	  std::map<std::string,G4SurfaceType>          OpticalTypeMap;
+	  std::map<std::string,G4OpticalSurfaceModel>  OpticalModelMap;
+	  std::map<std::string,G4OpticalSurfaceFinish> OpticalFinishMap;
+	  
+	  OpticalTypeMap["dielectric_metal"]=dielectric_metal;             // dielectric-metal interface
+	  OpticalTypeMap["dielectric_dielectric"]=dielectric_dielectric;   // dielectric-dielectric interface
+	  OpticalTypeMap["dielectric_LUT"]=dielectric_LUT;                 // dielectric-Look-Up-Table interface
+	  OpticalTypeMap["firsov"]=firsov;                                 // for Firsov Process
+	  OpticalTypeMap["x_ray"]=x_ray;                                   // for x-ray mirror process
+	  
+	  OpticalModelMap["glisur"]=glisur;                                // original GEANT3 model
+	  OpticalModelMap["unified"]=unified;                              // UNIFIED model
+	  OpticalModelMap["LUT"]=LUT;                                      // Look-Up-Table model
+	  
+	  OpticalFinishMap["polished"]=polished;                           // smooth perfectly polished surface
+	  OpticalFinishMap["polishedfrontpainted"]=polishedfrontpainted;   // smooth top-layer (front) paint
+	  OpticalFinishMap["polishedbackpainted"]=polishedbackpainted;     // same is 'polished' but with a back-paint
+	  OpticalFinishMap["ground"]=ground;                               // rough surface
+	  OpticalFinishMap["groundfrontpainted"]=groundfrontpainted;       // rough top-layer (front) paint
+	  OpticalFinishMap["groundbackpainted"]=groundbackpainted;         // same as 'ground' but with a back-paint  
+	  OpticalFinishMap["polishedlumirrorair"]=polishedlumirrorair;     // mechanically polished surface, with lumirror
+	  OpticalFinishMap["polishedlumirrorglue"]=polishedlumirrorglue;   // mechanically polished surface, with lumirror & meltmount
+	  OpticalFinishMap["polishedair"]=polishedair;                     // mechanically polished surface
+	  OpticalFinishMap["polishedteflonair"]=polishedteflonair;         // mechanically polished surface, with teflon
+	  OpticalFinishMap["polishedtioair"]=polishedtioair;               // mechanically polished surface, with tio paint
+	  OpticalFinishMap["polishedtyvekair"]=polishedtyvekair;           // mechanically polished surface, with tyvek
+	  OpticalFinishMap["polishedvm2000air"]=polishedvm2000air;         // mechanically polished surface, with esr film
+	  OpticalFinishMap["polishedvm2000glue"]=polishedvm2000glue;       // mechanically polished surface, with esr film & meltmount
+	  OpticalFinishMap["etchedlumirrorair"]=etchedlumirrorair;         // chemically etched surface, with lumirror
+	  OpticalFinishMap["etchedlumirrorglue"]=etchedlumirrorglue;       // chemically etched surface, with lumirror & meltmount
+	  OpticalFinishMap["etchedair"]=etchedair;                         // chemically etched surface
+	  OpticalFinishMap["etchedteflonair"]=etchedteflonair;             // chemically etched surface, with teflon
+	  OpticalFinishMap["etchedtioair"]=etchedtioair;                   // chemically etched surface, with tio paint
+	  OpticalFinishMap["etchedtyvekair"]=etchedtyvekair;               // chemically etched surface, with tyvek
+	  OpticalFinishMap["etchedvm2000air"]=etchedvm2000air;             // chemically etched surface, with esr film
+	  OpticalFinishMap["etchedvm2000glue"]=etchedvm2000glue;           // chemically etched surface, with esr film & meltmount  
+	  OpticalFinishMap["groundlumirrorair"]=groundlumirrorair;         // rough-cut surface, with lumirror
+	  OpticalFinishMap["groundlumirrorglue"]=groundlumirrorglue;       // rough-cut surface, with lumirror & meltmount
+	  OpticalFinishMap["groundair"]=groundair;                         // rough-cut surface
+	  OpticalFinishMap["groundteflonair"]=groundteflonair;             // rough-cut surface, with teflon
+	  OpticalFinishMap["groundtioair"]=groundtioair;                   // rough-cut surface, with tio paint
+	  OpticalFinishMap["groundtyvekair"]=groundtyvekair;               // rough-cut surface, with tyvek
+	  OpticalFinishMap["groundvm2000air"]=groundvm2000air;             // rough-cut surface, with esr film
+	  OpticalFinishMap["groundvm2000glue"]=groundvm2000glue;           // rough-cut surface, with esr film & meltmount
+	  
+	  
+	  G4SurfaceType OpticalType            = OpticalTypeMap[type];
+	  G4OpticalSurfaceModel OpticalModel   = OpticalModelMap[model];
+	  G4OpticalSurfaceFinish OpticalFinish = OpticalFinishMap[finish];
+	  if ((OpticalType==0)&&(strcmp(type,"dielectric_metal")!=0)) {
+	    G4cout<<"ERROR!  musrDetectorConstruction::Construct(): Optical type \""<<type<<"\" not found!"<<G4endl;
+	    G4cout << "      ==>  S T O P    F O R C E D"<<G4endl;
+	    G4cout<<"   "<<line;
+	    exit(1);
+	  }
+	  if ((OpticalModel==0)&&(strcmp(model,"glisur")!=0)) {
+	    G4cout<<"ERROR!  musrDetectorConstruction::Construct(): Optical surface model \""<<model<<"\" not found!"<<G4endl;
+	    G4cout << "      ==>  S T O P    F O R C E D"<<G4endl;
+	    G4cout<<"   "<<line;
+	    exit(1);
+	  }
+	  if ((OpticalFinish==0)&&(strcmp(finish,"polished")!=0)) {
+	    G4cout<<"ERROR!  musrDetectorConstruction::Construct(): Optical surface finish \""<<finish<<"\" not found!"<<G4endl;
+	    G4cout << "      ==>  S T O P    F O R C E D"<<G4endl;
+	    G4cout<<"   "<<line;
+	    exit(1);
+	  }
+
+
+	  optSurfTMP->SetType(OpticalType);
+	  optSurfTMP->SetFinish(OpticalFinish);
+	  optSurfTMP->SetModel(OpticalModel);
+	  // Assign the "Material properties table" if required by the user:
+	  G4cout<<"materialPropertiesTableName="<<materialPropertiesTableName<<G4endl;
+	  if (strcmp(materialPropertiesTableName,"Undefined")!=0) {
+	    G4MaterialPropertiesTable* MPT_tmp=NULL;
+	    itMPT = materialPropertiesTableMap.find(materialPropertiesTableName);
+	    if (itMPT==materialPropertiesTableMap.end()) {  // G4MaterialPropertiesTable of this name does not exist --> problem  
+	      G4cout<<"\n\n musrDetectorConstruction(): Material Properties Table \""<<materialPropertiesTableName
+		    <<"\" should be assigned to G4OpticalSurface \""<<optSurfaceName<<"\""<<G4endl;
+	      G4cout<<" but the table was not defined yet (by command /musr/command materialPropertiesTable )"<<G4endl;
+	      G4cout << "S T O P    F O R C E D"<<G4endl;
+	      G4cout << line << G4endl;
+	      exit(1);
+	    }
+	    else {
+	      MPT_tmp = itMPT->second;
+	    }
+	    optSurfTMP->SetMaterialPropertiesTable(MPT_tmp);
+	    G4cout<<optSurfTMP<<G4endl;
+	    optSurfTMP->GetMaterialPropertiesTable()->DumpTable();
+	  }
+	  G4cout<<"Optical surface \""<<optSurfaceName<<"\" created.  OpticalType="<<OpticalType
+		<<"  OpticalFinish="<<OpticalFinish<<"  OpticalModel="<<OpticalModel<<G4endl;
+	}
+      }
+
+      else if (strcmp(tmpString1,"OPSA")==0){   // optical photon signal analysis
+	if (musrParameters::boolG4OpticalPhotons) {
+	  musrScintSD* myMusrScintSD = musrScintSD::GetInstance();
+	  if (myMusrScintSD==NULL) {
+	    sprintf(eMessage,"musrDetectorConstruction.cc::Construct():  musrScintSD::GetInstance() is NULL - no musr/ScintSD set?");
+	    musrErrorMessage::GetInstance()->musrError(FATAL,eMessage,false);
+	  }
+	  char varName[100];
+	  float fVarValue;
+	  int   iVarValue;
+	  sscanf(&line[0],"%*s %*s %s",varName);
+	  if (strcmp(varName,"minNrOfDetectedPhotons")==0) {
+	    sscanf(&line[0],"%*s %*s %*s %d",&iVarValue);
+	    myMusrScintSD -> Set_OPSA_minNrOfDetectedPhotons(iVarValue);
+	  }
+	  else if (strcmp(varName,"signalSeparationTime")==0) {
+	    sscanf(&line[0],"%*s %*s %*s %g",&fVarValue);
+	    myMusrScintSD -> Set_OPSA_SignalSeparationTime(fVarValue*nanosecond);
+	  }
+	  else if (strcmp(varName,"photonFractions")==0) {
+	    float a, b, c, d, e;
+	    sscanf(&line[0],"%*s %*s %*s %g %g %g %g %g",&a, &b, &c, &d, &e);
+	    myMusrScintSD -> Set_OPSA_frac(a,b,c,d,e);
+	  }
+	  else if (strcmp(varName,"eventsForOPSAhistos")==0) {
+	    int i_eventID, i_detectorID;
+	    sscanf(&line[0],"%*s %*s %*s %d %d",&i_eventID,&i_detectorID);
+	    myMusrScintSD -> AddEventIDToMultimapOfEventIDsForOPSAhistos(i_eventID,i_detectorID);
+	  }
+	  else if (strcmp(varName,"OPSAhist")==0) {
+	    int nBins;
+	    float min, max;
+	    sscanf(&line[0],"%*s %*s %*s %d %g %g",&nBins,&min,&max);
+	    myMusrScintSD -> SetOPSAhistoBinning(nBins,min*nanosecond,max*nanosecond);
+	  }
+	}
+      }
+
       // cks:  Implementation of the Global Field (to allow overlapping fields) based
       //       on the Peter Gumplinger's implementation of G4BeamLine code into Geant4.
       //
@@ -823,7 +986,6 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	  int eventWeight;
 	  char tmpLogVolName[100];
 	  sscanf(&line[0],"%*s %*s %*s %s %d",tmpLogVolName,&eventWeight);
-	  musrSteppingAction* mySteppingAction = musrSteppingAction::GetInstance();
 	  mySteppingAction -> SetVolumeForMuonEventReweighting(G4String(tmpLogVolName),eventWeight);
 	}
 	else {
@@ -979,6 +1141,15 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	  if (strcmp(tmpString2,"det_VvvProcID")==0){musrRootOutput::store_det_VvvProcID = false;}
 	  if (strcmp(tmpString2,"det_VvvTrackID")==0){musrRootOutput::store_det_VvvTrackID = false;}
 	  if (strcmp(tmpString2,"det_VvvParticleID")==0){musrRootOutput::store_det_VvvParticleID = false;}
+	  if (strcmp(tmpString2,"odet_ID")==0)       {musrRootOutput::store_odet_ID = false;}
+	  if (strcmp(tmpString2,"odet_nPhot")==0)    {musrRootOutput::store_odet_nPhot = false;}
+	  if (strcmp(tmpString2,"odet_timeFirst")==0) {musrRootOutput::store_odet_timeFirst = false;}
+	  if (strcmp(tmpString2,"odet_timeA")==0)    {musrRootOutput::store_odet_timeA = false;}
+	  if (strcmp(tmpString2,"odet_timeB")==0)    {musrRootOutput::store_odet_timeB = false;}
+	  if (strcmp(tmpString2,"odet_timeC")==0)    {musrRootOutput::store_odet_timeC = false;}
+	  if (strcmp(tmpString2,"odet_timeD")==0)    {musrRootOutput::store_odet_timeD = false;}
+	  if (strcmp(tmpString2,"odet_timeE")==0)    {musrRootOutput::store_odet_timeE = false;}
+	  if (strcmp(tmpString2,"odet_timeLast")==0) {musrRootOutput::store_odet_timeLast = false;}
 	}
 	else if(strcmp(storeIt,"on")==0) {
 	  if (strcmp(tmpString2,"fieldIntegralBx")==0){musrRootOutput::store_fieldIntegralBx = true;}
@@ -991,12 +1162,13 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	  if ((musrRootOutput::store_fieldIntegralBx)||(musrRootOutput::store_fieldIntegralBy)||
 	      (musrRootOutput::store_fieldIntegralBz)||(musrRootOutput::store_fieldIntegralBz1)||
 	      (musrRootOutput::store_fieldIntegralBz2)||(musrRootOutput::store_fieldIntegralBz3) ) {
-	    musrSteppingAction::GetInstance()->SetCalculationOfFieldIntegralRequested(true);
+	    mySteppingAction -> SetCalculationOfFieldIntegralRequested(true);
 	  }
 	}
       }
       
-      else if (strcmp(tmpString1,"process")==0) {
+      else if ((strcmp(tmpString1,"process")==0)||(strcmp(tmpString1,"G4OpticalPhotons")==0)
+	       ||(strcmp(tmpString1,"materialPropertiesTable")==0)||(strcmp(tmpString1,"setMaterialPropertiesTable")==0)) {
 	;  // processes are interpreded later in  musrPhysicsList.cc
       }
 
@@ -1133,49 +1305,140 @@ void musrDetectorConstruction::DefineMaterials()
   
   new G4Material("ArgonGas", z= 18., a= 39.95*g/mole, density= 0.00000000001*mg/cm3);
 
-  if (musrParameters::boolG4OpticalPhotons) {
-    G4NistManager* man = G4NistManager::Instance();
-    G4Material* scintik = man->FindOrBuildMaterial("G4_PLASTIC_SC_VINYLTOLUENE");   
-    //    G4Material* scintik = G4Material::GetMaterial("G4_PLASTIC_SC_VINYLTOLUENE");
-    G4cout<<"scintik="<<scintik<<G4endl;
-    if (scintik!=NULL) {
-      const G4int nEntries = 14;
-      G4double PhotonEnergy[nEntries] = 
-	{ 2.695*eV,   2.75489*eV, 2.8175*eV,  2.88302*eV,                 // 460, 450, 440, 430 nm
-	  2.95167*eV, 3.02366*eV, 3.09925*eV, 3.17872*eV, 3.26237*eV,     // 420, 410, 400, 390, 380 nm
-	  3.30587*eV, 3.35054*eV, 3.44361*eV, 3.542*eV,   3.64618*eV };   // 375, 370, 360, 350, 340 nm
-      G4double RefractiveIndex[nEntries] = 
-	{ 1.58, 1.58, 1.58, 1.58, 
-	  1.58, 1.58, 1.58, 1.58, 1.58,
-	  1.58, 1.58, 1.58, 1.58, 1.58  };
-      G4double Absorption[nEntries] =
-	{ 8*cm, 8*cm, 8*cm, 8*cm, 
-	  8*cm, 8*cm, 8*cm, 8*cm, 8*cm,
-	  8*cm, 8*cm, 8*cm, 8*cm, 8*cm  };
-      G4double ScintilFast[nEntries] =
-	{ 0.01,  0.07, 0.15, 0.26,
-	  0.375, 0.52, 0.65, 0.80, 0.95,
-          1,     0.88, 0.44, 0.08, 0.01 };
-      G4double ScintilSlow[nEntries] =
-	{ 0.01,  0.07, 0.15, 0.26,
-	  0.375, 0.52, 0.65, 0.80, 0.95,
-          1,     0.88, 0.44, 0.08, 0.01 };
-      G4MaterialPropertiesTable* myMPT1 = new G4MaterialPropertiesTable();
-      myMPT1->AddProperty("RINDEX",        PhotonEnergy, RefractiveIndex, nEntries);
-      myMPT1->AddProperty("ABSLENGTH",     PhotonEnergy, Absorption,      nEntries);
-      myMPT1->AddProperty("FASTCOMPONENT", PhotonEnergy, ScintilFast,     nEntries);
-      myMPT1->AddProperty("SLOWCOMPONENT", PhotonEnergy, ScintilSlow,     nEntries);
-      myMPT1->AddConstProperty("SCINTILLATIONYIELD", 8400./MeV);
-      myMPT1->AddConstProperty("RESOLUTIONSCALE",1.0);
-      myMPT1->AddConstProperty("FASTTIMECONSTANT",1.6*ns);
-      myMPT1->AddConstProperty("SLOWTIMECONSTANT",1.6*ns);
-      myMPT1->AddConstProperty("YIELDRATIO",1.0);
-      scintik->SetMaterialPropertiesTable(myMPT1);
 
-      scintik->GetMaterialPropertiesTable()->DumpTable();
+  if (musrParameters::boolG4OpticalPhotons) {
+
+    FILE *fSteeringFile=fopen(parameterFileName.c_str(),"r");
+    if (fSteeringFile) {
+      char  line[5001];
+      while (!feof(fSteeringFile)) {
+	fgets(line,5000,fSteeringFile);
+	if ((line[0]!='#')&&(line[0]!='\n')&&(line[0]!='\r')) {
+	  char tmpString0[100]="Unset";
+	  sscanf(&line[0],"%s",tmpString0);
+	  if ( (strcmp(tmpString0,"/musr/ignore")!=0) &&(strcmp(tmpString0,"/musr/command")!=0) ) continue;
+	  
+	  
+	  char tmpString1[100]="Unset",tmpString2[100]="Unset";
+	  sscanf(&line[0],"%*s %s %s",tmpString1,tmpString2);
+	  if (strcmp(tmpString1,"materialPropertiesTable")==0){
+	    std::string materialPropertiesTableName=tmpString2;
+	    G4MaterialPropertiesTable* MPT_tmp;
+	    itMPT = materialPropertiesTableMap.find(materialPropertiesTableName);
+	    if (itMPT==materialPropertiesTableMap.end()) {  // G4MaterialPropertiesTable of this name does not exist yet --> create it  
+	      MPT_tmp = new G4MaterialPropertiesTable();
+	      materialPropertiesTableMap.insert ( std::pair<std::string,G4MaterialPropertiesTable*>(materialPropertiesTableName,MPT_tmp) );
+	    }
+	    else {
+	      MPT_tmp = itMPT->second;
+	    }
+	    
+	    char propertyName[100];
+	    int nEntries;
+	    sscanf(&line[0],"%*s %*s %*s %s %d",propertyName,&nEntries);
+	    std::cout<<" Optical Material Def:  MPT_tmp="<<MPT_tmp<<",  materialPropertiesTableName="<<materialPropertiesTableName
+		     <<",  propertyName="<<propertyName<<",  nEntries="<<nEntries<<std::endl;
+	    
+	    if (nEntries==0) {  // AddConstProperty
+	      float value;
+	      sscanf(&line[0],"%*s %*s %*s %*s %*d %g",&value);
+	      MPT_tmp->AddConstProperty(propertyName,value);
+	    }
+	    else {              // AddProperty
+	      char* pch = strstr(line,propertyName)+strlen(propertyName);
+	      float value;
+	      G4double photonEnergyArray[100];
+	      G4double valueArray[100];
+	      char dummy[100];
+	      sscanf(pch,"%s",dummy); char* pch2=strstr(pch,dummy)+strlen(dummy); pch = pch2;
+	      for (int i=0; i<nEntries; i++) {
+		sscanf(pch,"%g",&value);
+		//	      G4cout<<"  DDDDD var1="<<value<<"  &pch="<<&pch<<G4endl;
+		photonEnergyArray[i]=value;
+		sscanf(pch,"%s",dummy); char* pch2=strstr(pch,dummy)+strlen(dummy); pch = pch2;
+	      }
+	      for (int i=0; i<nEntries; i++) {
+		sscanf(pch,"%g",&value);
+		//	      G4cout<<"  DDDDD var2="<<value<<"  &pch="<<&pch<<G4endl;
+		valueArray[i]=value;
+		sscanf(pch,"%s",dummy); char* pch2=strstr(pch,dummy)+strlen(dummy); pch = pch2;
+	      }
+	      MPT_tmp->AddProperty(propertyName,photonEnergyArray,valueArray,nEntries);
+	    }
+	  }
+	  else if (strcmp(tmpString1,"setMaterialPropertiesTable")==0){
+	    char tmpString3[100]="Unset";
+	    sscanf(&line[0],"%*s %*s %*s %s",tmpString3);
+	    std::string materialPropertiesTableName=tmpString2;
+	    std::string materialName = tmpString3;
+	    G4NistManager* man = G4NistManager::Instance();
+	    G4Material* myMaterial = man->FindOrBuildMaterial(materialName);  
+
+	    G4MaterialPropertiesTable* MPT_tmp=NULL;
+	    itMPT = materialPropertiesTableMap.find(materialPropertiesTableName);
+	    if (itMPT==materialPropertiesTableMap.end()) {  // G4MaterialPropertiesTable of this name does not exist --> problem  
+	      G4cout<<"\n\n musrDetectorConstruction::DefineMaterials(): Material Properties Table \""<<materialPropertiesTableName
+		    <<"\" should be assigned to material \""<<materialName<<"\""<<G4endl;
+	      G4cout<<" but the table was not defined yet (by command /musr/command materialPropertiesTable )"<<G4endl;
+	      G4cout <<"  ===>  S T O P    F O R C E D"<<G4endl;
+	      G4cout<<line<<G4endl;
+	      exit(1);
+	    }
+	    else {
+	      MPT_tmp = itMPT->second;
+	    }
+	    myMaterial->SetMaterialPropertiesTable(MPT_tmp);
+	    G4cout<<myMaterial<<G4endl;
+	    myMaterial->GetMaterialPropertiesTable()->DumpTable();
+	  }
+	}
+      }
     }
-    G4cout<<"OK konec"<<G4endl;
-    exit(0);
+    fclose(fSteeringFile);
+
+    //    G4cout<<" 1eV = "<<eV<<G4endl;
+    //    G4cout<<" 1MeV = "<<MeV<<G4endl;
+    //    G4cout<<" 1ns = "<<ns<<G4endl;
+    //
+    //    G4NistManager* man = G4NistManager::Instance();
+    //    G4Material* scintik = man->FindOrBuildMaterial("G4_PLASTIC_SC_VINYLTOLUENE");   
+    //    G4cout<<"scintik="<<scintik<<G4endl;
+    //    if (scintik!=NULL) {
+    //      const G4int nEntries = 14;
+    //      G4double PhotonEnergy[nEntries] = 
+    //	{ 2.695*eV,   2.75489*eV, 2.8175*eV,  2.88302*eV,                 // 460, 450, 440, 430 nm
+    //	  2.95167*eV, 3.02366*eV, 3.09925*eV, 3.17872*eV, 3.26237*eV,     // 420, 410, 400, 390, 380 nm
+    //	  3.30587*eV, 3.35054*eV, 3.44361*eV, 3.542*eV,   3.64618*eV };   // 375, 370, 360, 350, 340 nm
+    //      G4double RefractiveIndex[nEntries] = 
+    //	{ 1.58, 1.58, 1.58, 1.58, 
+    //	  1.58, 1.58, 1.58, 1.58, 1.58,
+    //	  1.58, 1.58, 1.58, 1.58, 1.58  };
+    //      G4double Absorption[nEntries] =
+    //      	{ 8*cm, 8*cm, 8*cm, 8*cm, 
+    //      	  8*cm, 8*cm, 8*cm, 8*cm, 8*cm,
+    //      	  8*cm, 8*cm, 8*cm, 8*cm, 8*cm  };
+    //      G4double ScintilFast[nEntries] =
+    //	{ 0.01,  0.07, 0.15, 0.26,
+    //	  0.375, 0.52, 0.65, 0.80, 0.95,
+    //          1,     0.88, 0.44, 0.08, 0.01 };
+    //      G4double ScintilSlow[nEntries] =
+    //	{ 0.01,  0.07, 0.15, 0.26,
+    //	  0.375, 0.52, 0.65, 0.80, 0.95,
+    //          1,     0.88, 0.44, 0.08, 0.01 };
+    //      G4MaterialPropertiesTable* myMPT1 = new G4MaterialPropertiesTable();
+    //      myMPT1->AddProperty("RINDEX",        PhotonEnergy, RefractiveIndex, nEntries);
+    //      myMPT1->AddProperty("ABSLENGTH",     PhotonEnergy, Absorption,      nEntries);
+    //      myMPT1->AddProperty("FASTCOMPONENT", PhotonEnergy, ScintilFast,     nEntries);
+    //      myMPT1->AddProperty("SLOWCOMPONENT", PhotonEnergy, ScintilSlow,     nEntries);
+    //      myMPT1->AddConstProperty("SCINTILLATIONYIELD", 8400./MeV);
+    //      myMPT1->AddConstProperty("RESOLUTIONSCALE",1.0);
+    //      myMPT1->AddConstProperty("FASTTIMECONSTANT",1.6*ns);
+    //      myMPT1->AddConstProperty("SLOWTIMECONSTANT",1.6*ns);
+    //      myMPT1->AddConstProperty("YIELDRATIO",1.0);
+    //      scintik->SetMaterialPropertiesTable(myMPT1);
+
+    //      scintik->GetMaterialPropertiesTable()->DumpTable();
+    //    }
   }
 }
 
@@ -1240,6 +1503,25 @@ G4LogicalVolume* musrDetectorConstruction::FindLogicalVolume(G4String LogicalVol
   return NULL;
 }
 
+G4VPhysicalVolume* musrDetectorConstruction::FindPhysicalVolume(G4String PhysicalVolumeName) {
+  G4PhysicalVolumeStore* pPhysStore = G4PhysicalVolumeStore::GetInstance();
+  if (pPhysStore==NULL) {
+    G4cout<<"ERROR:  musrDetectorConstruction.cc:    G4PhysicalVolumeStore::GetInstance()  not found!"<<G4endl;
+  }
+  else {
+    for (unsigned int i=0; i<pPhysStore->size(); i++) {
+      G4VPhysicalVolume* pPhysVol=pPhysStore->at(i);
+      G4String iPhysName=pPhysVol->GetName(); 
+      if (iPhysName==PhysicalVolumeName) {
+	return pPhysVol;
+      }
+    }
+  }
+  G4cout<<"\n musrDetectorConstruction.cc::FindPhysicalVolume: \n          ERROR !!! Physical volume "
+	   <<PhysicalVolumeName<<" not found !!!"<<G4endl;
+  return NULL;
+}
+
 void musrDetectorConstruction::SetColourOfLogicalVolume(G4LogicalVolume* pLogVol,char* colour) {
   if (pLogVol!=NULL) {
     if      (strcmp(colour,"red"    )==0) {pLogVol->SetVisAttributes(G4Colour(1,0,0));}
diff --git a/src/musrEventAction.cc b/src/musrEventAction.cc
index 9130b9b..28cfb21 100644
--- a/src/musrEventAction.cc
+++ b/src/musrEventAction.cc
@@ -77,6 +77,8 @@ void musrEventAction::BeginOfEventAction(const G4Event* evt) {
 void musrEventAction::EndOfEventAction(const G4Event* evt)  {
   //  cout << ":." << flush;
   long thisEventNr = (long) evt->GetEventID();
+
+  //  musrSteppingAction::GetInstance()->DoAtTheEndOfEvent();
   
   // write out the root tree:
   musrRootOutput* myRootOutput = musrRootOutput::GetRootInstance();
diff --git a/src/musrPhysicsList.cc b/src/musrPhysicsList.cc
index bb31700..7dee52a 100644
--- a/src/musrPhysicsList.cc
+++ b/src/musrPhysicsList.cc
@@ -46,6 +46,8 @@
 #include "G4StepLimiter.hh"
 #include "G4UserSpecialCuts.hh"
 
+#include "G4OpticalPhysics.hh"
+
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 musrPhysicsList::musrPhysicsList():  G4VUserPhysicsList()
@@ -55,6 +57,26 @@ musrPhysicsList::musrPhysicsList():  G4VUserPhysicsList()
   cutForElectron = 0.1*mm;
   cutForMuon = 0.01*mm;
   SetVerboseLevel(0);
+
+  if (musrParameters::boolG4OpticalPhotons) {
+    
+
+    //    // * Optical Physics
+    //    G4OpticalPhysics* opticalPhysics = new G4OpticalPhysics();
+    //    RegisterPhysics( opticalPhysics );
+    //    
+    //    // adjust some parameters for the optical physics
+    //    opticalPhysics->SetWLSTimeProfile("delta");
+    //    
+    //    opticalPhysics->SetScintillationYieldFactor(1.0);
+    //    opticalPhysics->SetScintillationExcitationRatio(0.0);
+    //    
+    //    opticalPhysics->SetMaxNumPhotonsPerStep(100);
+    //    opticalPhysics->SetMaxBetaChangePerStep(10.0);
+    //    
+    //    opticalPhysics->SetTrackSecondariesFirst(true);
+  }
+
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -71,6 +93,13 @@ void musrPhysicsList::ConstructParticle()
   ConstructLeptons();
   ConstructMesons();
   ConstructBaryons();
+
+  if (musrParameters::boolG4OpticalPhotons) {
+    // optical photon
+    G4OpticalPhoton::OpticalPhotonDefinition();
+  }
+
+
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -234,6 +263,9 @@ void musrPhysicsList::ConstructProcess()
 // For a simple Muonium "scattering" when Mu hits solid materials
 #include "musrMuScatter.hh"
 
+// For optical photons
+#include "G4Scintillation.hh"
+
 // For models
 #include "G4ProcessTable.hh"
 #include "G4WentzelVIModel.hh"
@@ -273,6 +305,11 @@ void musrPhysicsList::ConstructEM()
 	G4ParticleDefinition* particleDefinition = G4ParticleTable::GetParticleTable() -> FindParticle(stringParticleName);
 	//      G4cout<<"particleDefinition of "<<stringParticleName<<" = "<<particleDefinition<<G4endl;
 	if (particleDefinition==NULL) {
+	  if ((stringParticleName=="opticalphoton")&& (!(musrParameters::boolG4OpticalPhotons))) {
+	    musrErrorMessage::GetInstance()->musrError(INFO,
+						       "musrPhysicsList: Ignoring optical photon process definition because G4OpticalPhotons in *.mac file is set to false",false);
+	    continue;
+	  }
 	  sprintf(eMessage,"musrPhysicsList:  Partile \"%s\" not found in G4ParticleTable when trying to find or assign process  \"%s\".",
 		  charParticleName,charProcessName);
 	  musrErrorMessage::GetInstance()->musrError(FATAL,eMessage,false);
@@ -294,6 +331,11 @@ void musrPhysicsList::ConstructEM()
 	  else if (stringProcessName=="G4LowEnergyRayleigh")        pManager->AddDiscreteProcess(new G4LowEnergyRayleigh);
 
 	  else if (stringProcessName=="G4CoulombScattering")        pManager->AddDiscreteProcess(new G4CoulombScattering);
+	  
+	  else if (stringProcessName=="G4OpAbsorption")             pManager->AddDiscreteProcess(new G4OpAbsorption);
+	  else if (stringProcessName=="G4OpRayleigh")               pManager->AddDiscreteProcess(new G4OpRayleigh);
+	  else if (stringProcessName=="G4OpBoundaryProcess")        pManager->AddDiscreteProcess(new G4OpBoundaryProcess);
+	  else if (stringProcessName=="G4OpWLS")                    pManager->AddDiscreteProcess(new G4OpWLS);
 	  else {
 	    sprintf(eMessage,"musrPhysicsList:  Process \"%s\" is not implemented in musrPhysicsList.cc for addDiscreteProcess.  It can be easily added.",
 		    charProcessName);
@@ -325,6 +367,7 @@ void musrPhysicsList::ConstructEM()
 	  //	  else if (stringProcessName=="G4hIonisation")        pManager->AddProcess(new G4hIonisation,nr1,nr2,nr3);
 	  //	  else if (stringProcessName=="G4hLowEnergyIonisation")      pManager->AddProcess(new G4hLowEnergyIonisation,nr1,nr2,nr3);
 	  else if (stringProcessName=="musrMuFormation")              pManager->AddProcess(new musrMuFormation,nr1,nr2,nr3);
+	  else if (stringProcessName=="G4Scintillation")              pManager->AddProcess(new G4Scintillation,nr1,nr2,nr3);
 	  // cks:  musrMuScatter could be uncommented here, but testing is needed, because Toni has some strange comments
 	  //       in his original "musrPhysicsList.cc about implementing musrMuScatter.
 	  //	  else if (stringProcessName=="musrMuScatter")       pManager->AddProcess(new musrMuScatter,nr1,nr2,nr3);
@@ -394,7 +437,6 @@ void musrPhysicsList::ConstructEM()
 	    G4MuMultipleScattering* mmm = (G4MuMultipleScattering*) processTable->FindProcess("muMsc",particleDefinition);
 	    mmm->AddEmModel(modelPriority, new G4UrbanMscModel93());
 	  }
-
 	  else {
 	    sprintf(eMessage,"musrPhysicsList:  Model \"%s\" is not implemented for \"%s\" in musrPhysicsList.cc for addModel.  It can be easily added.",
 		    charModelName,charProcessName);
diff --git a/src/musrRootOutput.cc b/src/musrRootOutput.cc
index c162522..4012aa6 100644
--- a/src/musrRootOutput.cc
+++ b/src/musrRootOutput.cc
@@ -147,6 +147,16 @@ G4bool musrRootOutput::store_fieldIntegralBz = false;
 G4bool musrRootOutput::store_fieldIntegralBz1 = false;
 G4bool musrRootOutput::store_fieldIntegralBz2 = false;
 G4bool musrRootOutput::store_fieldIntegralBz3 = false;
+G4bool musrRootOutput::store_odet_ID = true;
+G4bool musrRootOutput::store_odet_nPhot = true;
+G4bool musrRootOutput::store_odet_timeFirst = true;
+G4bool musrRootOutput::store_odet_timeA = true;
+G4bool musrRootOutput::store_odet_timeB = true;
+G4bool musrRootOutput::store_odet_timeC = true;
+G4bool musrRootOutput::store_odet_timeD = true;
+G4bool musrRootOutput::store_odet_timeE = true;
+G4bool musrRootOutput::store_odet_timeLast = true;
+
 
 G4int musrRootOutput::oldEventNumberInG4EqEMFieldWithSpinFunction=-1;
 
@@ -265,6 +275,19 @@ void musrRootOutput::BeginOfRunAction() {
     rootTree->Branch("save_poly",&save_poly,"save_poly[save_n]/D");
     rootTree->Branch("save_polz",&save_polz,"save_polz[save_n]/D");
   }
+  
+  if (store_odet_ID || store_odet_nPhot || store_odet_timeFirst || store_odet_timeA || store_odet_timeB ||
+      store_odet_timeC || store_odet_timeD || store_odet_timeE || store_odet_timeLast)
+                              {rootTree->Branch("odet_n",&odet_n,"odet_n/I");}
+  if (store_odet_ID)          {rootTree->Branch("odet_ID",&odet_ID,"odet_ID[odet_n]/I");}
+  if (store_odet_nPhot)       {rootTree->Branch("odet_nPhot",&odet_nPhot,"odet_nPhot[odet_n]/I");}
+  if (store_odet_timeFirst)   {rootTree->Branch("odet_timeFirst",&odet_timeFirst,"odet_timeFirst[odet_n]/D");}
+  if (store_odet_timeA)       {rootTree->Branch("odet_timeA",&odet_timeA,"odet_timeA[odet_n]/D");}
+  if (store_odet_timeB)       {rootTree->Branch("odet_timeB",&odet_timeB,"odet_timeB[odet_n]/D");}
+  if (store_odet_timeC)       {rootTree->Branch("odet_timeC",&odet_timeC,"odet_timeC[odet_n]/D");}
+  if (store_odet_timeD)       {rootTree->Branch("odet_timeD",&odet_timeD,"odet_timeD[odet_n]/D");}
+  if (store_odet_timeE)       {rootTree->Branch("odet_timeE",&odet_timeE,"odet_timeE[odet_n]/D");}
+  if (store_odet_timeLast)       {rootTree->Branch("odet_timeLast",&odet_timeLast,"odet_timeLast[odet_n]/D");}
 
   //  htest1 = new TH1F("htest1","The debugging histogram 1",50,-4.,4.);
   //  htest2 = new TH1F("htest2","The debugging histogram 2",50,0.,3.142);
@@ -317,7 +340,7 @@ void musrRootOutput::FillEvent() {
   htest5->Fill(atan2(posIniMomy,posIniMomx));
   htest6->Fill(atan2(sqrt(posIniMomx*posIniMomx+posIniMomy*posIniMomy),posIniMomz));
   if (weight>0.) {
-    if ( !((musrParameters::storeOnlyEventsWithHits)&&(det_n<=0)) ) {
+    if ( !((musrParameters::storeOnlyEventsWithHits)&&(det_n<=0)&&(odet_n<=0)) ) {
       rootTree->Fill();
     }
   }
@@ -349,7 +372,7 @@ void musrRootOutput::ClearAllRootVariables() {
   BzIntegral1 = -1000; BzIntegral2 = -1000; BzIntegral3 = -1000; 
   det_n=0;
   save_n=0;
-
+  odet_n=0;
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -498,3 +521,27 @@ void musrRootOutput::SetDetectorInfoVvv (G4int nDetectors,
     det_VvvParticleID[nDetectors]=particleID;
   }
 }
+
+
+void musrRootOutput::SetOPSAinfo    (G4int nDetectors, G4int ID, G4int nPhot, G4double timeFirst, G4double timeA, 
+				     G4double timeB, G4double timeC, G4double timeD, G4double timeE, G4double timeLast)
+{
+  if ((nDetectors<0)||(nDetectors>=(odet_nMax-1))) {
+    char message[200];  
+    sprintf(message,"musrRootOutput.cc::SetOPSAInfo: nDetectors %i is larger than det_nMax = %i",nDetectors,det_nMax);
+    musrErrorMessage::GetInstance()->musrError(SERIOUS,message,false);
+    return;
+  }
+  else {
+    odet_n=nDetectors+1; 
+    odet_ID[nDetectors]=ID;
+    odet_nPhot[nDetectors]=nPhot;
+    odet_timeFirst[nDetectors]=timeFirst/microsecond;
+    odet_timeA[nDetectors]=timeA/microsecond;
+    odet_timeB[nDetectors]=timeB/microsecond;
+    odet_timeC[nDetectors]=timeC/microsecond;
+    odet_timeD[nDetectors]=timeD/microsecond;
+    odet_timeE[nDetectors]=timeE/microsecond;
+    odet_timeLast[nDetectors]=timeLast/microsecond;
+  }
+}
diff --git a/src/musrScintSD.cc b/src/musrScintSD.cc
index 5cd5697..7fe042b 100644
--- a/src/musrScintSD.cc
+++ b/src/musrScintSD.cc
@@ -29,9 +29,14 @@
 #include "G4ios.hh"
 #include <algorithm>   // needed for the sort() function
 #include "G4VProcess.hh"  // needed for the degugging message of the process name
+#include "G4OpBoundaryProcess.hh"      // Optical photon process
+#include "G4RunManager.hh"
 #include "musrParameters.hh"
 #include "musrErrorMessage.hh"
 #include "musrSteppingAction.hh"
+//#include "TCanvas.h"
+#include "TMath.h"
+#include "TF1.h"
 #include <vector>
 
 //bool myREMOVEfunction (int i,int j) { return (i<j); }
@@ -49,17 +54,38 @@
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
+Double_t poissonf(Double_t* x, Double_t* par) {
+  return par[0]*TMath::Poisson(x[0],par[1]);
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+
 musrScintSD::musrScintSD(G4String name)
 :G4VSensitiveDetector(name)
 {
+  pointer=this;
   G4String HCname;
   collectionName.insert(HCname="scintCollection");
+  OPSA_minNrOfDetectedPhotons = 1;
+  OPSA_signalSeparationTime   = 10.;
+  OPSA_fracA = 0.01;
+  OPSA_fracB = 0.05;
+  OPSA_fracC = 0.10;
+  OPSA_fracD = 0.20;
+  OPSA_fracE = 0.5;
+  bool_multimapOfEventIDsForOPSAhistosEXISTS = false;
+  OPSAhistoNbin = 100;
+  OPSAhistoMin =0;
+  OPSAhistoMax = 10.;
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 musrScintSD::~musrScintSD(){ }
 
+musrScintSD* musrScintSD::pointer=NULL;
+musrScintSD* musrScintSD::GetInstance() {return pointer;}
+
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 void musrScintSD::Initialize(G4HCofThisEvent* HCE) {
@@ -78,6 +104,7 @@ void musrScintSD::Initialize(G4HCofThisEvent* HCE) {
   myStoreOnlyEventsWithHitInDetID = musrParameters::storeOnlyEventsWithHitInDetID;
   musrSteppingAction* myMusrSteppingAction = musrSteppingAction::GetInstance();
   boolIsVvvInfoRequested = myMusrSteppingAction->IsVvvInfoRequested();
+  myRootOutput = musrRootOutput::GetRootInstance();
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -91,7 +118,15 @@ G4bool musrScintSD::ProcessHits(G4Step* aStep,G4TouchableHistory*)
   }
 
   G4Track* aTrack = aStep->GetTrack();
-  G4String actualVolume=aTrack->GetVolume()->GetLogicalVolume()->GetName();
+  //  G4LogicalVolume* hitLogicalVolume = aTrack->GetVolume()->GetLogicalVolume();
+  //  G4String actualVolume=aTrack->GetVolume()->GetLogicalVolume()->GetName();
+  G4String hitLogicalVolumeName = aTrack->GetVolume()->GetLogicalVolume()->GetName();
+  G4String particleName=aTrack->GetDefinition()->GetParticleName();
+  //  if (particleName=="opticalphoton") {G4cout<<"UFON JE TU:  edep="<<edep<<G4endl; return false;}
+  if (particleName=="opticalphoton") {
+    ProcessOpticalPhoton(aStep);
+    return false;
+  }
 
   // If requested, store only the hit that happened first (usefull for some special studies, not for a serious simulation)
   if (myStoreOnlyTheFirstTimeHit) {
@@ -101,17 +136,21 @@ G4bool musrScintSD::ProcessHits(G4Step* aStep,G4TouchableHistory*)
       return false;
     }
   }
+  
+  
 
   musrScintHit* newHit = new musrScintHit();
-  newHit->SetParticleName (aTrack->GetDefinition()->GetParticleName());
+  //  newHit->SetParticleName (aTrack->GetDefinition()->GetParticleName());
+  newHit->SetParticleName(particleName);
   G4int particleID = aTrack->GetDefinition()->GetPDGEncoding();
   newHit->SetParticleID (particleID);
   newHit->SetEdep     (edep);
   newHit->SetPrePos   (aStep->GetPreStepPoint()->GetPosition());
   newHit->SetPostPos  (aStep->GetPostStepPoint()->GetPosition());
   newHit->SetPol      (aTrack->GetPolarization());
-  G4LogicalVolume* hitLogicalVolume = aTrack->GetVolume()->GetLogicalVolume();
-  newHit->SetLogVolName(hitLogicalVolume->GetName());
+  //  G4LogicalVolume* hitLogicalVolume = aTrack->GetVolume()->GetLogicalVolume();
+  newHit->SetLogVolName(hitLogicalVolumeName);
+  //  newHit->SetLogVolName(hitLogicalVolume->GetName());
   newHit->SetGlobTime(aTrack->GetGlobalTime());
   //  Warning - aStep->IsFirstStepInVolume() only available in Geant version >= 4.8.2 !
   //  newHit->SetFirstStepInVolumeFlag(aStep->IsFirstStepInVolume());
@@ -141,6 +180,94 @@ G4bool musrScintSD::ProcessHits(G4Step* aStep,G4TouchableHistory*)
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
+void musrScintSD::ProcessOpticalPhoton(G4Step* aStep) {
+  //  //Was the photon absorbed by the absorption process ?
+  //  const G4VProcess* process = aStep->GetPostStepPoint()->GetProcessDefinedStep();
+  //  G4String processName =  (process) ? process->GetProcessName() : "Unknown";
+  //  G4cout<<"processName="<<processName<<G4endl;
+  //  if (processName=="OpAbsorption") {
+  //    G4cout<<"\n ABSORPTION\n";
+  //  }
+
+
+  G4OpBoundaryProcessStatus boundaryStatus=Undefined;
+  static G4OpBoundaryProcess* boundaryProc=NULL;
+  //find the boundary process only once
+  if(!boundaryProc){
+    G4ProcessManager* pm = aStep->GetTrack()->GetDefinition()->GetProcessManager();
+    G4int nprocesses = pm->GetProcessListLength();
+    G4ProcessVector* pv = pm->GetProcessList();
+    for(G4int i=0; i<nprocesses; i++){
+      if((*pv)[i]->GetProcessName()=="OpBoundary"){
+	boundaryProc = (G4OpBoundaryProcess*)(*pv)[i];
+	break;
+      }
+    }
+  }
+  
+  boundaryStatus=boundaryProc->GetStatus();
+  
+  //Check to see if the partcile was actually at a boundary
+  //Otherwise the boundary status may not be valid
+  //Prior to Geant4.6.0-p1 this would not have been enough to check
+  if(aStep->GetPostStepPoint()->GetStepStatus()==fGeomBoundary){
+    //    G4cout<<"      boundaryStatus="<<boundaryStatus<<"   ";
+    G4String actualVolume = aStep->GetTrack()->GetVolume()->GetLogicalVolume()->GetName();
+    Int_t detID =  myRootOutput->ConvertVolumeToID(actualVolume);
+    //    G4cout<<" detID ="<<detID<<"  actualVolume="<<actualVolume;
+    optHitMapType::iterator iter = optHitMap.find(detID);
+    if (iter==optHitMap.end()) {  // optHitDetectorMapType does not exist ==> create it
+      optHitDetectorMapType* optHitDetectorMapTmp = new optHitDetectorMapType;
+      optHitMap.insert(std::pair<Int_t,optHitDetectorMapType*>(detID,optHitDetectorMapTmp));
+      iter = optHitMap.find(detID);
+    }
+    optHitDetectorMapType* optHitDetectorMap = (*iter).second;
+    //	  optHitDetectorMapType optHitDetectorMap = optHitMap[detID];
+    G4double tmpTime = aStep->GetPreStepPoint()->GetGlobalTime();
+    //    G4cout<<"   tmpTime="<<tmpTime;
+    //	  optHitDetectorMap->insert(std::pair<G4double,G4int>(tmpTime,boundaryStatus));
+    
+    if (boundaryStatus!=Detection) {
+      char message[200];
+      sprintf(message,"musrScintSD.cc::ProcessOpticalPhoton(): Optical photon boundary status is not Detection but %i",boundaryStatus);
+      musrErrorMessage::GetInstance()->musrError(FATAL,message,false);
+    }
+
+    //    switch(boundaryStatus){
+    //    case Absorption:
+    //      G4cout<<" AAAAAAAAAAAAAAAAAA  Absorption"<<G4endl;
+    //      break;
+    //    case Detection: //Note, this assumes that the volume causing detection
+      //is the photocathode because it is the only one with
+      //non-zero efficiency
+    optHitDetectorMap->insert(std::pair<G4double,G4int>(tmpTime,boundaryStatus));
+      //      G4cout<<"   Detection"<<G4endl;
+      //      break;
+      //    case FresnelReflection:
+      //    case LambertianReflection:
+      //    case LobeReflection:
+      //    case TotalInternalReflection:
+      //    case SpikeReflection:
+      //      G4cout<<"   Reflection"<<G4endl;
+      //      break;
+      //    case StepTooSmall:
+      //      G4cout<<"   StepTooSmall"<<G4endl;
+      //      break;
+      //    case NoRINDEX:
+      //      G4cout<<"   NoRINDEX"<<G4endl;
+      //      break;
+      //      //	  case :
+      //      //	    G4cout<<"   "<<G4endl;
+      //      //	    break;
+      //    default:
+      //      G4cout<<"   SomethingElse"<<G4endl;
+      //      break;
+      //    }
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
 void musrScintSD::EndOfEvent(G4HCofThisEvent*) {
   if (verboseLevel>1) { 
     G4cout<<"VERBOSE 2:  musrScintSD::EndOfEvent"<<G4endl;
@@ -149,8 +276,6 @@ void musrScintSD::EndOfEvent(G4HCofThisEvent*) {
 	   << " hits in the scint chambers: " << G4endl;
   } 
   
-  //  Positron_momentum_already_stored=0;
-  musrRootOutput* myRootOutput = musrRootOutput::GetRootInstance();
 
   G4int NbHits = scintCollection->entries(); 
   
@@ -367,6 +492,144 @@ void musrScintSD::EndOfEvent(G4HCofThisEvent*) {
 
     }
   }   //end "if (NbHits>0)"
+
+  // Analyse optical photons if they were produced
+  if (musrParameters::boolG4OpticalPhotons) EndOfEvent_OptiacalPhotons();
 }  
 
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void musrScintSD::EndOfEvent_OptiacalPhotons() {
+  //  G4double  kCarTolerance = G4GeometryTolerance::GetInstance()  ->GetSurfaceTolerance();
+  //  G4cout<<"  DEBUG 10:  kCarTolerance="<<kCarTolerance<<G4endl;
+
+  if (optHitMap.empty()) return;
+  
+  G4RunManager* fRunManager = G4RunManager::GetRunManager();
+  G4int eeeventID = fRunManager->GetCurrentEvent()->GetEventID();
+
+  for (optHitMapType::const_iterator it=optHitMap.begin() ; it != optHitMap.end(); it++ ) {
+    G4bool boolStoreThisOPSAhist = false;
+    G4int  OPSA_detID= it->first;
+    optHitDetectorMapType* optHitDetectorMap = it->second;
+
+    // Check whether OPSA histograming of times of optical photon detection is required for this eventID.
+    if (bool_multimapOfEventIDsForOPSAhistosEXISTS) {
+      if (multimapOfEventIDsForOPSAhistos.find(eeeventID)!=multimapOfEventIDsForOPSAhistos.end()) {
+	//  Now check whether the histogramming is required for the currently analysed detector.
+	std::pair<multimapOfEventIDsForOPSAhistos_Type::iterator,multimapOfEventIDsForOPSAhistos_Type::iterator> retOPSAhist;
+	multimapOfEventIDsForOPSAhistos_Type::iterator itOPSAhist;
+	retOPSAhist = multimapOfEventIDsForOPSAhistos.equal_range(eeeventID);
+	for (itOPSAhist = retOPSAhist.first; itOPSAhist!=retOPSAhist.second; itOPSAhist++) {
+	  //  Store histograms if the second parameters of eventsForOPSAhistos (i.e. detector ID) 
+          //      is set to 0 or corresponds to the currently analysed sensitive detector.
+	  if ( (itOPSAhist->second == 0) || (itOPSAhist->second == OPSA_detID) ) boolStoreThisOPSAhist = true;
+	}
+      }
+    }
+
+    if (optHitDetectorMap->empty()) continue;
+    optHitDetectorMapType::const_iterator it2_START = optHitDetectorMap->begin();
+    optHitDetectorMapType::const_iterator it2_STOP  = it2_START;
+    optHitDetectorMapType::const_iterator it2_LAST = optHitDetectorMap->end();  it2_LAST--;
+    Double_t time = -1000, lastTime = -1000;
+    G4int    OPSA_nPhot       = 0;
+    G4double OPSA_timeFirst   = -1000;
+    G4double OPSA_timeA       = -1000;
+    G4double OPSA_timeB       = -1000;
+    G4double OPSA_timeC       = -1000;
+    G4double OPSA_timeD       = -1000;
+    G4double OPSA_timeE       = -1000;
+    G4double OPSA_timeLast    = -1000;
+    G4int    iHistNr          = -1;
+    TF1* poiss = NULL;
+    for (optHitDetectorMapType::const_iterator it2 = optHitDetectorMap->begin(); it2 != optHitDetectorMap->end(); it2++ ) {
+      OPSA_nPhot++;
+      lastTime = time;
+      time = it2->first;
+      if (OPSA_nPhot==1) lastTime=time;  // First photon does not have a proper "lastTime" defined
+      if ( (it2==it2_LAST) || ((time-lastTime) > OPSA_signalSeparationTime)) {
+	// The iterator it2 reached last element of the map optHitDetectorMap       or
+	// the time difference between two subsequently detected photons is too big ==> divide the signal into more signals
+	it2_STOP = it2;
+	if (it2==it2_LAST) it2_STOP++;          // if we are at the end of optHitDetectorMap, 
+                                                //    the it2_STOP should point just behind the last map element
+	else OPSA_nPhot--;                      // if we split the optHitDetectorMap, however, the latest optical photon
+	                                        //    should have not be counted, because it belong to the next signal.
+	optHitDetectorMapType optHitDetectorSUBmap(it2_START, it2_STOP);
+	it2_START = it2_STOP;
+	if (OPSA_nPhot >= OPSA_minNrOfDetectedPhotons) {   // ignore hits with too low number of detected photons
+	  G4double OPSA_f_nPhot = OPSA_nPhot;
+	  G4int NA = int (OPSA_fracA * OPSA_f_nPhot + 0.5);
+	  G4int NB = int (OPSA_fracB * OPSA_f_nPhot + 0.5);
+	  G4int NC = int (OPSA_fracC * OPSA_f_nPhot + 0.5);
+	  G4int ND = int (OPSA_fracD * OPSA_f_nPhot + 0.5);
+	  G4int NE = int (OPSA_fracE * OPSA_f_nPhot + 0.5);
+	  Int_t nP=0;
+	  // Define OPSA histograms if required for this event
+	  if (boolStoreThisOPSAhist) {
+	    iHistNr++;
+	    char nameHist[100];  sprintf(nameHist,"OPSAhist_%d_%d_%d",eeeventID,OPSA_detID,iHistNr);
+	    char nameHistTitle[100]; sprintf(nameHistTitle,"OPSAhist_%d_%d_%d;time (ns);Nr of photons",eeeventID,OPSA_detID,iHistNr);
+	    OPSAhisto = new TH1D(nameHist, nameHistTitle, OPSAhistoNbin, OPSAhistoMin, OPSAhistoMax);
+	    poiss = new TF1("poiss",poissonf,0.,.5,2); // x in [0;300], 2
+	    poiss->SetParameter(0,1);
+	    poiss->SetParameter(1,1);
+	  }
+	  
+	  for (optHitDetectorMapType::const_iterator it3 = optHitDetectorSUBmap.begin(); it3 != optHitDetectorSUBmap.end(); it3++) {
+	    nP++;
+	    if (nP==1)          OPSA_timeFirst = it3->first;
+	    if (nP==NA)         OPSA_timeA     = it3->first;
+	    if (nP==NB)         OPSA_timeB     = it3->first;
+	    if (nP==NC)         OPSA_timeC     = it3->first;
+	    if (nP==ND)         OPSA_timeD     = it3->first;
+	    if (nP==NE)         OPSA_timeE     = it3->first;
+	    if (nP==OPSA_nPhot) OPSA_timeLast  = it3->first;
+	    if (boolStoreThisOPSAhist) OPSAhisto->Fill((it3->first)-OPSA_timeFirst+0.00000000001);
+	  }
+	  signalInfo* mySignalInfo = new signalInfo(OPSA_detID,OPSA_nPhot,OPSA_timeFirst,OPSA_timeA,OPSA_timeB,OPSA_timeC,
+						    OPSA_timeD,OPSA_timeE,OPSA_timeLast);
+	  OPSA_signal_Map.insert(std::pair<G4int,signalInfo*>(OPSA_nPhot,mySignalInfo) );
+	  if (boolStoreThisOPSAhist) {
+	    OPSAhisto->Fit(poiss,"Q");
+	    OPSAhisto->Write();
+	    //	    TCanvas* cc = new TCanvas();
+	    //	    OPSAhisto->Draw();
+	    //	    cc->Update();
+	  }
+	}
+	OPSA_nPhot       = 0;
+	OPSA_timeFirst   = -1000;
+	OPSA_timeA       = -1000;
+	OPSA_timeB       = -1000;
+	OPSA_timeC       = -1000;
+	OPSA_timeD       = -1000;
+	OPSA_timeE       = -1000;
+	OPSA_timeLast    = -1000;
+      }
+    }
+  }
+
+  // Now delete all optHitDetectorMap*
+  for (optHitMapType::const_iterator it=optHitMap.begin() ; it != optHitMap.end(); it++ ) {
+    delete (it->second);
+  }
+  optHitMap.clear();
+
+
+  //  Now store the information of all mySignalInfo*  to rootOutputFile
+  G4int nn=0;
+  for (OPSA_signal_MapType::reverse_iterator ritOPSA = OPSA_signal_Map.rbegin(); ritOPSA != OPSA_signal_Map.rend(); ritOPSA++) {
+    (ritOPSA->second) -> transferDataToRoot(myRootOutput,nn);
+    nn++;
+  }
+  //  Now delete all mySignalInfo*  from  OPSA_signal_Map
+  for (OPSA_signal_MapType::const_iterator itOPSA = OPSA_signal_Map.begin(); itOPSA != OPSA_signal_Map.end(); itOPSA++) {
+    delete (itOPSA->second);
+  }
+  OPSA_signal_Map.clear();
+
+}
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/src/musrSteppingAction.cc b/src/musrSteppingAction.cc
index aae160c..4ee6742 100644
--- a/src/musrSteppingAction.cc
+++ b/src/musrSteppingAction.cc
@@ -29,25 +29,35 @@
 #include "G4MagneticField.hh"          // needed for storing the magnetic field to the Root class
 #include "G4FieldManager.hh"           // ---------------||------------------
 #include "G4TransportationManager.hh"  // ---------------||------------------
+//#include "G4OpBoundaryProcess.hh"      // Optical photon process
 #include "musrErrorMessage.hh"
 #include "musrParameters.hh"
 #include "F04GlobalField.hh"
+//#include "TCanvas.h"
+//#include "TMath.h"
+//#include "TF1.h"
+//#include "G4GeometryTolerance.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//Double_t poissonf(Double_t* x, Double_t* par) {
+//  return par[0]*TMath::Poisson(x[0],par[1]);
+//}
+
+
 
 musrSteppingAction::musrSteppingAction()  { 
- pointer=this;
-
- boolIsAnySpecialSaveVolumeDefined = false;
- boolIsVvvInfoRequested = false;
- boolMuonEventReweighting = false;
- boolCalculateFieldIntegral = false;
- myRootOutput = musrRootOutput::GetRootInstance();
- if (myRootOutput == NULL) {
-   musrErrorMessage::GetInstance()->musrError(FATAL,
-    	  "musrSteppingAction::musrSteppingAction():  pointer to the musrRootOutput class not found! ==> EXECUTION STOPPED",true);
- }
- lastActualVolume="Unset";
+  pointer=this;
+  
+  boolIsAnySpecialSaveVolumeDefined = false;
+  boolIsVvvInfoRequested = false;
+  boolMuonEventReweighting = false;
+  boolCalculateFieldIntegral = false;
+  myRootOutput = musrRootOutput::GetRootInstance();
+  if (myRootOutput == NULL) {
+    musrErrorMessage::GetInstance()->musrError(FATAL,
+					       "musrSteppingAction::musrSteppingAction():  pointer to the musrRootOutput class not found! ==> EXECUTION STOPPED",true);
+  }
+  lastActualVolume="Unset";
 }
 
 musrSteppingAction::~musrSteppingAction()  { 
@@ -140,7 +150,7 @@ void musrSteppingAction::UserSteppingAction(const G4Step* aStep)  {
     }
     
 
-    // Save info about the old tracks, if the user wish to have Vvv info in the output Root Tree.
+    // Save info about the old tracks, if the user wishes to have Vvv info in the output Root Tree.
     if (boolIsVvvInfoRequested) {
       G4VPhysicalVolume* nextVolume = aTrack->GetNextVolume();
       if (nextVolume!=NULL) {
@@ -375,6 +385,7 @@ void musrSteppingAction::UserSteppingAction(const G4Step* aStep)  {
   }
 }
 
+
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 
@@ -472,3 +483,17 @@ G4bool  musrSteppingAction::AreTracksCommingFromSameParent(G4int trackID1, G4int
   //  G4cout<<"\t\t\t\ttrack1="<<track1<<"\ttrack2="<<track2<<G4endl;
   return false;
 }
+
+//  //Double_t musrSteppingAction::poissonf(Double_t* x, Double_t* par)                                         
+//  //{ 
+//Double_t musrSteppingAction::poissonf(Double_t* x, Double_t* par) {
+//  //  if (x<0)
+//  //    return 0;
+//  //  else if (x == 0.0)
+//  //    return 1./Math::Exp(par);
+//  //  else {
+//  //    Double_t lnpoisson = x*log(par)-par-LnGamma(x+1.);
+//  //    return Exp(lnpoisson);
+//  //  }
+//  return par[0]*TMath::Poisson(x[0],par[1]);
+//}