StRoot: StiRnD/Ist/StiIstDetectorBuilder.cxx Source File

00001 // $Id: StiIstDetectorBuilder.cxx,v 1.22 2010/10/20 20:05:06 fisyak Exp $
00002 // 
00003 // $Log: StiIstDetectorBuilder.cxx,v $
00004 // Revision 1.22  2010/10/20 20:05:06  fisyak
00005 // Move SROD and SBSP from StiSvtDetectorBuilder to StiStarDetectorBuilder to account configurations without SVT detector installed (bug #2025)
00006 //
00007 // Revision 1.21  2010/08/25 21:57:41  fisyak
00008 // Get rid off access to specfic detector tracking parameters which usage has been  disable since 2008/06/11
00009 //
00010 // Revision 1.20  2009/03/16 13:51:00  fisyak
00011 // Move out all Sti Chairs into StDetectorDb
00012 //
00013 // Revision 1.19  2009/02/06 21:26:30  wleight
00014 // UPGR15 Update
00015 //
00016 // Revision 1.18  2008/04/03 20:04:21  fisyak
00017 // Straighten out DB access via chairs
00018 //
00019 // Revision 1.17  2007/04/27 18:44:03  wleight
00020 // Corrected a problem with incorrect assignment of hit errors
00021 //
00022 // Revision 1.16  2007/04/23 14:42:10  wleight
00023 // Added new hit error calculator for outer half of 17cm layer
00024 //
00025 // Revision 1.14  2007/04/06 15:58:21  wleight
00026 // Changed some cout statements to LOG_INFO
00027 //
00028 // Revision 1.13  2006/12/14 22:01:47  wleight
00029 // Changed hit errors so that they are obtained from the database and are different for each layer
00030 //
00031 // Revision 1.12  2006/11/28 22:18:01  wleight
00032 // Changed hit errors to 60 microns for x and 1.9 mm for y
00033 //
00034 // Revision 1.11  2006/10/20 18:43:12  wleight
00035 // Changes to make perfect hits in the IST work with UPGR05
00040 #include <stdio.h>
00041 #include <map>
00042 #include <exception>
00043 using namespace std;
00044 #include <stdexcept>
00045 #include "StMessMgr.h"
00046 #include "StThreeVectorD.hh"
00047 
00048 //#include "Sti/Base/Messenger.h"
00049 #include "Sti/Base/Factory.h"
00050 #include "Sti/StiPlanarShape.h"
00051 #include "Sti/StiCylindricalShape.h"
00052 #include "Sti/StiMaterial.h"
00053 #include "Sti/StiPlacement.h"
00054 #include "Sti/StiDetector.h"
00055 #include "Sti/StiToolkit.h"
00056 #include "Sti/StiIsActiveFunctor.h"
00057 #include "Sti/StiNeverActiveFunctor.h"
00058 #include "StiRnD/Ist/StiIstIsActiveFunctor.h" 
00059 #include "StiRnD/Ist/StiIstDetectorBuilder.h" 
00060 #include "Sti/StiElossCalculator.h"
00061 #include "StiRnD/Ist/StiIstDetectorBuilder.h"
00062 #include "StDetectorDbMaker/StiIst1HitErrorCalculator.h"
00063 //#include "StSsdUtil/StSsdConfig.hh"
00064 //#include "StSsdUtil/StSsdGeometry.hh"
00065 //#include "StSsdUtil/StSsdWaferGeometry.hh"
00066 //#include "StSsdDbMaker/StSsdDbMaker.h"
00067 //#include "StSsdDbMaker/St_SsdDb_Reader.hh"
00068 
00069 StiIstDetectorBuilder::StiIstDetectorBuilder(bool active, const string & inputFile)
00070     : StiDetectorBuilder("Ist",active,inputFile), _siMat(0), _hybridMat(0)
00071 {
00072     // Hit error parameters : it is set to 20 microns, in both x and y coordinates 
00073     //_hitCalculator1.setName("ist1HitError");
00074     //_hitCalculator2.setName("ist2HitError");
00075     //_hitCalculator3.setName("ist3HitError");
00076 }
00077 
00078 StiIstDetectorBuilder::~StiIstDetectorBuilder()
00079 {} 
00080 
00081 void StiIstDetectorBuilder::loadDS(TDataSet& ds){
00082   cout<<"StiIstDetectorBuilder::loadDS(TDataSet& ds) -I- started: "<<endl;
00083   //_hitCalculator1.loadDS(ds);
00084   //_hitCalculator2.loadDS(ds);
00085   //_hitCalculator3.loadDS(ds);
00086 }
00087 
00088 
00089 void StiIstDetectorBuilder::buildDetectors(StMaker & source)
00090 {
00091     int nRows = 1 ;
00092     gMessMgr->Info() << "StiIstDetectorBuilder::buildDetectors() - I - Started "<<endm;
00093     //load(_inputFile, source);
00094     
00095     setNRows(nRows);
00096     if (StiVMCToolKit::GetVMC()) {useVMCGeometry();}
00097 }
00098 
00099 //________________________________________________________________________________
00100 void StiIstDetectorBuilder::useVMCGeometry() {
00101   cout << "StiIstDetectorBuilder::buildDetectors() -I- Use VMC geometry" << endl;
00102   SetCurrentDetectorBuilder(this);
00103 
00104   // Build the material map
00105   struct Material_t {
00106     const Char_t *name;
00107     StiMaterial    **p;
00108   };
00109   Material_t map[] = {
00110     {"AIR", &_gasMat},
00111     {"SILICON", &_siMat},
00112     {"SILICON", &_hybridMat}
00113   };
00114   Int_t M = sizeof(map)/sizeof(Material_t);
00115   for (Int_t i = 0; i < M; i++) {
00116     const TGeoMaterial *mat =  gGeoManager->GetMaterial(map[i].name); 
00117     if (! mat) continue;
00118     Double_t PotI = StiVMCToolKit::GetPotI(mat);
00119     *map[i].p = add(new StiMaterial(mat->GetName(),
00120                                     mat->GetZ(),
00121                                     mat->GetA(),
00122                                     mat->GetDensity(),
00123                                     mat->GetDensity()*mat->GetRadLen(),
00124                                     PotI));
00125   }
00126 
00127   // Build the volume map and loop over all found volumes
00128   const VolumeMap_t IstVolumes[] = 
00129         { 
00130           {"IBSS", "active silicon",  "HALL_1/CAVE_1/IBMO_1","",""}, 
00131           {"IBSP", "inactive silicon","HALL_1/CAVE_1/IBMO_1","",""}
00132         };
00133   Int_t NoIstVols = sizeof(IstVolumes)/sizeof(VolumeMap_t);
00134   gGeoManager->RestoreMasterVolume(); 
00135   gGeoManager->CdTop();
00136   for (Int_t i = 0; i < NoIstVols; i++) {
00137     gGeoManager->cd(IstVolumes[i].path); 
00138     TGeoNode *nodeT = gGeoManager->GetCurrentNode();
00139     if (! nodeT) continue;;
00140     StiVMCToolKit::LoopOverNodes(nodeT, IstVolumes[i].path, IstVolumes[i].name, MakeAverageVolume);
00141   }
00142 }
00143 
00144 void StiIstDetectorBuilder::AverageVolume(TGeoPhysicalNode *nodeP) {
00145 
00146   LOG_DEBUG << "StiDetectorBuilder::AverageVolume -I TGeoPhysicalNode\t" << nodeP->GetName() << endm;
00147   //cout << "StiDetectorBuilder::AverageVolume -I TGeoPhysicalNode\t" << nodeP->GetName() << endl;
00148 
00149   // Ugh, have to hardwire
00150   const Int_t NWAFERS = 12;
00151 
00152   TString nameP(nodeP->GetName());
00153 
00154   // decode detector ------------------------------
00155   nameP.ReplaceAll("HALL_1/CAVE_1/","");
00156   TString temp=nameP;
00157   temp.ReplaceAll("/IBMO_1/IBMY","");
00158   int q=temp.Index("_");
00159   temp.Replace(0,q+1,"");
00160   TString num0=temp(0,1);
00161   int layer=num0.Atoi();
00162   q=temp.Index("_");
00163   temp.Replace(0,q+1,"");
00164   TString num1=temp(0,2);
00165   if(!num1.IsDigit()) num1=temp(0,1);
00166   int ladder=num1.Atoi(); 
00167   q=temp.Index("_");
00168   temp.Replace(0,q+1,"");
00169   TString num2=temp(0,2);
00170   if(!num2.IsDigit()) num2=temp(0,1);
00171   int wafer=num2.Atoi();
00172   q=temp.Index("_");
00173   temp.Replace(0,q+1,"");
00174   TString num3=temp(0,1);
00175   //  int side=num3.Atoi();
00176   if(wafer!=1) return;
00177 
00178   // Check whether this is an active volume
00179   Bool_t ActiveVolume = kFALSE;
00180   if (nodeP->GetVolume()->GetMedium()->GetParam(0) == 1) {
00181     ActiveVolume = kTRUE;
00182   }
00183 
00184   // Material definitions
00185   TGeoMaterial *matP = nodeP->GetVolume()->GetMaterial();
00186   Double_t PotI = StiVMCToolKit::GetPotI(matP);
00187   static StiMaterial *matS = 0;
00188   if (! matS) matS = add(new StiMaterial(matP->GetName(),
00189                                          matP->GetZ(),
00190                                          matP->GetA(),
00191                                          matP->GetDensity(),
00192                                          matP->GetDensity()*matP->GetRadLen(),
00193                                          PotI));
00194   Double_t ionization = matS->getIonization();
00195   StiElossCalculator *ElossCalculator = new StiElossCalculator(matS->getZOverA(), 
00196                                                                ionization*ionization, 
00197                                                                matS->getA(), 
00198                                                                matS->getZ(),
00199                                                                matS->getDensity());
00200   // Extract volume geometry for this node
00201   TGeoBBox *box = (TGeoBBox *) nodeP->GetShape();
00202   StiShape *sh  = new StiPlanarShape(nodeP->GetVolume()->GetName(), // Name
00203                                      NWAFERS*box->GetDZ(),          // halfDepth
00204                                      box->GetDY(),                  // thickness
00205                                      box->GetDX());                 // halfWidth
00206   add(sh);
00207     
00208   // Position information
00209   TGeoHMatrix  *hmat   = nodeP->GetMatrix();
00210   Double_t     *xyz    = hmat->GetTranslation();
00211   Double_t     *rot    = hmat->GetRotationMatrix();
00212   StThreeVectorD centerVector(xyz[0],xyz[1],xyz[2]);
00213   StThreeVectorD normalVector(rot[1],rot[4],rot[7]);
00214   //cout<<"centerVector: "<<centerVector<<endl;
00215   //cout<<"normalVector: "<<normalVector<<endl;
00216 
00217   // Normalize normal vector, just in case....
00218   normalVector /= normalVector.magnitude();
00219 
00220   // Volume positioning
00221   StiPlacement *pPlacement = new StiPlacement;
00222   Double_t phi  = centerVector.phi();
00223   Double_t phiD = normalVector.phi();
00224   Double_t r    = centerVector.perp();
00225   pPlacement->setZcenter(0);
00226   pPlacement->setLayerRadius(r); 
00227   /*
00228   double layerAngleOffset=(wafer-NWAFERS/2.)/25.;
00229   cout<<"offset: "<<layerAngleOffset<<endl;
00230   double layerAngle=phi+layerAngleOffset;
00231   cout<<"layerAngle: "<<layerAngle<<endl;
00232   pPlacement->setLayerAngle(layerAngle);
00233   */
00234   if(nameP.Contains("IBSS")) pPlacement->setLayerAngle(phi);
00235   if(nameP.Contains("IBSP")) pPlacement->setLayerAngle(phi-.05);
00236   pPlacement->setRegion(StiPlacement::kMidRapidity);
00237   pPlacement->setNormalRep(phiD, r*TMath::Cos(phi-phiD), r*TMath::Sin(phi-phiD)); 
00238   assert(pPlacement);
00239   //cout<<"normal radius: "<<pPlacement->getNormalRadius()<<endl;
00240   //cout<<"normal ref angle: "<<pPlacement->getNormalRefAngle()<<endl;
00241    
00242   StiDetector *pDetector = getDetectorFactory()->getInstance();
00243   pDetector->setName(nameP.Data());
00244   pDetector->setIsOn(false);
00245   if (ActiveVolume) {
00246     pDetector->setIsActive(new StiIstIsActiveFunctor);
00247     cout<<"active volume: "<<nameP<<endl;
00248   }
00249   else {
00250     pDetector->setIsActive(new StiNeverActiveFunctor);
00251     cout<<"inactive volume: "<<nameP<<endl;
00252     layer++;
00253   }
00254   pDetector->setIsContinuousMedium(false);
00255   pDetector->setIsDiscreteScatterer(true);
00256   pDetector->setShape(sh);
00257   pDetector->setPlacement(pPlacement); 
00258   pDetector->setGas(GetCurrentDetectorBuilder()->getGasMat());
00259   pDetector->setMaterial(matS);
00260   pDetector->setElossCalculator(ElossCalculator);
00261   pDetector->setHitErrorCalculator(StiIst1HitErrorCalculator::instance());
00262 
00263   // Adding detector, note that no keys are set in IST!
00264   //add(ladder,wafer,pDetector);
00265   add(layer,ladder,pDetector);
00266 
00267   LOG_INFO << "StiIstDetectorBuilder: Added detector -I- " << pDetector->getName() << endm;
00268 
00269   // Whole bunch of debugging information
00270   /*
00271   Float_t rad2deg = 180.0/3.1415927;
00272   cout << "===>NEW:IST:pDetector:Name               = " << pDetector->getName()                       << endl;
00273   cout << "===>NEW:IST:pPlacement:NormalRefAngle    = " << pPlacement->getNormalRefAngle()*rad2deg    << endl;
00274   cout << "===>NEW:IST:pPlacement:NormalRadius      = " << pPlacement->getNormalRadius()              << endl;
00275   cout << "===>NEW:IST:pPlacement:NormalYoffset     = " << pPlacement->getNormalYoffset()             << endl;
00276   cout << "===>NEW:IST:pPlacement:CenterRefAngle    = " << pPlacement->getCenterRefAngle()*rad2deg    << endl;
00277   cout << "===>NEW:IST:pPlacement:CenterRadius      = " << pPlacement->getCenterRadius()              << endl;
00278   cout << "===>NEW:IST:pPlacement:CenterOrientation = " << pPlacement->getCenterOrientation()*rad2deg << endl;
00279   cout << "===>NEW:IST:pPlacement:LayerRadius       = " << pPlacement->getLayerRadius()               << endl;
00280   cout << "===>NEW:IST:pPlacement:LayerAngle        = " << pPlacement->getLayerAngle()*rad2deg        << endl;
00281   cout << "===>NEW:IST:pPlacement:Zcenter           = " << pPlacement->getZcenter()                   << endl;
00282   cout << "===>NEW:IST:pDetector:Layer             = " << layer                                       << endl;
00283   cout << "===>NEW:IST:pDetector:Ladder            = " << ladder                                      << endl;
00284   cout << "===>NEW:IST:pDetector:Wafer             = " << wafer                                       << endl;
00285   cout << "===>NEW:IST:pDetector:Side(defunct)     = " << side                                        << endl;
00286   cout << "===>NEW:IST:pDetector:Active?           = " << pDetector->isActive()                       << endl;
00287   */
00288 
00289 }
00290