StRoot: StiRnD/Hft/StiPixelDetectorBuilder.cxx Source File

00001 /*
00002  * $Id: StiPixelDetectorBuilder.cxx,v 1.27 2011/04/22 22:00:18 fisyak Exp $
00003  *
00004  * $Log: StiPixelDetectorBuilder.cxx,v $
00005  * Revision 1.27  2011/04/22 22:00:18  fisyak
00006  * warn off
00007  *
00008  * Revision 1.26  2010/08/25 21:57:41  fisyak
00009  * Get rid off access to specfic detector tracking parameters which usage has been  disable since 2008/06/11
00010  *
00011  * Revision 1.25  2009/03/16 13:51:00  fisyak
00012  * Move out all Sti Chairs into StDetectorDb
00013  *
00014  * Revision 1.24  2009/02/09 02:47:19  andrewar
00015  * UPGR15 update. Will break backward compatibility with older geometries.
00016  *
00017  * Revision 1.23  2008/04/03 20:04:20  fisyak
00018  * Straighten out DB access via chairs
00019  *
00020  * Revision 1.22  2007/10/20 00:16:27  fisyak
00021  * Active hit errors from DB
00022  *
00023  * Revision 1.21  2007/10/16 19:50:24  fisyak
00024  * rename Hft => Pxl, remove Hpd, Igt and Fst
00025  *
00026  * Revision 1.20  2007/05/16 15:02:57  andrewar
00027  * Removed couts in favor of LOG_INFO.
00028  *
00029  * Revision 1.19  2007/05/03 06:14:56  andrewar
00030  * Geometry fix to conform to StiHit:setGlobal() test.
00031  *
00032  * Revision 1.18  2007/03/30 02:14:19  andrewar
00033  * Removed some debug output.
00034  *
00035  * Revision 1.17  2006/11/30 16:37:19  andrewar
00036  * Removed call to dbase for tracking parameter loading for the review. Dynamic
00037  * access will be debugged and restored after the STAR review. Hit errors are
00038  * forced to 60um.
00039  *
00040  * Revision 1.16  2006/11/29 04:02:01  andrewar
00041  * Make use of pre-existing STAR DB inteface.
00042  *
00043  * Revision 1.15  2006/11/29 00:44:04  andrewar
00044  * Added call to get tracking parameters from DBase.
00045  *
00046  * Revision 1.14  2006/11/17 15:39:03  wleight
00047  * Changes to make PXL hits work with UPGR05 geometry
00048  *
00049  * Revision 1.13  2006/04/19 19:49:47  andrewar
00050  * Added call to setLayerAngle, needed for detector container sort.
00051  *
00052  * Revision 1.12  2006/02/23 00:22:54  andrewar
00053  * Set Detector Id to kPxlId, corrected Ist*pars -> Pixel*pars
00054  *
00055  * Revision 1.11  2006/02/17 21:39:32  andrewar
00056  * Added calls to StiDetector::setKey(key,val)
00057  *
00058  *
00059  */
00060 
00061 #include <stdio.h>
00062 #include <stdexcept>
00063 #include "Sti/StiPlanarShape.h"
00064 #include "Sti/StiCylindricalShape.h"
00065 #include "Sti/StiMaterial.h"
00066 #include "Sti/StiPlacement.h"
00067 #include "Sti/StiDetector.h"
00068 #include "Sti/Base/Factory.h"
00069 #include "Sti/StiToolkit.h"
00070 #include "StiPixelIsActiveFunctor.h"
00071 #include "Sti/StiNeverActiveFunctor.h"
00072 #include "Sti/StiElossCalculator.h"
00073 #include "StiPixelDetectorBuilder.h" 
00074 #include "StiPixelIsActiveFunctor.h"
00075 #include "StDetectorDbMaker/StiPixelHitErrorCalculator.h"
00076 #include "TDataSetIter.h"
00077 #include "tables/St_HitError_Table.h"
00078 #include "StEvent.h"
00079 #include "StEventTypes.h"
00080 
00081 StiPixelDetectorBuilder::StiPixelDetectorBuilder(bool active,
00082                                                  const string & inputFile)
00083   : StiDetectorBuilder("Pixel",active,inputFile)
00084 {
00085         //Parameterized hit error calculator.  Given a track (dip, cross, pt, etc)
00086         //returns average error once you actually want to do tracking, the results
00087         //depend strongly on the numbers below. 
00088 }
00089 
00090 StiPixelDetectorBuilder::~StiPixelDetectorBuilder()
00091 {}
00092 
00094 void StiPixelDetectorBuilder::buildDetectors(StMaker &source)
00095 {
00096 
00097   char name[50];
00098   LOG_INFO << "StiPixelDetectorBuilder::buildDetectors() -I- Started" << endm;
00099 
00100 
00101   unsigned int nRows=1;
00102 
00103   // 2 real rows, but we have detector elements and support elements. 
00104   setNRows(2);
00105   
00106   if (StiVMCToolKit::GetVMC()) {useVMCGeometry(); return;}
00107 
00108 
00109   //_gas is the gas that the pixel detector lives in
00110   _gasMat            = add(new StiMaterial("PixelAir",7.3, 14.61, 0.001205, 30420.*0.001205, 7.3*12.e-9));
00111   //_fcMaterial is the (average) material that makes up the detector elements.  Here I use ~silicon
00112   StiMaterial * material = add(new StiMaterial("PixelSi", 14.,  28.0855,   2.33,     21.82,   14.*12.*1e-9) );
00113 
00114 
00115   //Instantiate energy loss detector for si material  
00116   //const static double I2Ar = (15.8*18) * (15.8*18) * 1e-18; // GeV**2
00117   double ionization = material->getIonization();
00118   StiElossCalculator * siElossCalculator = new StiElossCalculator(material->getZOverA(),
00119                                                                   ionization*ionization,
00120                                                                   material->getA(),
00121                                                                   material->getZ(),
00122                                                                   material->getDensity());
00123   StiPlanarShape *pShape;
00124   for (unsigned int row=0; row<nRows; row++) 
00125     {
00126       pShape = new StiPlanarShape;
00127       if (!pShape) throw runtime_error("StiPixelDetectorBuilder::buildDetectors() - FATAL - pShape==0||ifcShape==0");
00128       sprintf(name, "Pixel/Layer_%d", row);
00129       pShape->setName(name);
00130       pShape->setThickness(0.0280); //cm 
00131       pShape->setHalfDepth( 20./2. );
00132       pShape->setHalfWidth(1.0);
00133       for(unsigned int sector = 0; sector<24; sector++) 
00134         {      
00135           StiPlacement *pPlacement = new StiPlacement;
00136           pPlacement->setZcenter(0.);
00137           double phi = phiForPixelSector(sector) + psiForPixelSector(sector);
00138           double r = radiusForPixelSector(sector)* cos(psiForPixelSector(sector)) - 0.0040; // note 40 microns offset
00139           double dY = radiusForPixelSector(sector)*sin(psiForPixelSector(sector));
00140           /*printf(" sector: %g phi: %g radius: %g normal: %g dY: %g\n",sector
00141                ,phi*180/3.1415
00142                << " radius:"<<radiusForPixelSector(sector)
00143                << " normal r:"<<r
00144                << "     dY:"<<dY<<endl;*/
00145           pPlacement->setNormalRep(phi, r, dY); 
00146           pPlacement->setLayerRadius(r);
00147           pPlacement->setLayerAngle(phi);
00148           pPlacement->setRegion(StiPlacement::kMidRapidity);
00149           sprintf(name, "Pixel/Layer_%d/Ladder_%d", row, sector);
00150           StiDetector *pDetector = _detectorFactory->getInstance();
00151           pDetector->setName(name);
00152           pDetector->setIsOn(true);
00153           pDetector->setIsActive(new StiPixelIsActiveFunctor);
00154           pDetector->setIsContinuousMedium(true);
00155           pDetector->setIsDiscreteScatterer(false);
00156           pDetector->setMaterial(material);
00157           pDetector->setGas(_gasMat);
00158           pDetector->setGroupId(kPxlId);
00159           pDetector->setShape(pShape);
00160           pDetector->setPlacement(pPlacement);
00161           pDetector->setHitErrorCalculator(StiPixelHitErrorCalculator::instance());
00162           pDetector->setElossCalculator(siElossCalculator);
00163           if (sector<18)
00164             {
00165               pDetector->setKey(1,1);
00166               pDetector->setKey(2,sector);
00167               add(1,sector,pDetector);
00168             }
00169           else
00170             {
00171               pDetector->setKey(1,0);
00172               pDetector->setKey(2,sector-18);
00173               add(0,(sector-18),pDetector);
00174             }
00175 
00176           //cout << "Setting detector: " << name << " with key values: "
00177           //     << pDetector->getKey(1) << " "  << pDetector->getKey(2) << endl;
00178         }
00179     }
00180   LOG_INFO << " -I- Done" << endl;
00181 }
00182 
00183 void StiPixelDetectorBuilder::useVMCGeometry() {
00184   LOG_INFO << "StiPixelDetectorBuilder::buildDetectors() -I- Use VMC geometry" 
00185        << endm;
00186   SetCurrentDetectorBuilder(this);
00187 
00188   // Build materials. In the Pixel detector we have two: Air for the mother volume,
00189   // silicon for both the detector and the ladder support. This will be updated with
00190   // more detailed support structures at the appropriate time.
00191   struct Material_t {
00192     const Char_t *name;
00193     StiMaterial    **p;
00194   };
00195   Material_t map[] = {
00196     {"AIR", &_gasMat},
00197     {"SILICON", &_fcMaterial} 
00198   };
00199   Int_t M = sizeof(map)/sizeof(Material_t);
00200   for (Int_t i = 0; i < M; i++) 
00201     {
00202       const TGeoMaterial *mat =  gGeoManager->GetMaterial(map[i].name); 
00203       if (! mat) continue;
00204       Double_t PotI = StiVMCToolKit::GetPotI(mat);
00205       *map[i].p = add(new StiMaterial(mat->GetName(),
00206                                     mat->GetZ(),
00207                                     mat->GetA(),
00208                                     mat->GetDensity(),
00209                                     mat->GetDensity()*mat->GetRadLen(),
00210                                     PotI));
00211     }
00212 
00213   //
00214   // Build volumes. Will be done from GEANT tables.
00215   //
00216 
00217   // Set volume name tree. Inactive volumes only here. Active volumes are declared in ::AverageVolume, called
00218   // through loop over StiDetectorBuilder::AverageVolume
00219   const VolumeMap_t PxlVolumes[] = 
00220     { 
00221       {"PLAC","Active ladder volume",   "HALL_1/CAVE_1/PXMO_1","",""},
00222       {"PLA1","Active ladder volume",   "HALL_1/CAVE_1/PXMO_1","",""},
00223       {"PLPS","Inactive ladder volume", "HALL_1/CAVE_1/PXMO_1","",""},
00224       {"PLP1","Inactive ladder volume", "HALL_1/CAVE_1/PXMO_1","",""}
00225     };
00226 
00227   Int_t NoPxlVols = sizeof(PxlVolumes)/sizeof(VolumeMap_t);
00228   TString pathT("HALL_1/CAVE_1/PXMO_1");
00229   gGeoManager->RestoreMasterVolume(); 
00230   gGeoManager->CdTop();
00231   TString path("");
00232   for (Int_t i = 0; i < NoPxlVols; i++) {
00233     gGeoManager->RestoreMasterVolume(); 
00234     gGeoManager->CdTop();
00235     if (gGeoManager->cd(pathT)) 
00236       {
00237         path = pathT;
00238         TGeoNode *nodeT = gGeoManager->GetCurrentNode();
00239         if (! nodeT) continue;;
00240         StiVMCToolKit::LoopOverNodes(nodeT, path, PxlVolumes[i].name, MakeAverageVolume);
00241       } 
00242     else {LOG_INFO << "StiPixelDetectorBuilder::useVMCGeometry skip node " 
00243                    << pathT.Data() << endm;}
00244   }
00245 }
00246 
00247 
00248 void StiPixelDetectorBuilder::AverageVolume(TGeoPhysicalNode *nodeP)
00249 {
00250   // Handle pathalogical input
00251   if (!nodeP) 
00252     {
00253       LOG_INFO << "StiPixelDetectorBuilder::AverageVolume -E- no TGeoPhysicalNode. "
00254                        << " Perhaps Pixel is turned on in tracking, but not present in simulation. Returning."
00255                        << endm;
00256      
00257       return;
00258     }
00259 
00260   // Note:
00261   // Volumes are currently all planes. I am coding this routine appropriately. Other
00262   // GEANT shapes (cylinder/sphere) are handled differently, and would require adding cases
00263   // if such are added to the PXL geometry.
00264   // AAR - Oct 31, 2006 
00265   
00266   TString nameP(nodeP->GetName());
00267   
00268   // Check whether this is an active volume
00269   Bool_t ActiveVolume = kFALSE;
00270   if (nodeP->GetVolume()->GetMedium()->GetParam(0) == 1) {
00271     ActiveVolume = kTRUE;
00272   }
00273   
00274   //Material definitions
00275   TGeoMaterial *matP   = nodeP->GetVolume()->GetMaterial(); 
00276   Double_t PotI = StiVMCToolKit::GetPotI(matP);
00277   static StiMaterial *matS = 0;
00278   if (! matS) matS = add(new StiMaterial(matP->GetName(),
00279                                          matP->GetZ(),
00280                                          matP->GetA(),
00281                                          matP->GetDensity(),
00282                                          matP->GetDensity()*matP->GetRadLen(),
00283                                          PotI));
00284   Double_t ionization = matS->getIonization();
00285   StiElossCalculator *ElossCalculator = new StiElossCalculator(matS->getZOverA(), 
00286                                                                ionization*ionization, 
00287                                                                matS->getA(), 
00288                                                                matS->getZ(),
00289                                                                matS->getDensity());
00290 
00291   //Extract volume geometry from TGeoPhysicalNode
00292   TGeoVolume   *volP   = nodeP->GetVolume();
00293   TGeoShape    *shapeP = nodeP->GetShape();  
00294   TGeoBBox     *box    = (TGeoBBox *) shapeP;
00295   StiShape     *sh     = new StiPlanarShape(volP->GetName(),        // Name
00296                                             box->GetDZ(),   // halfDepth
00297                                             box->GetDY(),           // thickness
00298                                             box->GetDX());          // halfWidth
00299 
00300 
00301   // position information
00302     TGeoHMatrix  *hmat   = nodeP->GetMatrix();  if (debug()) hmat->Print("");
00303     Double_t     *xyz    = hmat->GetTranslation();
00304     Double_t     *rot    = hmat->GetRotationMatrix();
00305     StThreeVectorD centerVector(xyz[0],xyz[1],xyz[2]);
00306     StThreeVectorD normalVector(rot[1],rot[4],rot[7]);
00307 
00308     //In ittf parlance, normalVector is the vector normal to the detector plane that passes through 0,0
00309     if( normalVector.magnitude() != 1. ) 
00310       {
00311         //printf("StiPixelDetectorBuilder -I- Normal vector is not a unit vector!\n");
00312         normalVector/=normalVector.magnitude();
00313       }
00314     normalVector *= centerVector.magnitude()*cos(normalVector.phi()-centerVector.phi());
00315     LOG_INFO <<"Setting detector normal angle: "
00316              << normalVector.phi()<<endm;
00317 
00318     double r = normalVector.perp();
00319     //Set the offset - this will be the difference between the two vectors. Set magnitude from vector difference...
00320     double dY = (normalVector - centerVector).magnitude();
00321     // set direction (+/-) from vector direction...
00322     // Remove auto-correction for change of ladder tilt. Assume dY is positive.
00323     // if ( (normalVector-centerVector).phi() < 0 ) dY=-1.*dY;
00324 
00325 
00326     StiPlacement *pPlacement = new StiPlacement;
00327     pPlacement->setZcenter(0);
00328     pPlacement->setLayerRadius(centerVector.perp()); //this is only used for ordering in detector container...
00329     pPlacement->setLayerAngle(centerVector.phi()); //this is only used for ordering in detector container...
00330     /*
00331     LOG_INFO << " -I- Setting detector center angle: " << centerVector.phi()
00332              << " offset: " << dY << " Ist style offset: "
00333              << centerVector.magnitude()*TMath::Sin(normalVector.phi() - centerVector.phi()) << endm;
00334     */
00335     pPlacement->setRegion(StiPlacement::kMidRapidity);
00336     pPlacement->setNormalRep(normalVector.phi(), r, dY); 
00337 
00338 
00339     // Add detector only if preceeding steps were successful
00340 
00341     if( sh ) add(sh);
00342     else 
00343     {
00344       LOG_INFO <<"StiPixelDetectorBuilder::AverageVolume() -E- StiPlanarShape build unsuccessful. Returning."
00345                <<endm;
00346       return;
00347     }
00348     if( pPlacement ) {}
00349     else
00350       {
00351         LOG_INFO <<"StiPixelDetectorBuilder::AverageVolume() -E- StiPlacement unsuccessful. "
00352                  << "Volume: "<<nameP.Data()<<". Returning."
00353                  <<endm;
00354         return;
00355       }
00356 
00357     //Build final detector object
00358     StiDetector *p =getDetectorFactory()->getInstance();
00359     if ( !p ) 
00360       {
00361         LOG_INFO <<"StiPixelDetectorBuilder::AverageVolume() -E- StiDetector pointer invalid." <<endm;
00362         return;
00363       }
00364     p->setName(nameP.Data());
00365 
00366     TString temp=nameP;
00367     temp.ReplaceAll("HALL_1/CAVE_1/PXMO_1/","");
00368     int q=temp.Index("_");
00369     temp.Replace(0,q+1,"");
00370     TString numlay=temp(0,1);
00371     int layer=numlay.Atoi();
00372     q=temp.Index("_");
00373     temp.Replace(0,q+1,"");
00374     TString numlad=temp(0,2);
00375     if(!numlad.IsDigit()) numlad=temp(0,1);
00376     int ladder=numlad.Atoi();
00377 
00378     p->setIsOn(false);
00379 
00380     p->setIsOn(false);
00381     if (ActiveVolume) {
00382       p->setIsActive(new StiPixelIsActiveFunctor);
00383     }
00384     else {
00385       p->setIsActive(new StiNeverActiveFunctor);
00386       layer=layer+2;
00387     }
00388 
00389     p->setIsContinuousMedium(false);
00390     p->setIsDiscreteScatterer(true);
00391     p->setShape(sh);
00392     p->setPlacement(pPlacement);
00393     p->setGas(GetCurrentDetectorBuilder()->getGasMat());
00394     if(!p->getGas()) LOG_INFO <<"gas not there!"<<endm;
00395     p->setMaterial(matS);
00396     p->setElossCalculator(ElossCalculator);
00397     p->setHitErrorCalculator(StiPixelHitErrorCalculator::instance());
00398     
00399     // Adding detector, not sure if setKey is necessary  
00400     //if(ActiveVolume){
00401     p->setKey(1, layer);
00402     p->setKey(2, ladder);
00403     add(layer,ladder, p);
00404     //}
00405     //else add(p);
00406     
00407     //LOG_INFO << "StiPixelDetectorBuilder: Added detector -I- " << p->getName() << endm;
00408     
00409     // Whole bunch of debugging information
00410     
00411     Float_t rad2deg = 180.0/3.1415927;
00412     LOG_DEBUG << "===>NEW:PIXEL:pDetector:Name               = " << p->getName()                               << endm;
00413     LOG_DEBUG << "===>NEW:PIXEL:pPlacement:NormalRefAngle    = " << pPlacement->getNormalRefAngle()*rad2deg    << endm;
00414     LOG_DEBUG << "===>NEW:PIXEL:pPlacement:NormalRadius      = " << pPlacement->getNormalRadius()              << endm;
00415     LOG_DEBUG << "===>NEW:PIXEL:pPlacement:NormalYoffset     = " << pPlacement->getNormalYoffset()             << endm;
00416     LOG_DEBUG << "===>NEW:PIXEL:pPlacement:CenterRefAngle    = " << pPlacement->getCenterRefAngle()*rad2deg    << endm;
00417     LOG_DEBUG << "===>NEW:PIXEL:pPlacement:CenterRadius      = " << pPlacement->getCenterRadius()              << endm;
00418     LOG_DEBUG << "===>NEW:PIXEL:pPlacement:CenterOrientation = " << pPlacement->getCenterOrientation()*rad2deg << endm;
00419     LOG_DEBUG << "===>NEW:PIXEL:pPlacement:LayerRadius       = " << pPlacement->getLayerRadius()               << endm;
00420     LOG_DEBUG << "===>NEW:PIXEL:pPlacement:LayerAngle        = " << pPlacement->getLayerAngle()*rad2deg        << endm;
00421     LOG_DEBUG << "===>NEW:PIXEL:pPlacement:Zcenter           = " << pPlacement->getZcenter()                   << endm;
00422     LOG_DEBUG << "===>NEW:PIXEL:pDetector:Layer              = " << layer                                      << endm;
00423     LOG_DEBUG << "===>NEW:PIXEL:pDetector:Ladder             = " << ladder                                     << endm;
00424     LOG_DEBUG << "===>NEW:PIXEL:pDetector:Active?            = " << p->isActive()                              << endm;
00425     
00426 
00427 }