6.9.2010 - Kamil Sedlak

1) bug found and corrected in the field-map extrapolation of the electric 
     and magnetic fields for the spin-rotator (3DEOpera and 3D and 
     symmetryType=1 or 2). The documentation was updated.
  2) Some new variables added to the musrSimAna - the documentation has not
     been updated yet.

doc/musrSim.tex

@@ -474,19 +474,23 @@ care.
 			$x$, $y$ and $z$).  In such cases, the user can specify this in the field map
 			file using the keyword ``symmetryType \emph{number}'', where the \emph{number}
 			specifies how the field should be extrapolated to other octants.  
-			The ``symmetryType 1'' case means that the planes of symmetry are (x,y) and
+			The ``symmetryType 1'' case can be used for a field between two electrodes
-			(x,z), i.e. if the field at point $(x,y,z)$ is $(F_x,F_y,F_z)$, the field in
+			(e.g.\ of a spin rotator), where the electrodes are parallel with the $(y,z)$ plane
-			different octants will look like this: 
+			(i.e.\ the main component of the field is along the $x$ axis).
-			$(-x,y,z) \rightarrow (F_x,F_y,F_z)$;
+			In this case:  if the field at point $(x,y,z)$ is $(F_x,F_y,F_z)$, the field in
-			$(x,-y,z) \rightarrow (F_x,-F_y,F_z)$;
+			different octants will look be: \\
-			$(-x,-y,z) \rightarrow (F_x,-F_y,F_z)$;
+			$(x,y,z) \rightarrow (F_x,F_y,F_z)$;\\
-			$(x,y,-z) \rightarrow (F_x,F_y,-F_z)$;
+			$(x,y,-z) \rightarrow (F_x,F_y,-F_z)$;\\
-			$(-x,y,-z) \rightarrow (F_x,F_y,-F_z)$;
+			$(x,-y,z) \rightarrow (F_x,-F_y,F_z)$;\\
-			$(x,-y,-z) \rightarrow (F_x,-F_y,-F_z)$;
+			$(x,-y,-z) \rightarrow (F_x,-F_y,-F_z)$;\\
-			$(-x,-y,-z) \rightarrow (F_x,-F_y,-F_z)$. \\
+			$(-x,y,z) \rightarrow (F_x,-F_y,-F_z)$;\\
-			Similar case is the ``symmetryType 2'', where the planes of symmetry are
+			$(-x,y,-z) \rightarrow (F_x,-F_y,F_z)$;\\
-			(x,y) and (y,z).
+			$(-x,-y,z) \rightarrow (F_x,F_y,-F_z)$;\\
-			These two symmetry types are realised in a the spin rotator oriented along the z axis.\\
+			$(-x,-y,-z) \rightarrow (F_x,F_y,F_z)$. \\
+			Similar case is the ``symmetryType 2'', where the main field component is along the
+			$y$ axis.
+			These two symmetry types are realised as electric and magnetic fields 
+			in a spin rotator oriented along the z axis.\\
 			Example of the beginning of the 3DBOpera-type field map file:\\
 			2   2   55\\
 			1 X\\

musrSimAna/musrAnalysis.cxx

@@ -985,6 +985,18 @@ 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;
+    detP_x = -1001.;
+    detP_y = -1001.;
+    detP_z = -1001.;
+    detP_time_start = -1001.;
+    detP_time_end   = -1001.;
+    detP_theta = -1001.;
+    detP_phi = -1001.;
+    detP_phi_MINUS_pos_Phi = -1001.;
+    detP_phi_MINUS_pos_Phi360 = -1001.;
+    detP_theta_MINUS_pos_Theta = -1001.;
+    detP_theta_MINUS_pos_Theta360 = -1001.;
+    detP_time_start_MINUS_muDecayTime = -1001.;
     pileup_eventID             = -1001;
     pileup_muDecayDetID_double = -1001;
     pileup_muDecayPosX         = -1000000000;
