31.8.2011 Kamil Sedlak

New geometry called "polycone" implemented in musrSim.
This required also an implementation of command
"/musr/command arrayDef arrayName N x1 x2 ... xN"
See the updated documentation musrSim.pdf.

doc/musrSim.tex

@@ -277,10 +277,10 @@ care.
         logical volume and physical volume.  Details can be found in {\sc Geant4} manual). \\
 	\begin{itemize}
 	\item \emph{solid} (string) can be one of the G4VSolid.cc particular types, presently ``tubs'', ``cons'',
-	             ``box'', ``trd'', ``sphere'', ``para'',
+	             ``box'', ``trd'', ``sphere'', ``para'', ``polycone''
 	             or it can be one of the specifically implemented solids by our program as ``uprofile'' 
 	             (an U-profiled bar), ``alcSupportPlate'' (shape specific to ALC support plate), ``tubsbox'' 
-	             (a tube with a rectangular hole along its axis), "tubsboxsegm" 
+	             (a tube with a rectangular hole along its axis), ``tubsboxsegm'' 
 	             (a volume that looks like an intersection of tube and box) and 
                      ``trd90y'' (a trd volume rotated by 90 degrees around $y$ axis in addition
                      to the rotation requested by \emph{matrixName}), 
@@ -289,6 +289,15 @@ care.
 		     in Geant4 are
                      presently supported, but it is relatively easy to implement a new kind of solids
                      in the musrDetectorConstruction.cc class.
+		     
+		     The ``polycone'' in Geant4 can be defined in two ways, which have to be 
+		     distinguished in the musrSim macro file:
+		     ``polyconeA'' corresponds to \\
+		     G4Polycone(solidName, phiStart, phiTotal, numZPlanes, zPlane[ ], rInner[ ], rOuter[ ])\\
+		     while  ``polyconeB'' corresponds to \\
+		     G4Polycone(solidName, phiStart, phiTotal, numRZ, r[ ], z[ ]),\\
+		     where \emph{zPlane}, \emph{rInner}, \emph{rOuter}, \emph{r}, \emph{z} are arrays, which
+		     have to be defined by the command ``/musr/command arrayDef''.
 	\item \emph{name} (string) stands for the name of the volume.  As the command 
 		     ``/musr/command construct'' constructs
 	             three kinds of Geant4 classes/volumes (the solid, logical volume and physical 
@@ -344,6 +353,9 @@ care.
 	vector (if the fourth float parameter behind the \emph{matrixName} is non-zero).
 	All angles are specified in degrees.
 
+\item{\bf /musr/command arrayDef \emph{arrayName} \emph{N} \emph{x$_1$} \emph{x$_2$} \ldots \emph{x$_N$}} \\
+        Defines an array with an assigned name \emph{arrayName}, which has \emph{N} elements.
+
 \item{\bf /musr/command region define  \emph{regionName} \emph{logicalVolume}}\\
 	The ``G4Region'' can be created using this command, and a logical volume of the
 	name \emph{logicalVolume} will be assigned to it.  If the G4Region of the name 

include/musrDetectorConstruction.hh

@@ -74,6 +74,8 @@ private:
   musrDetectorMessenger* detectorMessenger;  // pointer to the Messenger
 
   std::map<std::string,G4RotationMatrix*> pointerToRotationMatrix;
+  std::map<std::string,double*>           pointerToArray;
+  std::map<std::string,double*>::iterator iterArray;
   std::map<std::string,G4FieldManager*> pointerToField;
 
   std::map<std::string,G4MaterialPropertiesTable*> materialPropertiesTableMap;

musrSimAna/musrAnalysis.cxx

@@ -35,6 +35,8 @@ void musrAnalysis::Loop(char* runChar, char* v1190FileName, Int_t nrEvents)
 //by  b_branchname->GetEntry(ientry); //read only this branch
    if (fChain == 0) return;
 
+   std::cerr<<"musrSimAnalysis::Loop() :  Analysing run "<<v1190FileName<<std::endl;
+
    hInfo = new TH1D("hInfo","Different values that need to be passed to ploting program",10000,1.,10001.);
    ReadInInputParameters(v1190FileName);
    CreateHistograms();
@@ -707,7 +709,9 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
 	int iChn2 = itCoinc->second;
 	counterMapType::const_iterator itTwo = allCounterMap.find(abs(iChn2));
         if (itTwo==allCounterMap.end()) { 
-	  std::cout<<" Pointer to coincidence counter ("<<iChn2<<") not found!  This should never happen!!! ==> S T O P"<<std::endl;
+	  std::cout<<" Pointer to coincidence counter ("<<iChn2<<") not found!    Perhaps the counter nr. "<<abs(iChn2)
+		   <<" was not defined ?"<<std::endl;
+	  std::cout<<" Serious problem !!! ==> S T O P"<<std::endl;
 	  exit(1);
 	}
 	else {

src/musrDetectorConstruction.cc

@@ -43,6 +43,7 @@
 #include "G4NistManager.hh"
 #include "G4Box.hh"
 #include "G4Cons.hh"
+#include "G4Polycone.hh"
 #include "G4LogicalVolume.hh"
 #include "G4PVPlacement.hh"
 #include "G4Tubs.hh"
@@ -165,7 +166,43 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	}
       }
 
