26.8.2011 Kamil Sedlak

- implementataion of volume "TubeWithHolePlusTubeHole"
 - some other changes in musrSimAna

doc/musrSimAna.tex

@@ -334,6 +334,12 @@ All events should/have to be (?) saved in the Root tree
 	   \item[wght]                \ldots weight of the event.
 	   \item[det\_m0edep]         \ldots energy deposited in the M-counter that gives the muon signal.
 	   \item[det\_posEdep]        \ldots energy deposited in the P-counter that gives the positron signal.
+	   \item[muIniPosR]           \ldots $\sqrt{ {\rm muIniPosX}^2 + {\rm muIniPosY}^2}$.
+	   \item[muIniMomTrans]       \ldots $\sqrt{ {\rm muIniMomX}^2 + {\rm muIniMomY}^2}$.
+	   \item[muTargetPol\_Theta]   \ldots theta angle of the muon spin when muon enters target (-180,180 deg).
+	   \item[muTargetPol\_Theta360]\ldots theta angle of the muon spin when muon enters target (0,360 deg).
+	   \item[muTargetPol\_Phi]    \ldots phi angle of the muon spin when muon enters target (-180,180 deg).
+	   \item[muTargetPol\_Phi360] \ldots phi angle of the muon spin when muon enters target (0,360 deg).
 	   \item[pos\_Momentum]       \ldots magnitude of the momentum of the decay positron (``generated'', not measurable variable).
 	   \item[pos\_Trans\_Momentum] \ldots transverse momentum of the decay positron.
 	   \item[pos\_Radius]         \ldots positron radius calculated from the decay positron momentum and magnetic 

musrSimAna/musrAnalysis.cxx

@@ -575,12 +575,20 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
 	  funct -> SetParameter(1,p1);
 	  funct -> SetParameter(2,p2);
 	}
+	else if (strcmp(functionName,"gaus")==0) {
+	  std::cout<<"Gausssssssss"<<std::endl;
+	  funct = new TF1("gaus","gaus");
+	  funct -> SetParameter(0,p0);
+	  funct -> SetParameter(1,p1);
+	  funct -> SetParameter(2,p2);
+	  std::cout<<"GausssssssssGausssssssss"<<std::endl;
+	}
 	else {
 	  std::cout<<"musrAnalysis::ReadInInputParameters:  function \""<<functionName<<"\" not defined! ==> S T O P"<<std::endl;
 	  exit(1);
 	}
 	
