28.9.2011 Kamil Sedlak

1) implemented a new command /gun/turtleMomentumScalingFactor in the musrSim
   (see the manual)
2) some small changes in the musrSimAna

doc/musrSim.tex

@@ -820,6 +820,13 @@ This way the smearing of the vertex as well as beam tilt/pitch are propagated th
 	{\it turtleSmearingFactor} = 0 ~~... the muon beam momentum will be set to the constant value of {\it turtleMomentumP0}.\\
 	{\it turtleSmearingFactor} = 200 ... the muon beam will have two times broader distribution compared to the original TURTLE file.
 
+\item{\bf /gun/turtleMomentumScalingFactor \emph{scalingFactor} }\\
+        Modify the turtle momentum by this scaling factor -- e.g.\ if the user wants to shift
+	the mean momentum from 28\,MeV/c to 14\,MeV/c, this scaling factor has to be set to 0.5.
+	Note that also the momentum bite will be reduced by the same scaling factor.
+	This scaling is done before the command ``/gun/turtleMomentumBite'' (if both are
+	requested).
+
 \item{\bf /gun/turtleFirstEventNr \emph{lineNumberOfTurtleFile} }\\
 	Set the line number that should be taken as the first event from the TURTLE input file.
 	This option is needed when the user wants to reproduce the simulation of an event

include/musrPrimaryGeneratorAction.hh

@@ -83,6 +83,7 @@ class musrPrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
     void SetOrReadTheRandomNumberSeeds(G4Event* anEvent);
     void SetTurtleMomentumBite (G4ThreeVector smearingParam)
                                {turtleMomentumBite=true; turtleMomentumP0=smearingParam[0]*MeV; turtleSmearingFactor=smearingParam[1]*0.01;}
+    void SetTurtleMomentumScalingFactor(G4double momentumScaling) {turtleMomentumScalingFactor=momentumScaling;}
     void SetPrimaryParticule(G4String particleName);
     static G4String GetPrimaryName();                
 
@@ -127,6 +128,7 @@ class musrPrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
     G4bool   turtleMomentumBite;
     G4double turtleMomentumP0;
     G4double turtleSmearingFactor;
+    G4double turtleMomentumScalingFactor;
     void     swapTheAxisInTurtle(float& x_x, float& x_xprime, float& y_y, float& y_yprime);
 
 public: 

include/musrPrimaryGeneratorMessenger.hh

@@ -72,6 +72,7 @@ class musrPrimaryGeneratorMessenger: public G4UImessenger
     G4UIcmdWithADoubleAndUnit* setTurtleZ0Cmd;
     G4UIcmdWithAString*        setTurtleInterpretAxesCmd;
     G4UIcmdWith3Vector*        setTurtleMomentumBite;
+    G4UIcmdWithADouble*        setTurtleMomentumScalingFactor;
     G4UIcmdWithAnInteger*      setTurtleEventNrCmd;
 };
 #endif

musrSimAna/musrAnalysis.cxx

@@ -272,6 +272,9 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
 	    hInfo->Fill(21, rot_ref_frequency);  // value may be overwritten - just the last rot. ref. frame will be saved
 	    hInfo->Fill(22, rot_ref_phase);      // value may be overwritten - just the last rot. ref. frame will be saved
 	  }
+	  else if (strcmp(furtherOption,"correctexpdecay")==0) {
+	    myTH->SetExpDecayCorrection();
+	  }
 	}
 	else {
 	  Double_t* pointerToVariable2 = iter2->second;
@@ -587,6 +590,21 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
 	  funct -> SetParameter(0,p0);
 	  funct -> SetParameter(1,p1);
 	}
