15.9.2011 Kamil Sedlak

1) Bug removed from the simulation of optical photons with
    the option "APDcellsEffectRequested" -- namely, when any
    simulated photon hit a given APD cell, the cell became "blind"
    to any other photon, irrespective of the arrival time of the
    next photon. However, we need that if the "second" photon
    (simulated at later phase) enters the APD cell before the
    "first" one, that the time of photon arrival changes to the
    earlierst photon.
 2) A new keyword "sampleID ..." (defining which volume is the
    sample volume) was added to the musrSimAna.  The manual
    was updated accordingly.

doc/musrSimAna.tex

@@ -435,6 +435,11 @@ All events should/have to be (?) saved in the Root tree
      This defines a group of detectors, where \emph{group} stands for ``B'' (backward), 
      ``F'' (forward), ``U'' (up), ``D'' (down), ``L'' (left) and ``R'' (right) detectors.
      This grouping is used in the definition of some conditions.
+   \item{\bf sampleID \emph{ID$_1$ ID$_2$ \ldots ID$_n$} } \\
+     Defines which volume (or volumes, if there are more) is the sample.  Typically, the
+     sample is just one volume, but sometimes there is a smaller volume inside of the sample,
+     to which a field is applied, so the small volume also has to be considered as the sample.
+     This information is needed for the condition ``muonDecayedInSample\_gen''.
    \item{\bf  setSpecialAnticoincidenceTimeWindow  \emph{detectorID}  \emph{timeMin}  \emph{timeMax} \emph{unit}} \\
      This command sets a special anti-coincidence time window for a detector \emph{detectorID}.
      Normally, the anti-coincidence time window is defined by {\tt VCOINCIDENCEW}, and is the same for all anti-coincidence

musrSimAna/musrAnalysis.cxx

@@ -474,6 +474,23 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
 	  nscan = sscanf(pch,"%d",&N1);
 	} while (nscan==1);
       }
