StRoot: Sti/StiDetectorBuilder.cxx Source File

00001 #include "Stiostream.h"
00002 #include "Sti/Base/Factory.h"
00003 #include "Sti/StiDetector.h"
00004 #include "Sti/StiPlanarShape.h"
00005 #include "StiCylindricalShape.h"
00006 #include "Sti/StiPlacement.h"
00007 #include "Sti/StiMaterial.h"
00008 #include "Sti/StiDetectorBuilder.h"
00009 #include "Sti/StiToolkit.h"
00010 #include "Sti/StiNeverActiveFunctor.h"
00011 #include "Sti/StiElossCalculator.h"
00012 #include "StDetectorDbMaker/StiDefaultTrackingParameters.h"
00013 #include "StThreeVector.hh"
00014 #include "StMaker.h"
00015 #include "StThreeVectorD.hh"
00016 #include "TMath.h"
00017 StiDetectorBuilder* StiDetectorBuilder::fCurrentDetectorBuilder = 0;
00018 int StiDetectorBuilder::_debug = 0;
00019 StiDetectorBuilder::StiDetectorBuilder(const string & name,bool active, const string & inputFile)
00020   : Named(name+"Builder"),
00021     _groupId(-1),
00022     _active(active),
00023     _detectorFactory( StiToolkit::instance()->getDetectorFactory() ),
00024     _trackingParameters(0),
00025     _inputFile(inputFile),
00026     _gasMat(0)
00027 {
00028   cout << "StiDetectorBuilder::StiDetectorBuilder() - INFO - Instantiating builder named:"<<name<<endl;
00029   fCurrentDetectorBuilder = this;
00030 }
00031 
00032 StiDetectorBuilder::~StiDetectorBuilder()
00033 {}
00034 
00035 bool StiDetectorBuilder::hasMore() const 
00036 {
00037   //cout<<"StiDetectorBuilder::hasMore() - INFO - Started"<<endl;
00038   return mDetectorIterator != mDetectorMap.end();
00039 } // hasMore()
00040 
00041 StiDetector * StiDetectorBuilder::next()
00042 {
00043   //cout<<"StiDetectorBuilder::hasMore() - INFO - Started"<<endl;
00044   if (mDetectorIterator != mDetectorMap.end())
00045     return (mDetectorIterator++)->second;
00046   else 
00047     return 0;
00048 } // next()
00049 
00050 StiMaterial* StiDetectorBuilder::findMaterial(const string& szName) const
00051 {
00052   materialMap::const_iterator where = mMaterialMap.find(NameMapKey(szName));
00053   return (where!= mMaterialMap.end()) ? (*where).second : 0;
00054 } // findMaterial()
00055 
00056 StiShape* StiDetectorBuilder::findShape(const string& szName) const
00057 {
00058   shapeMap::const_iterator where = mShapeMap.find(NameMapKey(szName));
00059   return (where!=mShapeMap.end()) ? (*where).second: 0;
00060 } // findShape()
00061 
00062 StiDetector* StiDetectorBuilder::findDetector(const string& szName) const
00063 {
00064   detectorMap::const_iterator where = mDetectorMap.find(NameMapKey(szName));
00065   return (where!=mDetectorMap.end()) ? (*where).second: 0;
00066 } // findDetector()
00067 
00068 StiMaterial * StiDetectorBuilder::add(StiMaterial *material)
00069 {  
00070   NameMapKey key(material->getName());
00071   mMaterialMap.insert( materialMapValType(key,material) );
00072   return material;
00073 }
00074 
00075 StiShape * StiDetectorBuilder::add(StiShape *shape)
00076 {
00077   NameMapKey key(shape->getName());
00078   mShapeMap.insert( shapeMapValType(key, shape) );
00079         return shape;
00080 }
00081 
00082 StiDetector * StiDetectorBuilder::add(unsigned int row, unsigned int sector, StiDetector *detector)
00083 {
00084   setNSectors(row,sector+1);
00085   _detectors[row][sector] = detector;
00086   if (_debug || sector == 0) {
00087     cout << "StiDetectorBuilder::add(" << row << "," << sector << ") detector ";
00088     if (detector) cout << detector->getName();
00089     else          cout << " NULL ??";
00090     cout <<endl;
00091   }
00092   return add(detector);
00093 }
00094 
00098 StiDetector * StiDetectorBuilder::add(StiDetector *detector)
00099 {
00100   NameMapKey key(detector->getName());
00101   mDetectorMap.insert( detectorMapValType(key, detector) );
00102   //complete the building of this detector element
00103   // in the base class nothing is actually done
00104   // but ROOT stuff is built in the drawable version of this class.
00105   detector->setGroupId(_groupId);
00106   detector->setTrackingParameters(StiDefaultTrackingParameters::instance());
00107   return detector;
00108 }
00109 
00110 void StiDetectorBuilder::build(StMaker& source)
00111 {
00112   buildDetectors(source);
00113   mDetectorIterator = mDetectorMap.begin();
00114 }
00115 
00116 //________________________________________________________________________________
00117 void StiDetectorBuilder::buildDetectors(StMaker& source)
00118 {}
00119 //________________________________________________________________________________
00120 void StiDetectorBuilder::AverageVolume(TGeoPhysicalNode *nodeP) {
00121   if (debug()) {cout << "StiDetectorBuilder::AverageVolume -I TGeoPhysicalNode\t" << nodeP->GetName() << endl;}
00122   TGeoVolume   *volP   = nodeP->GetVolume();
00123   TGeoMaterial *matP   = volP->GetMaterial(); if (debug()) matP->Print("");
00124   TGeoShape    *shapeP = nodeP->GetShape();   if (debug()) {cout << "New Shape\t"; StiVMCToolKit::PrintShape(shapeP);}
00125   TGeoHMatrix  *hmat   = nodeP->GetMatrix();  if (debug()) hmat->Print("");
00126   Double_t PotI = StiVMCToolKit::GetPotI(matP);
00127   StiMaterial *matS = add(new StiMaterial(matP->GetName(),
00128                                                 matP->GetZ(),
00129                                                 matP->GetA(),
00130                                                 matP->GetDensity(),
00131                                                 matP->GetDensity()*matP->GetRadLen(),
00132                                                 PotI));
00133   Double_t ionization = matS->getIonization();
00134   StiElossCalculator *ElossCalculator = 
00135     new StiElossCalculator(matS->getZOverA(), ionization, matS->getA(), matS->getZ(),matS->getDensity());
00136   StiShape     *sh     = findShape(volP->GetName());
00137   Double_t     *xyz    = hmat->GetTranslation();
00138   Double_t     *rot    = hmat->GetRotationMatrix();
00139   Double_t      Phi    = 0;
00140   //  Double_t xc,yc,zc,rc,rn, nx,ny,nz,yOff;
00141   StiPlacement *pPlacement = 0;
00142   do {//only once
00143     if (!shapeP->TestShapeBit(TGeoShape::kGeoTube))     break;
00144     TGeoTube *shapeC = (TGeoTube *)shapeP;
00145     Double_t Rmax = shapeC->GetRmax();
00146     Double_t Rmin = shapeC->GetRmin();
00147     Double_t delta=fabs(xyz[0])+fabs(xyz[1]);
00148     if (delta>0.1*Rmin)                                 break;
00149     Double_t dZ   = shapeC->GetDz();
00150     Double_t radius = (Rmin + Rmax)/2;
00151     Double_t dPhi = 2*TMath::Pi();
00152     Double_t dR   = Rmax - Rmin;
00153     dR = TMath::Min(0.2*dZ, dR);
00154     if (dR < 0.1) dR = 0.1;
00155     int Nr = (int) ((Rmax - Rmin)/dR);
00156     if (Nr <= 0) Nr = 1;
00157     dR = (Rmax - Rmin)/Nr;
00158 
00159     if (shapeP->TestShapeBit(TGeoShape::kGeoTubeSeg)) {
00160       TGeoTubeSeg *shapeS = (TGeoTubeSeg *) shapeP;
00161       Double_t gloV[3];
00162       Double_t Phi1 = TMath::DegToRad()*shapeS->GetPhi1();
00163       Double_t Phi2 = TMath::DegToRad()*shapeS->GetPhi2();
00164       if (Phi2<Phi1) Phi2+=M_PI*2;
00165       Double_t PhiM = (Phi2+Phi1)/2;
00166       dPhi        = (Phi2-Phi1);
00167       Double_t locV[3]={cos(PhiM),sin(PhiM),0};
00168       hmat->LocalToMaster(locV,gloV);                  
00169       Phi = atan2(gloV[1],gloV[0]);
00170     }
00171     
00172     for (Int_t ir = 0; ir < Nr; ir++) {
00173       TString Name(volP->GetName());
00174       if (ir) {Name += "__";Name += ir;}
00175       sh     = findShape(Name.Data());
00176       if (! sh) {// I assume that the shape name is unique
00177         sh = new StiCylindricalShape(Name.Data(),   // Name
00178                                      dZ,              // halfDepth
00179                                      dR,              // thickness
00180                                      Rmin + (ir+1)*dR,// outerRadius
00181                                      dPhi);           // openingAngle
00182         add(sh);
00183       }
00184       pPlacement = new StiPlacement;
00185       pPlacement->setZcenter(xyz[2]);
00186       pPlacement->setLayerRadius(Rmin + (ir+0.5)*dR);
00187       pPlacement->setLayerAngle(Phi);
00188       pPlacement->setRegion(StiPlacement::kMidRapidity);
00189       pPlacement->setNormalRep(Phi,radius, 0); 
00190     }
00191   } while(0);
00192 
00193   if (!pPlacement)  {// BBox
00194     shapeP->ComputeBBox();
00195     TGeoBBox *box = (TGeoBBox *)shapeP;
00196     Double_t dx = box->GetDX();
00197     Double_t dy = box->GetDY();
00198     Double_t dz = box->GetDZ();
00199     StThreeVectorD centerVector(xyz[0],xyz[1],xyz[2]);
00200     Double_t r = centerVector.perp();
00201     Double_t phi  = centerVector.phi();
00202     Int_t ix = -1;
00203     Double_t distMax = 0;
00204     Double_t dist;
00205     Bool_t swap = kFALSE;
00206     for (Int_t i = 0; i < 3; i++) {
00207       StThreeVectorD normalVector(rot[i],rot[i+3],rot[i+6]);
00208       Double_t phiD = normalVector.phi();
00209       dist = r*TMath::Cos(phi-phiD);
00210       swap = kFALSE;
00211       if (dist < -1e-7) {normalVector *= -1; swap = kTRUE;}
00212       phiD = normalVector.phi();
00213       dist = r*TMath::Cos(phi-phiD);
00214       if (! swap) dist += 1;
00215       if (dist > distMax) {distMax = dist; ix = i; }
00216     }
00217     assert( ix != -1);
00218     StThreeVectorD normalVector(rot[ix],rot[ix+3],rot[ix+6]);
00219     Double_t prod = centerVector*normalVector;
00220     if (prod < -1e-7) normalVector *= -1;
00221     Double_t phiD = normalVector.phi();
00222     if (ix == 0) {
00223       dx = dy;
00224       dy = box->GetDX();
00225     }
00226     if (ix == 2) {
00227       dy = dz;
00228       dz = box->GetDY();
00229     }
00230     if (! sh) {
00231       sh = new StiPlanarShape(volP->GetName(),// Name
00232                               dz,             // halfDepth
00233                             2*dy,             // thickness
00234                               dx);            // halfWidth
00235       add(sh);
00236     }
00237     pPlacement = new StiPlacement;
00238     pPlacement->setZcenter(xyz[2]);
00239     pPlacement->setLayerRadius(r); //this is only used for ordering in detector container...
00240     pPlacement->setLayerAngle(phi); //this is only used for ordering in detector container...
00241     pPlacement->setRegion(StiPlacement::kMidRapidity);
00242     pPlacement->setNormalRep(phiD, r*TMath::Cos(phi-phiD), r*TMath::Sin(phi-phiD)); 
00243   }
00244   assert(pPlacement);
00245   StiDetector *pDetector = getDetectorFactory()->getInstance();
00246   TString nameP(nodeP->GetName());
00247   nameP.ReplaceAll("HALL_1/CAVE_1/","");
00248   nameP.Resize(30); nameP.Strip();
00249   pDetector->setName(nameP.Data());
00250   pDetector->setIsOn(false);
00251   pDetector->setIsActive(new StiNeverActiveFunctor);
00252   pDetector->setIsContinuousMedium(false);
00253   pDetector->setIsDiscreteScatterer(true);
00254   pDetector->setShape(sh);
00255   pDetector->setPlacement(pPlacement); 
00256   pDetector->setGas(GetCurrentDetectorBuilder()->getGasMat());
00257   pDetector->setMaterial(matS);
00258   pDetector->setElossCalculator(ElossCalculator);
00259   Int_t layer = getNRows();
00260   add(layer+1,0,pDetector); 
00261   cout << "StiDetectorBuilder::AverageVolume build detector " << pDetector->getName() << " at layer " << layer << endl;
00262 }
00266 unsigned int  StiDetectorBuilder::getNSectors(unsigned int row) const
00267 {
00268   assert(row<_detectors.size());
00269   return _detectors[row].size();
00270 }
00271 
00272 
00273 StiDetector * StiDetectorBuilder::getDetector(unsigned int row, unsigned int sector) const
00274 {
00275   assert(row<_detectors.size());
00276   assert(sector<_detectors[row].size());
00277   return _detectors[row][sector];
00278 }
00279 
00280 void StiDetectorBuilder::setDetector(unsigned int row, unsigned int sector, StiDetector *detector)
00281 {
00282   setNSectors(row+1,sector+1);
00283    _detectors[row][sector] = detector;
00284 }