StRoot: StiTpc/StiTpcDetectorBuilder.cxx Source File

00001 #include <assert.h>
00002 #include <stdio.h>
00003 #include "Stiostream.h"
00004 #include <stdexcept>
00005 #include "StDbUtilities/StTpcLocalCoordinate.hh"
00006 #include "StDbUtilities/StTpcCoordinateTransform.hh"
00007 #include "StTpcDb/StTpcDb.h"
00008 #include "Sti/StiPlanarShape.h"
00009 #include "Sti/StiCylindricalShape.h"
00010 #include "Sti/StiMaterial.h"
00011 #include "Sti/StiPlacement.h"
00012 #include "Sti/StiDetector.h"
00013 #include "Sti/Base/Factory.h"
00014 #include "Sti/StiToolkit.h"
00015 #include "Sti/StiIsActiveFunctor.h"
00016 #include "Rtypes.h"
00017 #include "Stiostream.h"
00018 #include "Sti/StiNeverActiveFunctor.h"
00019 #include "StDetectorDbMaker/StiTpcInnerHitErrorCalculator.h"
00020 #include "StDetectorDbMaker/StiTpcOuterHitErrorCalculator.h"
00021 #include "StiTpcDetectorBuilder.h"
00022 #include "StiTpcIsActiveFunctor.h"
00023 #include "Sti/StiElossCalculator.h"
00024 #include "StDetectorDbMaker/StDetectorDbTpcRDOMasks.h"
00025 #include "StDetectorDbMaker/St_tpcPadPlanesC.h"
00026 #include "StDbUtilities/StCoordinates.hh"
00027 #include "StTpcDb/StTpcDb.h"
00028 #include "StMatrixD.hh"
00029 #include "StDetectorDbMaker/St_tpcAnodeHVavgC.h"
00030 #include "StDetectorDbMaker/St_tpcPadGainT0C.h"
00031 //#define TPC_IDEAL_GEOM
00032 
00033 StiTpcDetectorBuilder::StiTpcDetectorBuilder(Bool_t active, const string & inputFile)
00034   : StiDetectorBuilder("Tpc",active,inputFile), _fcMaterial(0){}
00035 
00036 StiTpcDetectorBuilder::~StiTpcDetectorBuilder() {}
00037 
00046 void StiTpcDetectorBuilder::buildDetectors(StMaker&source)
00047 {
00048   cout << "StiTpcDetectorBuilder::buildDetectors() -I- Started" << endl;
00049   if (!gStTpcDb)
00050     throw runtime_error("StiTpcDetectorBuilder::buildDetectors() -E- gStTpcDb==0");
00051   useVMCGeometry();
00052   cout << "StiTpcDetectorBuilder::buildDetectors() -I- Done" << endl;
00053 }
00054 //________________________________________________________________________________
00055 void StiTpcDetectorBuilder::useVMCGeometry() {
00056   Int_t debug = 0;
00057 
00058   if (debug>1) StiVMCToolKit::SetDebug(1);
00059   cout << "StiTpcDetectorBuilder::buildDetectors() -I- Use VMC geometry" << endl;
00060   SetCurrentDetectorBuilder(this);
00061   const VolumeMap_t TpcVolumes[] = {
00062     //  {"TPCE","the TPC system in STAR","HALL_1/CAVE_1/TPCE_1","",""},
00063     //  {"TPCW","the TPC supporting endcap Wheel","HALL_1/CAVE_1/TPCE_1/TPCW_1-2/*","",""},
00064     //  {"TPEA","one endcap placed in TPC","HALL_1/CAVE_1/TPCE_1/TPEA_1-2/*","",""},
00065     //  {"TPCM","the Central Membrane placed in TPC","HALL_1/CAVE_1/TPCE_1/TPCM_1","",""},
00066     //  {"TOFC","outer field cage - fill it with insulating gas already","HALL_1/CAVE_1/TPCE_1/TOFC_1/*","",""},
00067     {"TIFC","Inner Field Cage","HALL_1/CAVE_1/TPCE_1/TIFC_1","",""},
00068     {"TIFC","Inner Field Cage","HALL_1/CAVE_1/TpcRefSys_1/TPCE_1/TIFC_1","",""},
00069     {"TOFC","Inner Field Cage","HALL_1/CAVE_1/TPCE_1/TOFC_1","",""},
00070     {"TOFC","Inner Field Cage","HALL_1/CAVE_1/TpcRefSys_1/TPCE_1/TOFC_1","",""},
00071     //  {"TPGV","the Gas Volume placed in TPC","HALL_1/CAVE_1/TPCE_1/TPGV_1-2/*","",""},
00072     //  {"TPSS","a division of gas volume corresponding to a supersectors","HALL_1/CAVE_1/TPCE_1/TPGV_1-2/TPSS_1-12/*","",""},
00073     {"TPAD","inner pad row","HALL_1/CAVE_1/TPCE_1/TPGV_%d/TPSS_%d/TPAD_%d","tpc",""},// <+++
00074     {"TPA1","outer pad row","HALL_1/CAVE_1/TPCE_1/TPGV_%d/TPSS_%d/TPA1_%d","tpc",""},
00075     {"tpad","all pad rows","/HALL_1/CAVE_1/TpcRefSys_1/TPCE_1/TpcSectorWhole_%d/TpcGas_1/TpcPadPlane_%d/tpad_%d","tpc"} // VMC
00076   };
00077   Bool_t newRefSystem = kTRUE;
00078   TString path("HALL_1/CAVE_1/TpcRefSys_1/TPCE_1");
00079   if (! gGeoManager->cd(path)) newRefSystem = kFALSE;
00080 
00081   // change to +1 instead of +2 to remove the ofc.
00082   UInt_t nRows = St_tpcPadPlanesC::instance()->numberOfRows();// Only sensitive detectors
00083   setNRows(nRows);
00084   UInt_t row;
00085 #if 0
00086   Int_t NoStiSectors = 24;
00087 #else
00088   Int_t NoStiSectors = 12;
00089 #endif
00090   for (row = 0; row < nRows; row++) setNSectors(row,NoStiSectors);
00091   // Get Materials
00092   TGeoVolume *volT = gGeoManager->GetVolume("TPAD"); 
00093   if (! volT) volT = gGeoManager->GetVolume("tpad"); 
00094   assert (volT);
00095   TGeoMaterial *mat = volT->GetMaterial(); assert(mat); if (debug>1) mat->Print();
00096   Double_t PotI = StiVMCToolKit::GetPotI(mat); if (debug>1) cout << "PotI " << PotI << endl;
00097   _gasMat = add(new StiMaterial(mat->GetName(),
00098                                 mat->GetZ(),
00099                                 mat->GetA(),
00100                                 mat->GetDensity(),
00101                                 mat->GetDensity()*mat->GetRadLen(),
00102                                 PotI));
00103   Double_t ionization = _gasMat->getIonization();
00104   StiElossCalculator *gasElossCalculator =  new StiElossCalculator(_gasMat->getZOverA(), ionization*ionization,
00105                                                                    _gasMat->getA(), _gasMat->getZ(), _gasMat->getDensity());
00106   StDetectorDbTpcRDOMasks *s_pRdoMasks = StDetectorDbTpcRDOMasks::instance();
00107   StiPlanarShape *pShape;
00108   //Active TPC padrows
00109   //  Double_t radToDeg = 180./3.1415927;
00110   StTpcCoordinateTransform transform(gStTpcDb);
00111   StMatrixD  local2GlobalRotation;
00112   StMatrixD  unit(3,3,1);
00113   StThreeVectorD RowPosition;
00114   UInt_t nInnerPadrows = St_tpcPadPlanesC::instance()->numberOfInnerRows();
00115   for(row = 0; row<45; row++)    {
00116     //Nominal pad row information.
00117     // create properties shared by all sectors in this padrow
00118     float fRadius = St_tpcPadPlanesC::instance()->radialDistanceAtRow(row+1);
00119     TString name(Form("Tpc/Padrow_%d", row));
00120     pShape = new StiPlanarShape;
00121     if (!pShape)
00122       throw runtime_error("StiTpcDetectorBuilder::buildDetectors() - FATAL - pShape==0||ofcShape==0");
00123     Double_t dZ = 0;
00124     if(row < nInnerPadrows) {
00125       pShape->setThickness(St_tpcPadPlanesC::instance()->innerSectorPadLength());
00126       dZ = St_tpcPadPlanesC::instance()->innerSectorPadPlaneZ();
00127     }
00128     else {
00129       pShape->setThickness(St_tpcPadPlanesC::instance()->outerSectorPadLength());
00130       dZ = St_tpcPadPlanesC::instance()->outerSectorPadPlaneZ();
00131     }
00132     pShape->setHalfDepth(dZ*24/NoStiSectors);
00133     pShape->setHalfWidth(St_tpcPadPlanesC::instance()->PadPitchAtRow(row+1) * St_tpcPadPlanesC::instance()->numberOfPadsAtRow(row+1) / 2.);
00134     pShape->setName(name.Data()); if (debug>1) cout << *pShape << endl;
00135     for(UInt_t sector = 0; sector<getNSectors(); sector++) {
00136       //Retrieve position and orientation of the TPC pad rows from the database.
00137       StTpcLocalSectorDirection  dirLS[3];
00138       dirLS[0] = StTpcLocalSectorDirection(1.,0.,0.,sector+1,row+1);
00139       dirLS[1] = StTpcLocalSectorDirection(0.,1.,0.,sector+1,row+1);
00140       dirLS[2] = StTpcLocalSectorDirection(0.,0.,1.,sector+1,row+1);
00141       local2GlobalRotation = unit;
00142       for (Int_t i = 0; i < 3; i++) {
00143         //      if (debug>1) cout << "dirLS\t" << dirLS[i] << endl;
00144 #ifndef TPC_IDEAL_GEOM
00145         StTpcLocalDirection        dirL;
00146         StTpcLocalSectorAlignedDirection  dirLSA;
00147         transform(dirLS[i],dirLSA);//   if (debug>1) cout << "dirLSA\t" << dirLSA << endl;
00148         transform(dirLSA,dirL);    //   if (debug>1) cout << "dirL\t" << dirL << endl;
00149         StGlobalDirection          dirG;
00150         transform(dirL,dirG);//      if (debug>1) cout << "dirG\t" << dirG << endl;
00151 #else
00152         StTpcLocalDirection  dirG;
00153         transform(dirLS[i],dirG);
00154 #endif
00155         local2GlobalRotation(i+1,1) = dirG.position().x();
00156         local2GlobalRotation(i+1,2) = dirG.position().y();
00157         local2GlobalRotation(i+1,3) = dirG.position().z();
00158       }
00159       //      if (debug>1) cout << "Local2GlobalRotation = " << local2GlobalRotation << endl;
00160       Double_t y  = transform.yFromRow(row+1);
00161       StTpcLocalSectorCoordinate  lsCoord(0., y, dZ, sector+1, row+1);// if (debug>1) cout << lsCoord << endl;
00162 #ifndef TPC_IDEAL_GEOM
00163       StTpcLocalSectorAlignedCoordinate lsCoordA;
00164       transform(lsCoord,lsCoordA);//                       if (debug>1) cout << lsCoordA << endl;
00165       StGlobalCoordinate  gCoord;
00166       transform(lsCoordA, gCoord);//                       if (debug>1) cout << gCoord << endl;
00167 #else  // Ideal geom
00168       StTpcLocalCoordinate gCoord;
00169       transform(lsCoord, gCoord);
00170 #endif
00171       //unit vector normal to the pad plane
00172       StThreeVectorD centerVector(gCoord.position().x(),gCoord.position().y(),gCoord.position().z());
00173       StThreeVectorD normalVector(local2GlobalRotation(2,1),
00174                                   local2GlobalRotation(2,2),
00175                                   local2GlobalRotation(2,3));
00176       Double_t prod = centerVector*normalVector;
00177       if (prod < 0) normalVector *= -1;
00178       Double_t phi  = centerVector.phi();
00179       Double_t phiD = normalVector.phi();
00180       Double_t r = centerVector.perp();
00181       StiPlacement *pPlacement = new StiPlacement;
00182       Double_t zc = 0;
00183       if (NoStiSectors != 12) zc = centerVector.z();
00184       pPlacement->setZcenter(zc);
00185       pPlacement->setLayerRadius(fRadius);
00186       pPlacement->setLayerAngle(phi);
00187       pPlacement->setRegion(StiPlacement::kMidRapidity);
00188       pPlacement->setNormalRep(phiD, r*TMath::Cos(phi-phiD), r*TMath::Sin(phi-phiD));
00189       name = Form("Tpc/Padrow_%d/Sector_%d", row, sector);
00190       // fill in the detector object and save it in our vector
00191       StiDetector *pDetector = _detectorFactory->getInstance();
00192       pDetector->setName(name.Data());
00193       pDetector->setIsOn(true);
00194 
00195       Int_t iRdo  = s_pRdoMasks->rdoForPadrow(row+1);
00196       Bool_t west = s_pRdoMasks->isOn(sector+1, iRdo);
00197       Bool_t east = s_pRdoMasks->isOn( 24-(sector+1)%12, iRdo);
00198 
00199 #if 0
00200       if (row==12) {
00201         east = false;
00202         west = false;
00203       }
00204 #endif
00205       if (west) {
00206         Int_t sec = sector+1;
00207         west = St_tpcAnodeHVavgC::instance()->livePadrow(sec,row+1) &&
00208                St_tpcPadGainT0C::instance()->livePadrow(sec,row+1);
00209       }
00210       if (east) {
00211         Int_t sec = 24-(sector+1)%12;
00212         east = St_tpcAnodeHVavgC::instance()->livePadrow(sec,row+1) &&
00213                St_tpcPadGainT0C::instance()->livePadrow(sec,row+1);
00214       }
00215       pDetector->setIsActive(new StiTpcIsActiveFunctor(_active,west,east));
00216       pDetector->setIsContinuousMedium(true);
00217       pDetector->setIsDiscreteScatterer(false);
00218       pDetector->setMaterial(_gasMat);
00219       pDetector->setGas(_gasMat);
00220       pDetector->setShape(pShape);
00221       pDetector->setPlacement(pPlacement);
00222       if (row < nInnerPadrows)
00223         pDetector->setHitErrorCalculator(StiTpcInnerHitErrorCalculator::instance());
00224       else
00225         pDetector->setHitErrorCalculator(StiTpcOuterHitErrorCalculator::instance());
00226       pDetector->setElossCalculator(gasElossCalculator);
00227       pDetector->setKey(1,row);
00228       pDetector->setKey(2,sector);
00229       add(row,sector,pDetector); if (debug>1) cout << *pDetector << endl;
00230     }// for sector
00231   }// for row
00232   for (Int_t i = 0; i < 4; i++) {
00233     gGeoManager->RestoreMasterVolume();
00234     gGeoManager->CdTop();
00235     TGeoNode *nodeT = gGeoManager->GetCurrentNode();
00236     path = TpcVolumes[i].path;
00237     if (  newRefSystem && ! path.Contains("TpcRefSys")) continue;
00238     if (! newRefSystem &&   path.Contains("TpcRefSys")) continue;
00239     if (! gGeoManager->cd(path)) continue;
00240     nodeT = gGeoManager->GetCurrentNode();
00241     if (! nodeT) continue;
00242     path = gGeoManager->GetPath();
00243     StiVMCToolKit::LoopOverNodes(nodeT, path, TpcVolumes[i].name, MakeAverageVolume);
00244   }
00245   cout << "StiTpcDetectorBuilder::buildDetectors() -I- Done" << endl;
00246 }