@@ -997,18 +1009,31 @@ void musrAnalysis::FillHistograms(Int_t iiiEntry) {
       Long64_t dataBinMax = timeBin0+dataWindowBinMax;
       pCounterHitExistsForThisEventID = PositronCounterHit(eventID,dataBinMin,dataBinMax,timeBin1,timeBin2,posEntry,idetP,idetP_ID,idetP_edep,doubleHitP);
       //cDEL      if (pCounterHitExistsForThisEventID) std::cout<<"        timeBin1-timeBin2 = "<<timeBin1<<"-"<<timeBin2<<"="<<timeBin1-timeBin2 <<std::endl;
-      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)
+      if (pCounterHitExistsForThisEventID&&(posEntry>0)) {
+	// Get variables of the detector hit corresponding to the positron candidate:
 	fChain->GetEntry(posEntry);
-	pileup_eventID             = eventID;
+	Double_t detP_edep = det_edep[idetP];
-	pileup_muDecayDetID_double = muDecayDetID;
+	if (detP_edep!=idetP_edep) {std::cout<<"KIKS: detP_edep="<<detP_edep<<";   idetP_edep="<<idetP_edep<<";    idetP= "<<idetP<<std::endl; exit(1);}
-	pileup_muDecayPosX         = muDecayPosX;
+	detP_x = det_x[idetP];
-	pileup_muDecayPosY         = muDecayPosY;
+	detP_y = det_y[idetP];
-	pileup_muDecayPosZ         = muDecayPosZ;
+	detP_z = det_z[idetP];
-	pileup_muDecayPosR         = sqrt(muDecayPosX*muDecayPosX+muDecayPosY*muDecayPosY);
+	detP_time_start = det_time_start[idetP];
-	//	if (pileup_muDecayDetID_double==-1000) {
+	detP_time_end   = det_time_end[idetP];
-	//	  std::cout<<"DEBUG:  pileup_muDecayDetID_double==-1000, posEntry="<<posEntry<<", eventID="<<eventID<<", idetP_edep="<<idetP_edep<<", idetP="<<idetP<<",  idetP_ID="<<idetP_ID<<std::endl;
+	detP_phi = myAtan2(detP_y,detP_x) * 180./pi;
-	//	}
+	detP_theta = myAtan2(sqrt(detP_y*detP_y+detP_x+detP_x),detP_z) * 180./pi;
+	//
+	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_muDecayPosX         = muDecayPosX;
+	  pileup_muDecayPosY         = muDecayPosY;
+	  pileup_muDecayPosZ         = muDecayPosZ;
+	  pileup_muDecayPosR         = sqrt(muDecayPosX*muDecayPosX+muDecayPosY*muDecayPosY);
+	  //	fChain->GetEntry(iiiEntry);
+	}
 	fChain->GetEntry(iiiEntry);
       }
     }