-      
+      // Define array, which is needed for G4Polycone definition
+      else if (strcmp(tmpString1,"arrayDef")==0){
+	const int maxArrayLength=100;
+	double* newArray = new double[maxArrayLength];
+	char arrayName[100]="Unset";
+	int  nValues;
+        sscanf(&line[0],"%*s %*s %s %d",arrayName,&nValues);
+	if (nValues>maxArrayLength) {
+	  G4cout<<"musrDetectorConstruction: arrayDef: currently the array can not have more than "<<maxArrayLength<<" elements."<<G4endl;
+	  G4cout<<"      Please either increse the \"maxArrayLength\" variable in musrDetectorConstruction.cc or use shorter arrays."<<G4endl;
+	  G4cout << "S T O P    F O R C E D"<<G4endl;
+	  G4cout<<"   "<<line<<G4endl;
+	  exit(1);
+	}
+	else {
+	  char NAME1[100];
+	  double pp1;
+	  char* pch = strstr(line,arrayName)+strlen(arrayName);
+	  for (int iji=0; iji<nValues; iji++) {
+	    sscanf(pch,"%s",NAME1);
+	    char *pch1 = strstr(pch,NAME1)+strlen(NAME1);
+	    pch=pch1;
+	    sscanf(pch,"%lf",&pp1);
+	    newArray[iji]=pp1;
+	  }
+	}
+	G4cout<<"New array defined: "<<arrayName<<"=(";
+	for (int iji=0; iji<nValues; iji++) {
+	  G4cout<<newArray[iji]; 
+	  if (iji<nValues-1) {G4cout<<",";}
+	  else             {G4cout<<")";}
+	}
+	G4cout<<G4endl;
+	pointerToArray[arrayName]=newArray;
+      }
+	
+     
       else if (strcmp(tmpString1,"construct")==0){
 	double x1=0,x2=0,x3=0,x4=0,x5=0,x6=0,x7=0,x8=0,x9=0,x10=0,x11,x12; 
 	double posx,posy,posz;
@@ -217,6 +254,74 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	  solidName+=name;
 	  solid = new G4Sphere(solidName,x1*mm,x2*mm,x3*deg,x4*deg,x5*deg,x6*deg);
 	}