+      else if (strcmp(tmpString0,"sampleID")==0) {
+	int nscan; int N1; char NAME[100];
+	char *pch = line + strlen("sampleID");
+	do {
+	  nscan = sscanf(pch,"%d",&N1);
+	  SampleDetIDList.push_back(N1);
+	  nscan = sscanf(pch,"%s",NAME);
+	  char* pch2 = strstr(pch ,NAME)+strlen(NAME);
+	  pch=pch2;
+	  nscan = sscanf(pch,"%d",&N1);
+	} while (nscan==1);
+	std::cout<<"SampleID :";
+	for (std::list <Int_t>::const_iterator it=SampleDetIDList.begin(); it!=SampleDetIDList.end(); ++it) {
+	  std::cout<<" "<< *it;
+	}
+	std::cout<<std::endl;
+      }
       else if  (strcmp(tmpString0,"setSpecialAnticoincidenceTimeWindow")==0) {
 	int channelNumber;
 	double min, max;
@@ -823,10 +840,12 @@ void musrAnalysis::SaveHistograms(char* runChar, char* v1190FileName) {
 
   Long64_t numberOfMuonCandidates = mCounter->GetNumberOfMuonCandidates();
   Long64_t numberOfMuonCandidatesAfterVK = mCounter->GetNumberOfMuonCandidatesAfterVK();
+  Long64_t numberOfMuonCandidatesAfterVKandDoubleHitRemoval = mCounter->GetNumberOfMuonCandidatesAfterVKandDoubleHitRemoval();
   Double_t durationOfExperiment = (numberOfRewinds*rewindTimeBins*tdcresolution+currentTime)/1000000;
   hInfo->Fill(5,(Double_t) numberOfMuonCandidates);  // number of "triggers"
   hInfo->Fill(4,(Double_t) numberOfMuonCandidatesAfterVK);  // number of "triggers"
   hInfo->Fill(3,durationOfExperiment);
+  hInfo->Fill(8,(Double_t) numberOfMuonCandidatesAfterVKandDoubleHitRemoval);
   //==============================
   //  Write out the histograms
   std::cout<<"musrAnalysis::SaveHistograms()"<<std::endl;
@@ -1086,7 +1105,9 @@ void musrAnalysis::FillHistograms(Int_t iiiEntry) {
     
     // CALCULATE CONDITIONS
     alwaysTrue = true;
-    muonDecayedInSample_gen = muDecayDetID==201;
+    // Find whether the muon decayed in the sample;  the sample may consist of more than just one volume.  The ID
+    // of the sample volume (or volumes) is defined in the *.v1190 file after the keyword "sampleID":
+    muonDecayedInSample_gen = (find(SampleDetIDList.begin(), SampleDetIDList.end(), muDecayDetID)) != SampleDetIDList.end() ;
     muonTriggered_gen = muM0Time > -1000;
     muonTriggered_det = mCounterHitExistsForThisEventID;
     positronHit_det   = pCounterHitExistsForThisEventID;

musrSimAna/musrAnalysis.hh

@@ -409,6 +409,7 @@ private:
   std::list <Int_t> D_posCounterList;
   std::list <Int_t> L_posCounterList;
   std::list <Int_t> R_posCounterList;
+  std::list <Int_t> SampleDetIDList;
   //  std::list <Int_t>::iterator posCounterList_Iterator;
 };
 

musrSimAna/musrCounter.cxx

@@ -32,6 +32,7 @@ musrCounter::musrCounter(int CHANNEL_NR, char CHANNEL_NAME[200], char CHANNEL_TY
   doubleHitN=0;
   numberOfMuonCandidates=0;
   numberOfMuonCandidatesAfterVK=0;
+  numberOfMuonCandidatesAfterVKandDoubleHitRemoval=0;
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -194,6 +195,8 @@ Bool_t musrCounter::GetNextGoodMuon(Int_t evtID, Long64_t timeBinMin, Long64_t&
       goto endOfThisHit;   
     }
 
+    numberOfMuonCandidatesAfterVKandDoubleHitRemoval++;
+  
     kEntry   = (it->second)->eventEntry;
     idet     = (it->second)->det_i;
     idetID   = (it->second)->det_id;

musrSimAna/musrCounter.hh

@@ -64,7 +64,7 @@ class musrCounter {
     void myPrintThisCounter(Int_t evtID,  Int_t detail=2);
     Long64_t GetNumberOfMuonCandidates(){return numberOfMuonCandidates;}
     Long64_t GetNumberOfMuonCandidatesAfterVK(){return numberOfMuonCandidatesAfterVK;}
-
+    Long64_t GetNumberOfMuonCandidatesAfterVKandDoubleHitRemoval(){return numberOfMuonCandidatesAfterVKandDoubleHitRemoval;}
 
   private:
   //    static musrCounter* pointerToAnalysis;
@@ -91,6 +91,7 @@ class musrCounter {
   Int_t doubleHitN;
   Long64_t numberOfMuonCandidates;
   Long64_t numberOfMuonCandidatesAfterVK;
+  Long64_t numberOfMuonCandidatesAfterVKandDoubleHitRemoval;
 };
 
 #endif

src/musrScintSD.cc

@@ -269,18 +269,25 @@ void musrScintSD::ProcessOpticalPhoton(G4Step* aStep) {
     // whether a given cell already fired - if so, do not store the photon.
     G4bool APDcellAlreadyFired = false;
     G4int APDcellID = 0;
+    if (APDcellsTimeVariationRequested)  tmpTime += G4RandGauss::shoot(0,APDcellsTimeVariationSigma);
+
     if (APDcellsEffectRequested) {
       APDcellID = FindAPDcellID(aStep);
       //      G4cout<<"Cell ID="<<APDcellID<<G4endl;
-      for (optHitDetectorMapType::const_iterator it2 = optHitDetectorMap->begin(); it2 != optHitDetectorMap->end(); it2++ ) {
+      for (optHitDetectorMapType::iterator it2 = optHitDetectorMap->begin(); it2 != optHitDetectorMap->end(); it2++ ) {
 	if ((it2->second)==APDcellID) {
+	  //  G4cout<<"Already fired cell ="<<APDcellID<<G4endl;
 	  APDcellAlreadyFired = true;
-	  //	  G4cout<<"HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH already fired cell ="<<APDcellID<<G4endl;
+	  //  Now lets check whether the new photon came earlier than the already saved one.
+	  //  If so, delete the previous photon and store the time of the new one instead:
+	  if (tmpTime<(it2->first)) {
+	    optHitDetectorMap->erase(it2);
+	    optHitDetectorMap->insert(std::pair<G4double,G4int>(tmpTime,APDcellID));
+	  }
 	  break;
 	}
       }
     }
-    if (APDcellsTimeVariationRequested)  tmpTime += G4RandGauss::shoot(0,APDcellsTimeVariationSigma);
     if (!APDcellAlreadyFired) optHitDetectorMap->insert(std::pair<G4double,G4int>(tmpTime,APDcellID));
     //    G4cout<<"  tmpTime="<<tmpTime<<G4endl;
       //      G4cout<<"   Detection"<<G4endl;