// ------------------------------------------------------------------------- // ----- CbmStsHitProducerIdeal source file ----- // ----- Created 10/01/06 by V. Friese ----- // ------------------------------------------------------------------------- #include "TClonesArray.h" #include "TArrayD.h" #include "TVector2.h" #include "TGeoManager.h" #include "CbmRootManager.h" #include "PndMvdStripHitProducer.h" #include "PndMvdMCPoint.h" #include "CbmRunAna.h" #include "CbmRuntimeDb.h" #include "CbmGeoNode.h" #include "CbmRuntimeDb.h" #include "CbmGeoNode.h" #include "CbmGeoVector.h" #include "PndStringVector.h" #include "PndMvdCalcStrip.h" #include "PndMvdDigiStrip.h" // ----- Default constructor ------------------------------------------- PndMvdStripHitProducer::PndMvdStripHitProducer() : CbmTask("MVD Strip Hit Producer") { fBranchName = "MVDPoint"; // stripHits = 0; ftopPitch = 0.; fbotPitch = 0.; forient = 0.; fskew = 0.; fTopAnchor = TVector2(0,0); fBotAnchor = TVector2(0,0); fnrTopFE = 0; fnrBotFE = 0; fnrFECh = 0; fthreshold = 0.; fnoise = 0.; fOverrideParams = false; // fHitArray = new TClonesArray("PndMvdHit"); // fStripArray = new TClonesArray("PndMvdStripHit"); } // ------------------------------------------------------------------------- PndMvdStripHitProducer::PndMvdStripHitProducer(Double_t topPitch, Double_t botPitch, Double_t ori, Double_t skew, TVector2 topAnchor, TVector2 botAnchor, Int_t nrTopFE, Int_t nrBotFE, Int_t nrFECh, Double_t threshold, Double_t noise) : CbmTask("MVD Strip Hit Producer") { fBranchName = "MVDPoint"; ftopPitch = topPitch; fbotPitch = botPitch; forient = ori; fskew = skew; fTopAnchor = topAnchor; fBotAnchor = botAnchor; fnrTopFE = nrTopFE; fnrBotFE = nrBotFE; fnrFECh = nrFECh; fthreshold = threshold; fnoise = noise; fOverrideParams = true; // stripHits = 0; std::cout << "MVD Strip Hit Producer initiated" << std::endl; } // ------------------------------------------------------------------------- // ----- Destructor ---------------------------------------------------- PndMvdStripHitProducer::~PndMvdStripHitProducer() { delete fGeoH; } // ------------------------------------------------------------------------- // ----- Initialization of Parameter Containers ------------------------- void PndMvdStripHitProducer::SetParContainers() { // Get Base Container CbmRunAna* ana = CbmRunAna::Instance(); CbmRuntimeDb* rtdb=ana->GetRuntimeDb(); fDigiPar = (PndMvdStripDigiPar*)(rtdb->getContainer("PndMvdStripDigiPar")); fGeoMappingPar = (PndMvdGeoMappingPar*)(rtdb->getContainer("PndMvdGeoMappingPar")); } InitStatus PndMvdStripHitProducer::ReInit() { CbmRunAna* ana = CbmRunAna::Instance(); CbmRuntimeDb* rtdb=ana->GetRuntimeDb(); fDigiPar=(PndMvdStripDigiPar*)(rtdb->getContainer("PndMvdStripDigiPar")); fGeoMappingPar = (PndMvdGeoMappingPar*)(rtdb->getContainer("PndMvdGeoMappingPar")); return kSUCCESS; } // ----- Public method Init -------------------------------------------- InitStatus PndMvdStripHitProducer::Init() { CbmRunAna* ana = CbmRunAna::Instance(); CbmRootManager* ioman = CbmRootManager::Instance(); fGeoH = new PndMvdGeoHandling(gGeoManager); if ( ! ioman ) { std::cout << "-E- PndMvdStripHitProducer::Init: " << "RootManager not instantiated!" << std::endl; return kFATAL; } fPointArray = (TClonesArray*) ioman->GetObject(fBranchName); if ( ! fPointArray ) { std::cout << "-W- PndMvdStripHitProducer::Init: " << "No MVDPoint array!" << std::endl; return kERROR; } // Create and register output array // fHitArray = new TClonesArray("PndMvdHit"); // ioman->Register("MVDHit", "MVD", fHitArray, kTRUE); // Create and register output array fStripArray = new TClonesArray("PndMvdDigiStrip"); ioman->Register("MVDStripHit", "MVD", fStripArray, kTRUE); // Create and register parameter array // fStripArray = new TClonesArray("PndMvdDigiPar"); // ioman->Register("MVDDigiParam", "MVD", fDigiPar, kTRUE); // SetParContainers(); std::cout << "-I- PndMvdStripHitProducer: Initialisation successfull" << std::endl; std::cout << "-I- Test PndMvdStripHitProducer: detIdPairArray " << (fGeoMappingPar->getPairs())->GetEntries() << std::endl; if (!fDigiPar){ std::cout<<"-E- PndMvdStripHitProducer: DigiPar Container does not exist!"<SetTopPitch(ftopPitch); fDigiPar->SetBotPitch(fbotPitch); fDigiPar->SetSkew(fskew); fDigiPar->SetOrient(forient); fDigiPar->SetTopAnchor(fTopAnchor); fDigiPar->SetBotAnchor(fBotAnchor); fDigiPar->SetNrFECh(fnrFECh); fDigiPar->SetNrTopFE(fnrTopFE); fDigiPar->SetNrBotFE(fnrBotFE); fDigiPar->SetThreshold(fthreshold); fDigiPar->SetNoise(fnoise); fDigiPar->setChanged(); fDigiPar->setInputVersion(ana->GetRunId(),1); } fDigiPar->Print(); return kSUCCESS; } // ------------------------------------------------------------------------- // ----- Public method Exec -------------------------------------------- void PndMvdStripHitProducer::Exec(Option_t* opt) { // Reset output array // if ( ! fHitArray ) // Fatal("Exec", "No HitArray"); //fHitArray->Clear(); fStripArray->Clear(); // fFeStripArray->Clear(); // Declare some variables PndMvdMCPoint *point = NULL; Int_t detID = 0, // Detector ID trackID = 0; // Track index // Loop over PndMvdMCPoints Int_t nPoints = fPointArray->GetEntriesFast(); if (fVerbose > 0){ std::cout<<" Nr of Points: "<GetTopPitch(), fDigiPar->GetOrient(), fDigiPar->GetNrTopFE() * fDigiPar->GetNrFECh() , fDigiPar->GetNrFECh(), fDigiPar->GetTopAnchor(), fDigiPar->GetThreshold(), fDigiPar->GetNoise() ); PndMvdCalcStrip StripCalcBot( fDigiPar->GetBotPitch(), fDigiPar->GetOrient() + fDigiPar->GetSkew(), fDigiPar->GetNrBotFE() * fDigiPar->GetNrFECh(), fDigiPar->GetNrFECh(), fDigiPar->GetBotAnchor(), fDigiPar->GetThreshold(), fDigiPar->GetNoise() ); for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) { point = (PndMvdMCPoint*) fPointArray->At(iPoint); if (!((fGeoH->GetPath(point->GetDetName())).Contains("Strip"))) { if (fVerbose > 1) std::cout<<"detector name does not contain 'Strip'"< 2) std::cout<<" DetName : "<GetPath(point->GetDetName())<GetPath(point->GetDetName())).Contains("Disk")) { if (fVerbose > 1) std::cout<<"detector name does contain 'Disk'"< 2) std::cout<<" DetName : "<GetPath(point->GetDetName())<GetPath(point->GetDetName())).Contains("SideChips")) continue; if (fVerbose > 1){ std::cout<<"***** Strip Hit ******"<GetPath(point->GetDetName())< 1){ std::cout << "****Global Point: " << std::endl; point->Print(""); } CbmGeoVector posInL, posOutL, meanPos, meanPosL; GetLocalHitPoints(point, posInL, posOutL); if (fVerbose > 1){ std::cout << "Local Hits: " << std::endl; posInL.Print(); posOutL.Print(); } meanPosL = posInL + posOutL; meanPosL /= 2; if (fVerbose > 1){ meanPosL.Print(); std::cout << "Energy: " << point->GetEnergyLoss() << std::endl; } TVector3 size = GetSensorDimensions(fGeoH->GetPath(point->GetDetName()).Data()); if (fVerbose > 2) std::cout<< "Sensor Size : x="< topStrips = StripCalcTop.GetStrips(posInL.getX(), posInL.getY(), posInL.getZ(), posOutL.getX(), posOutL.getY(), posOutL.getZ(), point->GetEnergyLoss()*1E9); if (topStrips.size() != 0) { // stripHits += topStrips.size(); if (fVerbose > 1) std::cout << "SensorStrips: " << std::endl; trackID = point->GetTrackID(); detID = point->GetDetectorID(); // if (fVerbose > 2) std::cout<<"Number of Hits: "<::const_iterator kit=topStrips.begin(); kit!= topStrips.end(); ++kit) { new ((*fStripArray)[iStrip]) PndMvdDigiStrip(iPoint,detID,fGeoH->GetPath(point->GetDetName()), StripCalcTop.CalcFEfromStrip(kit->GetIndex()), StripCalcTop.CalcChannelfromStrip(kit->GetIndex()),kit->GetCharge()); dynamic_cast(fStripArray->operator[](iStrip))->SetMCID(trackID); if (fVerbose > 1) std::cout << *kit << std::endl; iStrip++; } } // Bottom Side if (fVerbose > 1) std::cout << "Bottom Side: " << std::endl; StripCalcBot.SetVerboseLevel(fVerbose); std::vector botStrips = StripCalcBot.GetStrips(posInL.getX(), posInL.getY(), posInL.getZ(), posOutL.getX(), posOutL.getY(), posOutL.getZ(), point->GetEnergyLoss()*1E9); if (botStrips.size() != 0){ // stripHits += botStrips.size(); if (fVerbose > 1) std::cout << " SensorStrips: " << std::endl; trackID = point->GetTrackID(); detID = point->GetDetectorID(); for(std::vector::const_iterator kit=botStrips.begin(); kit!= botStrips.end(); ++kit) { new ((*fStripArray)[iStrip]) PndMvdDigiStrip(iPoint,detID,fGeoH->GetPath(point->GetDetName()), StripCalcBot.CalcFEfromStrip(kit->GetIndex()) + fDigiPar->GetNrTopFE(), StripCalcBot.CalcChannelfromStrip(kit->GetIndex()),kit->GetCharge()); dynamic_cast( (*fStripArray)[iStrip])->SetMCID(trackID); if (fVerbose > 1) std::cout << *kit << std::endl; iStrip++; } } } // Loop over MCPoints // Event summary std::cout << "-I- PndMvdStripHitProducer: " << nPoints << " PndMvdMCPoints, " << iStrip << " Hits created."<< std::endl; } void PndMvdStripHitProducer::GetLocalHitPoints(PndMvdMCPoint* myPoint, CbmGeoVector& myHitIn, CbmGeoVector& myHitOut) { if (fVerbose > 1) std::cout << "GetLocalHitPoints" << std::endl; TGeoHMatrix trans = GetTransformation(fGeoH->GetPath(myPoint->GetDetName()).Data()); Double_t posIn[3]; Double_t posOut[3]; Double_t posInLocal[3]; Double_t posOutLocal[3]; posIn[0] = myPoint->GetX(); posIn[1] = myPoint->GetY(); posIn[2] = myPoint->GetZ(); posOut[0] = myPoint->GetXOut(); posOut[1] = myPoint->GetYOut(); posOut[2] = myPoint->GetZOut(); if (fVerbose > 1){ for (Int_t i = 0; i < 3; i++) std::cout << "posIn "<< i << ": " << posIn[i] << std::endl; trans.Print(""); } trans.MasterToLocal(posIn, posInLocal); trans.MasterToLocal(posOut, posOutLocal); if (fVerbose > 1) { for (Int_t i = 0; i < 3; i++){ std::cout << "posInLocal "<< i << ": " << posInLocal[i] << std::endl; std::cout << "posOutLocal "<< i << ": " << posOutLocal[i] << std::endl; } } //posIn/OutLocal have the center of the coordinate system in the center of the shape //typically sensors have their coordinate system centered at the lower left corner TVector3 offset = GetSensorDimensions(fGeoH->GetPath(myPoint->GetDetName()).Data()); posInLocal[0] += offset.x(); posInLocal[1] += offset.y(); //posInLocal[2] += offset.z(); posOutLocal[0] += offset.x(); posOutLocal[1] += offset.y(); //posOutLocal[2] += offset.z(); myHitIn.setVector(posInLocal); myHitOut.setVector(posOutLocal); } TGeoHMatrix PndMvdStripHitProducer::GetTransformation(std::string detName) const { gGeoManager->cd(detName.c_str()); TGeoHMatrix* transMat = gGeoManager->GetCurrentMatrix(); if (fVerbose > 1) transMat->Print(""); return *transMat; } TVector3 PndMvdStripHitProducer::GetSensorDimensions(std::string detName) const { gGeoManager->cd(detName.c_str()); TGeoVolume* actVolume = gGeoManager->GetCurrentVolume(); TGeoBBox* actBox = (TGeoBBox*)(actVolume->GetShape()); TVector3 result; result.SetX(actBox->GetDX()); result.SetY(actBox->GetDY()); result.SetZ(actBox->GetDZ()); //result.Dump(); return result; } // ------------------------------------------------------------------------- ClassImp(PndMvdStripHitProducer)