-	musrTH* tmpHistograms;
+	musrTH* tmpHistograms=NULL;
 	for(std::list<musrTH*>::const_iterator it = listOfAllHistograms1D.begin(); it != listOfAllHistograms1D.end(); ++it) {
 	  if ((*it)->IsThisHistoNamed(histoName)) {
 	    tmpHistograms = *it;
@@ -1005,6 +1013,13 @@ void musrAnalysis::FillHistograms(Int_t iiiEntry) {
     //    }
 
     // CALCULATE VARIABLES
+    // Initial muon
+    muIniPosR = sqrt(muIniPosX*muIniPosX+muIniPosY*muIniPosY);
+    muIniMomTrans = sqrt(muIniMomX*muIniMomX+muIniMomY*muIniMomY);
+    muTargetPol_Theta = myAtan2(sqrt(muTargetPolX*muTargetPolX+muTargetPolY*muTargetPolY),muTargetPolZ) * 180./pi;
+    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;
     // Initial positron
     pos_Trans_Momentum = sqrt(posIniMomX*posIniMomX+posIniMomY*posIniMomY);
     pos_Momentum = sqrt(pos_Trans_Momentum*pos_Trans_Momentum+posIniMomZ*posIniMomZ);

musrSimAna/musrAnalysis.hh

@@ -193,6 +193,12 @@ public :
    Double_t        wght;
    Double_t        det_m0edep;
    Double_t        det_posEdep;
+   Double_t        muIniPosR;
+   Double_t        muIniMomTrans;
+   Double_t        muTargetPol_Theta;
+   Double_t        muTargetPol_Theta360;
+   Double_t        muTargetPol_Phi;
+   Double_t        muTargetPol_Phi360;
    Double_t        pos_Trans_Momentum;
    Double_t        pos_Momentum;
    Double_t        pos_Radius;
@@ -440,6 +446,12 @@ musrAnalysis::musrAnalysis(TTree *tree)
   variableMap["wght"]=&wght;
   variableMap["det_m0edep"]=&det_m0edep;
   variableMap["det_posEdep"]=&det_posEdep;
+  variableMap["muIniPosR"]=&muIniPosR;
+  variableMap["muIniMomTrans"]=&muIniMomTrans;
+  variableMap["muTargetPol_Theta"]=&muTargetPol_Theta;
+  variableMap["muTargetPol_Theta360"]=&muTargetPol_Theta360;
+  variableMap["muTargetPol_Phi"]=&muTargetPol_Phi;
+  variableMap["muTargetPol_Phi360"]=&muTargetPol_Phi360;
   variableMap["pos_Trans_Momentum"]=&pos_Trans_Momentum;
   variableMap["pos_Momentum"]=&pos_Momentum;
   variableMap["pos_Radius"]=&pos_Radius;

musrSimAna/musrTH.cxx

@@ -234,7 +234,8 @@ void  musrTH::FitHistogramsIfRequired(Double_t omega) {
 
   Double_t ppp[100];
   std::cout<<"    Initial parameter setting: ";
-  Int_t n_ppp = funct->GetNumberFreeParameters();
+  //  Int_t n_ppp = funct->GetNumberFreeParameters();
+  Int_t n_ppp = funct->GetNpar();
   for (Int_t i=0; i<n_ppp; i++) {ppp[i]=funct->GetParameter(i); std::cout<<ppp[i]<<", ";}
   std::cout<<std::endl;
       

src/musrDetectorConstruction.cc

@@ -330,6 +330,37 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	  G4Transform3D transform(rot,zTrans);
 	  solid =  new G4SubtractionSolid(solidName, solidOuterDetTube, solidInnerDetTube, transform);
 	}
+	else if (strcmp(tmpString2,"TubeWithHolePlusTubeHole")==0) {
+	  // Create a tube with a partial hole inside it:        ----------
+          //                                                     |        |
+          //                                                     |     _8_| 7
+          //                                                     |    | 
+          //                                                     |    |___  
+          //                                                     |        | 6
+          //                                                     |        |
+          //                                                     ----------
+          //                                                            
+          //                                                     ----------
+          //                                                     |        |
+          //                                                     |     ___|
+          //                                                     |    |
+          //                                                     |    |___
+          //                                                     |        |
+          //                                                     ----------
+	  // First 5 parameters as for the outer tube, the 6th, 7th and 8th define the second hole.
+	  sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %lf %lf %lf %lf %s",
+		 name,&x1,&x2,&x3,&x4,&x5,&x6,&x7,&x8,material);
+	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*g %*g %*g %*s %lf %lf %lf %s %s %s %d %s",
+		 &posx,&posy,&posz,mothersName,rotMatrix,sensitiveDet,&volumeID,actualFieldName);
+	  solidName+=name;
+	  G4double roundingErr=0.01*mm;  // to avoid some displaying problems of the subtracted volumes
+	  G4Tubs* solidInnerDetTube = new G4Tubs("SolidInnerDetTube",x6*mm-roundingErr,x7*mm+roundingErr,x8/2*mm+roundingErr,x4*deg,x5*deg);
+	  G4Tubs* solidOuterDetTube = new G4Tubs("SolidOuterDetTube",x1*mm,x2*mm,x3*mm,x4*deg,x5*deg);
+	  G4RotationMatrix rot(0,0,0);
+	  G4ThreeVector zTrans(0,0,(x8/2.-x3)*mm);
+	  G4Transform3D transform(rot,zTrans);
+	  solid =  new G4SubtractionSolid(solidName, solidOuterDetTube, solidInnerDetTube, transform);
+	}
 	else if (strcmp(tmpString2,"tubsbox")==0){
 	  // Create a tube, from which center a box is cut out.  x1=box half-width; x2,x3,x4,x5 define the tube.
 	  sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %lf %s %lf %lf %lf %s %s",
@@ -347,7 +378,8 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	  //           x11, x12 are the half-withs of the ractangular opening in the collimator
 	  sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %s",
 		 name,&x1,&x2,&x3,&x4,&x5,&x6,&x7,&x8,&x9,&x10,&x11,&x12,material);
-	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*g %*g %*g %*g %*g %*g %*g %*s %lf %lf %lf %s %s %s %d %s",&posx,&posy,&posz,mothersName,rotMatrix,sensitiveDet,&volumeID,actualFieldName);
+	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*g %*g %*g %*g %*g %*g %*g %*s %lf %lf %lf %s %s %s %d %s",
+		 &posx,&posy,&posz,mothersName,rotMatrix,sensitiveDet,&volumeID,actualFieldName);
 	  solidName+=name;
 	  G4Box* solidDetBox = new G4Box("SolidDetBox",x1*mm,x2*mm,x3*mm);
 	  G4Box* solidHole   = new G4Box("SolidDetBox",x11*mm,x2*mm+0.1,x12*mm);
@@ -420,7 +452,6 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	  solid =  new G4SubtractionSolid(solidName,solidGPDmHolder,solidGPDcutOut,rot,trans);
 	  //	  solid =  new G4SubtractionSolid(solidName,solidGPDcutOut,solidGPDmHolder,rot,trans);
 	}
-
 	else ReportGeometryProblem(line);
 
 	G4ThreeVector position = G4ThreeVector (posx*mm,posy*mm,posz*mm);