+	else if (strcmp(functionName,"TFieldCosPLUSbg")==0) {
+	  funct = new TF1("TFieldCosPLUSbg","[3]*(1+[2]*cos(x*[0]+[1]))+[4]*exp(x/2.19703)");
+	  funct -> SetParameter(0,p0);
+	  funct -> SetParameter(1,p1);
+	  funct -> SetParameter(2,p2);
+	  funct -> SetParameter(3,p3);
+	  funct -> SetParameter(4,p4);
+	}
+	else if (strcmp(functionName,"TFieldCos")==0) {
+	  funct = new TF1("TFieldCos","[3]*(1+[2]*cos(x*[0]+[1]))");
+	  funct -> SetParameter(0,p0);
+	  funct -> SetParameter(1,p1);
+	  funct -> SetParameter(2,p2);
+	  funct -> SetParameter(3,p3);
+	}
 	else if (strcmp(functionName,"rotFrameTime20")==0) {
 	  //	  funct = new TF1("rotFrameTime20","[2]*exp(-x/2.19703)*cos(x*[0]+[1]) ");
 	  funct = new TF1("rotFrameTime20","[2]*cos(x*[0]+[1]) ");

musrSimAna/musrTH.cxx

@@ -5,15 +5,16 @@
 musrTH::musrTH(char* dimension, char* histoName, char* histoTitle, Int_t nrOfBinsX, Double_t minBoundaryX, Double_t maxBoundaryX, Int_t nrOfBinsY, Double_t minBoundaryY, Double_t maxBoundaryY, Double_t* variableName1, Double_t* variableName2) {
 //musrTH::musrTH(char* dimension, char* histoName, char* histoTitle, Int_t nrOfBins, Double_t minBoundary, Double_t maxBoundary, Double_t* variableName1, Double_t* varibleName2) {
   std::cout<<" Defining "<<dimension<<" histogram "<<histoName<<" \""<<histoTitle<<"\" "
-           <<nrOfBinsX<<" "<<minBoundaryX<<" "<<maxBoundaryX
+           <<nrOfBinsX<<" "<<minBoundaryX<<" "<<maxBoundaryX<<" "
-	   <<nrOfBinsY<<" "<<minBoundaryY<<" "<<maxBoundaryY
+	   <<nrOfBinsY<<" "<<minBoundaryY<<" "<<maxBoundaryY<<" "
-	   <<" "<<variableName1<<" "<<variableName2<<std::endl;
+	   <<variableName1<<" "<<variableName2<<std::endl;
 
   strcpy(histogramName,histoName);
   variableToBeFilled_X = variableName1;
   variableToBeFilled_Y = variableName2;
   funct = NULL;
   bool_rotating_reference_frame=false;
+  bool_exp_decay_correction=false;
   rot_ref_frequency=0;
   rot_ref_phase=0;
   strcpy(funct_option,"");
@@ -49,6 +50,11 @@ void musrTH::FillTH1D(Double_t vaha, Bool_t* cond){
     if (bool_rotating_reference_frame) {
       //      Double_t var = *variableToBeFilled_X;
       Double_t waha = vaha*exp((*variableToBeFilled_X)/2.19703)*cos(2*musrAnalysis::pi*rot_ref_frequency*(*variableToBeFilled_X)+rot_ref_phase);
+      //      std::cout<<"rot_ref_frequency="<<rot_ref_frequency<<std::endl;
+      if (cond[i])  histArray1D[i]->Fill(*variableToBeFilled_X,waha);
+    }
+    else if (bool_exp_decay_correction) {
+      Double_t waha = vaha*exp((*variableToBeFilled_X)/2.19703);
       if (cond[i])  histArray1D[i]->Fill(*variableToBeFilled_X,waha);
     }
     else {
@@ -219,16 +225,18 @@ void  musrTH::ListHistograms() {
 //==============================================================================================
 void  musrTH::FitHistogramsIfRequired(Double_t omega) {
   if (funct==NULL) return;
-  if (bool_rotating_reference_frame) omega = fabs(omega) - 2*musrAnalysis::pi*rot_ref_frequency;
+  if (bool_rotating_reference_frame) omega = fabs(omega) - 2*musrAnalysis::pi*fabs(rot_ref_frequency);
   std::cout<<"============================================================================================================"<<std::endl;
   std::cout<<"Fitting \""<<funct->GetName()<<"\",   funct_xMin="<<funct_xMin<<"   funct_xMax="<<funct_xMax<<"    omega="<<omega<<std::endl;
   if (strcmp(funct->GetName(),"funct1")==0) {funct->FixParameter(0,omega);}
   if (strcmp(funct->GetName(),"funct2")==0) {funct->FixParameter(0,omega);}
   if (strcmp(funct->GetName(),"funct3")==0) {funct->FixParameter(0,omega);}
   if (strcmp(funct->GetName(),"funct4")==0) {funct->FixParameter(0,omega);}
+  if (strcmp(funct->GetName(),"TFieldCosPLUSbg")==0) {funct->FixParameter(0,omega);}
+  if (strcmp(funct->GetName(),"TFieldCos")==0) {funct->FixParameter(0,omega);}
   if (strcmp(funct->GetName(),"rotFrameTime20")==0)  {
     if (funct->GetParameter(0)==0) {
-      funct->SetParameter(0,omega); std::cout<<" FUNKCE rotFrameTime20"<<"omega initialy at "<<omega<<std::endl;
+      funct->SetParameter(0,omega); std::cout<<" FUNKCE rotFrameTime20, omega initialy at "<<omega<<std::endl;
     }
   }
 

musrSimAna/musrTH.hh

@@ -43,6 +43,7 @@ public:
   void  FitHistogramsIfRequired(Double_t omega);
   void  SetRotatingReferenceFrame(Double_t frequency, Double_t phase) {bool_rotating_reference_frame=true; 
                                                                        rot_ref_frequency=frequency; rot_ref_phase=phase;}
+  void  SetExpDecayCorrection() {bool_exp_decay_correction=true;}
   void  ListHistograms();
 
 private:
@@ -57,6 +58,7 @@ private:
   char   funct_option[100];
   Bool_t bool_rotating_reference_frame;
   Double_t rot_ref_frequency, rot_ref_phase;
+  Bool_t bool_exp_decay_correction;
   //  Double_t N0_FromLastFit;
   //  Double_t BinWidth_FromLastFit;
   //  Double_t funct_p0, funct_p1, funct_p2, funct_p3, funct_p4, funct_p5, funct_p6, funct_p7, funct_p8, funct_p9;

src/musrPrimaryGeneratorAction.cc

@@ -71,7 +71,7 @@ musrPrimaryGeneratorAction::musrPrimaryGeneratorAction(
    muonDecayTimeMin(-1), muonDecayTimeMax(-1), muonMeanLife(2197.03*ns),
    takeMuonsFromTurtleFile(false), z0_InitialTurtle(0),
    numberOfGeneratedEvents(0), 
-   turtleMomentumBite(false), turtleMomentumP0(0.), turtleSmearingFactor(0.)
+   turtleMomentumBite(false), turtleMomentumP0(0.), turtleSmearingFactor(0.), turtleMomentumScalingFactor(1.)
 						      //, firstCall(true)
 {
   //create a messenger for this class
@@ -175,7 +175,7 @@ void musrPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
       yangle = yAngleTmp*mrad;
       x      = xTmp*cm + (z0-z0_InitialTurtle)*tan(xangle) ;       // usually z0 is negative
       y      = yTmp*cm + (z0-z0_InitialTurtle)*tan(yangle) ;       // z0_InitialTurtle is the z0 at whith the turtle file was generated.
-      p      = pTmp*GeV;
+      p      = pTmp*GeV*turtleMomentumScalingFactor;
       // add some offset, if requested:
       x      = x + x0;
       y      = y + y0;

src/musrPrimaryGeneratorMessenger.cc

@@ -171,6 +171,12 @@ musrPrimaryGeneratorMessenger::musrPrimaryGeneratorMessenger(musrPrimaryGenerato
   setTurtleMomentumBite->SetGuidance("   broader beam.  The third parameter is dummy.");
   setTurtleMomentumBite->SetParameterName("mes_turtleMomentumP0","mes_turtleSmearingFactor","mes_turtleSmearingDummy",true,true);
 
+  setTurtleMomentumScalingFactor = new G4UIcmdWithADouble("/gun/turtleMomentumScalingFactor",this);
+  setTurtleMomentumScalingFactor->SetGuidance(" Modify the turtle momentum by this scaling factor - e.g. if the user wants to shift");
+  setTurtleMomentumScalingFactor->SetGuidance("   the mean momentum from 28 MeV/c to 14 MeV/c, this scaling factor has to be set to 0.5 .");
+  setTurtleMomentumScalingFactor->SetGuidance("   Note that also the momentum bite will be reduced by the same scaling factor.");
+  setTurtleMomentumScalingFactor->SetParameterName("mes_turtleMomentumScalingFactor",true);
+
 
   setTurtleEventNrCmd = new G4UIcmdWithAnInteger("/gun/turtleFirstEventNr",this);
   setTurtleEventNrCmd->SetGuidance("Set the line number that should be taken as the first event from the turtle input file.");
@@ -208,6 +214,7 @@ musrPrimaryGeneratorMessenger::~musrPrimaryGeneratorMessenger()
   delete setTurtleZ0Cmd;
   delete setTurtleInterpretAxesCmd;
   delete setTurtleMomentumBite;
+  delete setTurtleMomentumScalingFactor;
   delete setTurtleEventNrCmd;
 }
 
@@ -265,6 +272,8 @@ void musrPrimaryGeneratorMessenger::SetNewValue(G4UIcommand * command,G4String n
     { musrAction->SetTurtleInterpretAxes(newValue);}
   if( command == setTurtleMomentumBite)
     { musrAction-> SetTurtleMomentumBite(setTurtleMomentumBite->GetNew3VectorValue(newValue)); }
+  if( command == setTurtleMomentumScalingFactor)
+    { musrAction-> SetTurtleMomentumScalingFactor(setTurtleMomentumScalingFactor->GetNewDoubleValue(newValue));}
   if( command == setTurtleEventNrCmd)
     { musrAction->SetTurtleInputFileToEventNo(setTurtleEventNrCmd->GetNewIntValue(newValue));}
 }