+	else if ((strcmp(tmpString2,"polyconeA")==0)||(strcmp(tmpString2,"polyconeB")==0)){
+	  G4bool polyconeA=true; if (strcmp(tmpString2,"polyconeB")==0) polyconeA=false;
+	  char zPlane[100]; double* zPLANE;
+	  char rInner[100]; double* rINNER;
+	  char rOuter[100]; double* rOUTER;
+	  int numZPlanes;
+	  if (polyconeA) {
+	    sscanf(&line[0],"%*s %*s %*s %s %lf %lf %d %s %s %s %s %lf %lf %lf %s %s",
+		   name,&x1,&x2,&numZPlanes,zPlane,rInner,rOuter,material,&posx,&posy,&posz,mothersName,rotMatrix);
+	    sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*d %*s %*s %*s %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
+	  }
+	  else {  // polyconeB
+	    sscanf(&line[0],"%*s %*s %*s %s %lf %lf %d %s %s %s %lf %lf %lf %s %s",
+		   name,&x1,&x2,&numZPlanes,zPlane,rInner,material,&posx,&posy,&posz,mothersName,rotMatrix);
+	    sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*d %*s %*s %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
+	  }
+	  solidName+=name;
+	  G4cout<<name<<G4endl;
+	  //	  solid = new G4Sphere(solidName,x1*mm,x2*mm,x3*deg,x4*deg,x5*deg,x6*deg);
+	  //	  double* zPLANE = pointerToArray[zPlane];
+	  iterArray = pointerToArray.find(zPlane);
+	  if (iterArray==pointerToArray.end()) {  // array does not exist
+	    G4cout<<"\n\n musrDetectorConstruction(): array \""<<zPlane<<"\" requested for the volume \""
+		  <<name<<"\",\n  but the array was not found -> was it defined?"<<G4endl;
+	    G4cout << "S T O P    F O R C E D"<<G4endl;
+	    G4cout << line << G4endl;
+	    exit(1);
+	  }
+	  else {
+	    zPLANE = iterArray->second;
+	  }
+	  G4cout<<"zPLANE="; for (int iiiii=0; iiiii<numZPlanes; iiiii++) {G4cout<<zPLANE[iiiii]<<",";} G4cout<<G4endl;
+	  //
+	  iterArray = pointerToArray.find(rInner);
+	  if (iterArray==pointerToArray.end()) {  // array does not exist
+	    G4cout<<"\n\n musrDetectorConstruction(): array \""<<rInner<<"\" requested for the volume \""
+		  <<name<<"\",\n  but the array was not found -> was it defined?"<<G4endl;
+	    G4cout << "S T O P    F O R C E D"<<G4endl;
+	    G4cout << line << G4endl;
+	    exit(1);
+	  }
+	  else {
+	    rINNER = iterArray->second;
+	  }
+	  G4cout<<"rINNER="; for (int iiiii=0; iiiii<numZPlanes; iiiii++) {G4cout<<rINNER[iiiii]<<",";} G4cout<<G4endl;
+	  //
+	  if (polyconeA) {
+	    iterArray = pointerToArray.find(rOuter);
+	    if (iterArray==pointerToArray.end()) {  // array does not exist
+	      G4cout<<"\n\n musrDetectorConstruction(): array \""<<rOuter<<"\" requested for the volume \""
+		    <<name<<"\",\n  but the array was not found -> was it defined?"<<G4endl;
+	      G4cout << "S T O P    F O R C E D"<<G4endl;
+	      G4cout << line << G4endl;
+	      exit(1);
+	    }
+	    else {
+	      rOUTER = iterArray->second;
+	    }
+	    G4cout<<"rOUTER="; for (int iiiii=0; iiiii<numZPlanes; iiiii++) {G4cout<<rOUTER[iiiii]<<",";} G4cout<<G4endl; 
+	  }
+	  //
+	  if (polyconeA) {
+	    solid = new G4Polycone(solidName,x1*deg,x2*deg,numZPlanes,zPLANE,rINNER,rOUTER);
+	  }
+	  else {  //  polyconeB
+	   solid = new G4Polycone(solidName,x1*deg,x2*deg,numZPlanes,zPLANE,rINNER);
+	  } 
+	}
 	else if (strcmp(tmpString2,"para")==0){  // NOT YET TESTED
 	  sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %lf %lf %s %lf %lf %lf %s %s",
 		 name,&x1,&x2,&x3,&x4,&x5,&x6,material,&posx,&posy,&posz,mothersName,rotMatrix);