@@ -1024,12 +1049,22 @@ void musrAnalysis::FillHistograms(Int_t iiiEntry) {
     muTargetPol_Theta360 = (muTargetPol_Theta<0) ? muTargetPol_Theta+360. : muTargetPol_Theta;
     muTargetPol_Phi   = myAtan2(muTargetPolY,muTargetPolX) * 180./pi;
     muTargetPol_Phi360=  (muTargetPol_Phi<0) ? muTargetPol_Phi+360. : muTargetPol_Phi;
+    muDecayPol_Theta = myAtan2(sqrt(muDecayPolX*muDecayPolX+muDecayPolY*muDecayPolY),muDecayPolZ) * 180./pi;
+    muDecayPol_Theta360 = (muDecayPol_Theta<0) ? muDecayPol_Theta+360. : muDecayPol_Theta;
+    muDecayPol_Phi   = myAtan2(muDecayPolY,muDecayPolX) * 180./pi;
+    muDecayPol_Phi360=  (muDecayPol_Phi<0) ? muDecayPol_Phi+360. : muDecayPol_Phi;
     // Initial positron
     pos_Trans_Momentum = sqrt(posIniMomX*posIniMomX+posIniMomY*posIniMomY);
     pos_Momentum = sqrt(pos_Trans_Momentum*pos_Trans_Momentum+posIniMomZ*posIniMomZ);
     pos_Radius = pos_Trans_Momentum/(-BFieldAtDecay_Bz)/0.3;
-    pos_Theta = myAtan2(pos_Trans_Momentum,posIniMomZ);
+    pos_Theta = myAtan2(pos_Trans_Momentum,posIniMomZ) * 180./pi;
-    pos_Phi   = myAtan2(posIniMomY,posIniMomX);
+    pos_Theta360 = (pos_Theta<0) ? pos_Theta+360. : pos_Theta;
+    pos_Phi      = myAtan2(posIniMomY,posIniMomX) * 180./pi;
+    pos_Phi360   = (pos_Phi<0) ? pos_Phi+360. : pos_Phi;
+    pos_Theta_MINUS_muDecayPol_Theta = deltaAngle(pos_Theta,muDecayPol_Theta);
+    pos_Theta_MINUS_muDecayPol_Theta360 = (pos_Theta_MINUS_muDecayPol_Theta<0) ? pos_Theta_MINUS_muDecayPol_Theta+360 : pos_Theta_MINUS_muDecayPol_Theta;
+    pos_Phi_MINUS_muDecayPol_Phi = deltaAngle(pos_Phi,muDecayPol_Phi);
+    pos_Phi_MINUS_muDecayPol_Phi360 = (pos_Phi_MINUS_muDecayPol_Phi<0) ? pos_Phi_MINUS_muDecayPol_Phi+360 : pos_Phi_MINUS_muDecayPol_Phi;
     // Detector info
     det_m0edep = (mCounterHitExistsForThisEventID) ? idetM_edep : -1000.;
     //   det_time0 = timeBin0*tdcresolution;
@@ -1040,6 +1075,13 @@ void musrAnalysis::FillHistograms(Int_t iiiEntry) {
     gen_time10 = muDecayTime-muM0Time;
     det_time10_MINUS_gen_time10 = (det_time10 - gen_time10)/picosecond;
     det_posEdep = (pCounterHitExistsForThisEventID) ?  idetP_edep : -1000.;
+    det_time1_MINUS_muDecayTime = (timeBin1*tdcresolution-muDecayTime)/picosecond;
+    if ((detP_time_start>-998.)&&(muDecayTime>-998.)) detP_time_start_MINUS_muDecayTime = (detP_time_start - muDecayTime)/picosecond;
+    detP_phi_MINUS_pos_Phi = deltaAngle(detP_phi,pos_Phi);
+    //  std::cout<<"detP_phi_MINUS_pos_Phi="<<detP_phi_MINUS_pos_Phi<<std::endl;
+    detP_phi_MINUS_pos_Phi360 = (detP_phi_MINUS_pos_Phi<0) ? detP_phi_MINUS_pos_Phi+360 : detP_phi_MINUS_pos_Phi;
+    detP_theta_MINUS_pos_Theta = deltaAngle(detP_theta,pos_Theta);
+    detP_theta_MINUS_pos_Theta360 = (detP_theta_MINUS_pos_Theta<0) ? detP_theta_MINUS_pos_Theta+360 : detP_theta_MINUS_pos_Theta;
     //    std::cout<<eventID<<"  det_time10="<<det_time10<<" t1="<<timeBin1*tdcresolution<<" t0="<<timeBin0*tdcresolution<<"  gen_time10="<<gen_time10<<std::endl;
     
     // CALCULATE CONDITIONS
@@ -1347,3 +1389,15 @@ Double_t musrAnalysis::myAtan2(Double_t y, Double_t x) {
   if ( (x==0) && (y==0) ) return -1000.;
   return atan2(y,x);
 }
+//================================================================
+Double_t musrAnalysis::deltaAngle(Double_t alpha, Double_t beta) {
+  // Calculates the difference between angle alpha and beta.
+  // The angles alpha and beta are in degrees.
+  // The difference will be in the range of (-180,180> degrees.
+  if ((alpha<-998.)||(beta<-998.)) {return -1001.;}   // one of the angles was undefined;
+  Double_t delta = alpha - beta;
+  if (delta<=-180) {delta+=360;}
+  else {if (delta>180)  delta-=360;}
+  return delta;
+}
+//================================================================

musrSimAna/musrAnalysis.hh

@@ -199,11 +199,21 @@ public :
    Double_t        muTargetPol_Theta360;
    Double_t        muTargetPol_Phi;
    Double_t        muTargetPol_Phi360;
+   Double_t        muDecayPol_Theta;
+   Double_t        muDecayPol_Theta360;
+   Double_t        muDecayPol_Phi;
+   Double_t        muDecayPol_Phi360;
    Double_t        pos_Trans_Momentum;
    Double_t        pos_Momentum;
    Double_t        pos_Radius;
    Double_t        pos_Theta;
+   Double_t        pos_Theta360;
    Double_t        pos_Phi;
+   Double_t        pos_Phi360;
+   Double_t        pos_Theta_MINUS_muDecayPol_Theta;
+   Double_t        pos_Theta_MINUS_muDecayPol_Theta360;
+   Double_t        pos_Phi_MINUS_muDecayPol_Phi;
+   Double_t        pos_Phi_MINUS_muDecayPol_Phi360;
   //   Double_t        det_time0;
   //   Double_t        get_time0;
   //   Double_t        det_time1;
@@ -212,6 +222,19 @@ public :
    Double_t        gen_time10;
    Double_t        det_time10_MINUS_gen_time10;
    Double_t        det_time20;
+   Double_t        det_time1_MINUS_muDecayTime;
+   Double_t        detP_x;
+   Double_t        detP_y;
+   Double_t        detP_z;
+   Double_t        detP_time_start;
+   Double_t        detP_time_end;
+   Double_t        detP_theta;
+   Double_t        detP_phi;
+   Double_t        detP_phi_MINUS_pos_Phi;
+   Double_t        detP_phi_MINUS_pos_Phi360;
+   Double_t        detP_theta_MINUS_pos_Theta;
+   Double_t        detP_theta_MINUS_pos_Theta360;
+   Double_t        detP_time_start_MINUS_muDecayTime;
    Double_t        pileup_eventID;
    Double_t        pileup_muDecayDetID_double;
    Double_t        pileup_muDecayPosX;
@@ -279,6 +302,7 @@ public :
    void MyPrintTree();
    void MyPrintConditions();
    Double_t myAtan2(Double_t y, Double_t x);
+   Double_t deltaAngle(Double_t alpha, Double_t beta);
 
    typedef std::map<int,Double_t> phaseShiftMapType;
    phaseShiftMapType phaseShiftMap;
@@ -452,11 +476,21 @@ musrAnalysis::musrAnalysis(TTree *tree)
   variableMap["muTargetPol_Theta360"]=&muTargetPol_Theta360;
   variableMap["muTargetPol_Phi"]=&muTargetPol_Phi;
   variableMap["muTargetPol_Phi360"]=&muTargetPol_Phi360;
+  variableMap["muDecayPol_Theta"]=&muDecayPol_Theta;
+  variableMap["muDecayPol_Theta360"]=&muDecayPol_Theta360;
+  variableMap["muDecayPol_Phi"]=&muDecayPol_Phi;
+  variableMap["muDecayPol_Phi360"]=&muDecayPol_Phi360;
   variableMap["pos_Trans_Momentum"]=&pos_Trans_Momentum;
   variableMap["pos_Momentum"]=&pos_Momentum;
   variableMap["pos_Radius"]=&pos_Radius;
   variableMap["pos_Theta"]=&pos_Theta;
+  variableMap["pos_Theta360"]=&pos_Theta360;
   variableMap["pos_Phi"]=&pos_Phi;
+  variableMap["pos_Phi360"]=&pos_Phi360;
+  variableMap["pos_Theta_MINUS_muDecayPol_Theta"]=&pos_Theta_MINUS_muDecayPol_Theta;
+  variableMap["pos_Theta_MINUS_muDecayPol_Theta360"]=&pos_Theta_MINUS_muDecayPol_Theta360;
+  variableMap["pos_Phi_MINUS_muDecayPol_Phi"]=&pos_Phi_MINUS_muDecayPol_Phi;
+  variableMap["pos_Phi_MINUS_muDecayPol_Phi360"]=&pos_Phi_MINUS_muDecayPol_Phi360;
   //  variableMap["det_time0"]=&det_time0;
   //  variableMap["gen_time0"]=&gen_time0;
   //  variableMap["det_time1"]=&det_time1;
@@ -464,6 +498,19 @@ 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["det_time1_MINUS_muDecayTime"]=&det_time1_MINUS_muDecayTime;
+  variableMap["detP_x"]=&detP_x;
+  variableMap["detP_y"]=&detP_x;
+  variableMap["detP_z"]=&detP_x;
+  variableMap["detP_time_start"]=&detP_time_start;
+  variableMap["detP_time_end"]=&detP_time_end;
+  variableMap["detP_theta"]=&detP_theta;
+  variableMap["detP_phi"]=&detP_phi;
+  variableMap["detP_phi_MINUS_pos_Phi"]=&detP_phi_MINUS_pos_Phi;
+  variableMap["detP_phi_MINUS_pos_Phi360"]=&detP_phi_MINUS_pos_Phi360;
+  variableMap["detP_theta_MINUS_pos_Theta"]=&detP_theta_MINUS_pos_Theta;
+  variableMap["detP_theta_MINUS_pos_Theta360"]=&detP_theta_MINUS_pos_Theta360;
+  variableMap["detP_time_start_MINUS_muDecayTime"]=&detP_time_start_MINUS_muDecayTime;
   variableMap["pileup_eventID"]=&pileup_eventID;
   variableMap["pileup_muDecayDetID"]=&pileup_muDecayDetID_double;
   variableMap["pileup_muDecayPosX"]=&pileup_muDecayPosX;

src/musrPrimaryGeneratorAction.cc

@@ -366,9 +366,13 @@ void musrPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
     //    G4double decaytime       = -muonMeanLife*log(1-randomVal);
     //
     //  The following code is numerically more stable compared to the commented lines above:
-    G4double expMin    = exp(-muonDecayTimeMin/muonMeanLife);
+    G4double decaytime;
-    G4double expMax    = exp(-muonDecayTimeMax/muonMeanLife);
+    if (muonDecayTimeMin==muonDecayTimeMax) {decaytime=muonDecayTimeMin;}
-    G4double decaytime = -muonMeanLife * log(G4UniformRand()*(expMax-expMin)+expMin);
+    else {
+      G4double expMin    = exp(-muonDecayTimeMin/muonMeanLife);
+      G4double expMax    = exp(-muonDecayTimeMax/muonMeanLife);
+      decaytime = -muonMeanLife * log(G4UniformRand()*(expMax-expMin)+expMin);
+    }
     //    G4cout<<"decaytime="<<decaytime/ns<<"ns."<< G4endl;
     generatedMuon->SetProperTime(decaytime);
   }

src/musrTabulatedElementField.cc

@@ -570,10 +570,12 @@ void musrTabulatedElementField::addFieldValue3D(const G4double point[4],
     if (symmetryType!=0) {
       switch(symmetryType) {
       case (1): 
-	B[1]*=y_sign;  B[2]*=z_sign;
+	//	B[1]*=y_sign;  B[2]*=z_sign;
+	B[1]*=y_sign*x_sign; B[2]*=z_sign*x_sign;
 	break;
       case (2): 
-	B[0]*=x_sign;  B[2]*=z_sign;
+	//	B[0]*=x_sign;  B[2]*=z_sign;
+	B[0]*=x_sign*y_sign;  B[2]*=z_sign*y_sign;
 	break;
       case (101):
 	if (xySwap) {