// ------------------------------------------------------------------------- // ----- CbmStsHitProducerIdeal source file ----- // ----- Created 10/01/06 by V. Friese ----- // ------------------------------------------------------------------------- #include "TClonesArray.h" #include "TArrayD.h" #include "TGeoManager.h" #include "CbmRootManager.h" #include "MvdStripHitProducer.h" #include "MvdHit.h" #include "MvdPoint.h" #include "CbmRunAna.h" #include "CbmRuntimeDb.h" #include "CbmGeoNode.h" #include "CbmRuntimeDb.h" #include "CbmGeoNode.h" #include "CbmGeoVector.h" #include "StringVector.h" #include "MvdCalcStrip.h" #include "MvdCalcFePixel.h" #include "MvdDigiPixel.h" // ----- Default constructor ------------------------------------------- MvdStripHitProducer::MvdStripHitProducer() : CbmTask("MVD Strip Hit Producer") { fBranchName = "MVDPoint"; stripHits = 0; // fHitArray = new TClonesArray("MvdHit"); // fStripArray = new TClonesArray("MvdPixelHit"); } // ------------------------------------------------------------------------- MvdStripHitProducer::MvdStripHitProducer(Double_t lx, Double_t ly, Double_t skew, Double_t threshold, Double_t noise) : CbmTask("MVD Strip Hit Producer") { fBranchName = "MVDPoint"; flx = lx; fly = ly; fskew = skew; fthreshold = threshold; fnoise = noise; stripHits = 0; std::cout << "MVD Strip Hit Producer initiated" << std::endl; } // ------------------------------------------------------------------------- // ----- Destructor ---------------------------------------------------- MvdStripHitProducer::~MvdStripHitProducer() { } // ------------------------------------------------------------------------- // ----- Initialization of Parameter Containers ------------------------- void MvdStripHitProducer::SetParContainers() { // Get Base Container CbmRunAna* ana = CbmRunAna::Instance(); CbmRuntimeDb* rtdb=ana->GetRuntimeDb(); fDigiPar = (MvdDigiPar*)(rtdb->getContainer("MvdDigiPar")); fGeoMappingPar = (MvdGeoMappingPar*)(rtdb->getContainer("MvdGeoMappingPar")); } InitStatus MvdStripHitProducer::ReInit() { CbmRunAna* ana = CbmRunAna::Instance(); CbmRuntimeDb* rtdb=ana->GetRuntimeDb(); fDigiPar=(MvdDigiPar*)(rtdb->getContainer("MvdDigiPar")); fGeoMappingPar = (MvdGeoMappingPar*)(rtdb->getContainer("MvdGeoMappingPar")); return kSUCCESS; } // ----- Public method Init -------------------------------------------- InitStatus MvdStripHitProducer::Init() { CbmRunAna* ana = CbmRunAna::Instance(); CbmRootManager* ioman = CbmRootManager::Instance(); if ( ! ioman ) { std::cout << "-E- MvdStripHitProducer::Init: " << "RootManager not instantiated!" << std::endl; return kFATAL; } fPointArray = (TClonesArray*) ioman->GetObject(fBranchName); if ( ! fPointArray ) { std::cout << "-W- MvdStripHitProducer::Init: " << "No MVDPoint array!" << std::endl; return kERROR; } // Create and register output array fHitArray = new TClonesArray("MvdHit"); ioman->Register("MVDHit", "MVD", fHitArray, kTRUE); // Create and register output array fStripArray = new TClonesArray("MvdDigiPixel"); ioman->Register("MVDStripHit", "MVD", fStripArray, kTRUE); fFeStripArray = new TClonesArray("MvdDigiPixel"); ioman->Register("MVDFePixelHit", "MVD", fFeStripArray, kTRUE); std::cout << "-I- MvdStripHitProducer: Intialisation successfull" << std::endl; std::cout << "-I- Test: " << fGeoMappingPar->detIdPairArray->GetEntries() << std::endl; fDigiPar->dimX = flx; fDigiPar->dimY = fly; fDigiPar->skew = fskew; fDigiPar->threshold = fthreshold; fDigiPar->noise = fnoise; fDigiPar->setInputVersion(ana->GetRunId(),1); fDigiPar->setChanged(); return kSUCCESS; } // ------------------------------------------------------------------------- // ----- Public method Exec -------------------------------------------- void MvdStripHitProducer::Exec(Option_t* opt) { // Reset output array if ( ! fHitArray ) Fatal("Exec", "No HitArray"); fHitArray->Clear(); fStripArray->Clear(); fFeStripArray->Clear(); // Declare some variables MvdPoint *point = NULL; Int_t detID = 0, // Detector ID trackID = 0; // Track index // Loop over MvdPoints Int_t nPoints = fPointArray->GetEntriesFast(); Int_t iStrip = 0; Int_t iFeStrip = 0; for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) { point = (MvdPoint*) fPointArray->At(iPoint); if ( ! point){ std::cout<< "No Point!" << std::endl; continue; } if (fVerbose > 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; } MvdCalcStrip StripCalc(flx, fly, fskew, fthreshold, fnoise); // Calculate a cluster of Strips fired std::vector myStrips = StripCalc.GetStrips(posInL.getX(), posInL.getY(), posInL.getZ(), posOutL.getX(), posOutL.getY(), posOutL.getZ(), point->GetEnergyLoss()*1E9); if (myStrips.size() == 0) { std::cout << "Deposited charge below threshold" << std::endl; continue; } // else { stripHits += myStrips.size(); if (fVerbose > 1) std::cout << "SensorStrips: " << std::endl; for(uint i = 0; i < myStrips.size(); i++) { myStrips[i].SetDetName(point->GetDetName().Data()); if (fVerbose > 1) std::cout << myStrips[i] << std::endl; } // TODO read Parameters from Database // Calculate channel numbers MvdCalcFePixel FeCalc(76, 84, 10); // this is a cluster, ideal finding std::vector myFeStrips = FeCalc.CalcFEHits(myStrips); if (fVerbose > 1) { std::cout << "FeStrips: " << myFeStrips.size() << std::endl; for(uint i = 0; i < myFeStrips.size(); i++){ std::cout << myFeStrips[i] << std::endl; } } // MCTrack ID trackID = point->GetTrackID(); detID = point->GetDetectorID(); //std::cout << "Local Coordinates from Point:" << std::endl; //std::cout << point->GetXInLocal()<< " " << point->GetYInLocal() // << " " << point->GetZInLocal() << std::endl; //std::cout << point->GetXOutLocal() << " " << point->GetYOutLocal() // << " " << point->GetZOutLocal() << std::endl; MvdHit hitCoord = CalcGlobalPoint(myStrips); TVector3 MCCoord(point->GetX() + point->GetXOut(), point->GetY()+point->GetYOut(), point->GetZ()+point->GetZOut()); MCCoord *= 0.5; TVector3 pos = hitCoord.GetPosition(); if (fVerbose > 1) { std::cout << "Calculated Hit: " << pos.X() << " " << pos.Y() << " " << pos.Z() << " " << hitCoord.GetCharge() << "\n" << "MC Hit:" << MCCoord.X() << " " << MCCoord.Y() << " " << MCCoord.Z() << "\n" << "Difference: " << (MCCoord - pos).Mag() << "\n" << "ChargeDifference: " << (hitCoord.GetCharge() * 3.61 - point->GetEnergyLoss()*1E9) << std::endl; } hitCoord.SetTrackID(trackID); hitCoord.SetDetectorID(detID); hitCoord.SetDetName(point->GetDetName()); hitCoord.SetRefIndex(iPoint); new ((*fHitArray)[iPoint])MvdHit(hitCoord); for (uint iPix = 0; iPix < myStrips.size(); iPix++) { //if (myStrips[iPix].GetCharge() > 0) new ((*fStripArray)[iStrip++]) MvdDigiPixel(iPoint, detID, point->GetDetName(), myStrips[iPix].GetFE(), myStrips[iPix].GetCol(), myStrips[iPix].GetRow(), myStrips[iPix].GetCharge()); } for (uint iPix = 0; iPix < myFeStrips.size(); iPix++) { //if (myStrips[iPix].GetCharge() > 0) new ((*fFeStripArray)[iFeStrip++]) MvdDigiPixel(iPoint, detID, point->GetDetName(), myFeStrips[iPix].GetFE(), myFeStrips[iPix].GetCol(), myFeStrips[iPix].GetRow(), myFeStrips[iPix].GetCharge()); } // } // } // Loop over MCPoints // Event summary std::cout << "-I- MvdStripHitProducer: " << nPoints << " MvdPoints, " << stripHits << " Hits created." << " " << iStrip<< std::endl; } void MvdStripHitProducer::GetLocalHitPoints(MvdPoint* myPoint, CbmGeoVector& myHitIn, CbmGeoVector& myHitOut) { if (fVerbose > 1) std::cout << "GetLocalHitPoints" << std::endl; TGeoHMatrix trans = GetTransformation(myPoint->GetDetName().Data()); double posIn[3]; double posOut[3]; double posInLocal[3]; double 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 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 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(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); } MvdHit MvdStripHitProducer::CalcGlobalPoint(std::vector strips) { Double_t col = 0, row = 0, charge = 0; Double_t colB = 0, rowB = 0; int count = 0; if (strips.size() == 1){ if (strips[0].GetCharge() > 0){ col = strips[0].GetCol(); if (col < 0) col -= 0.5; else col += 0.5; row = strips[0].GetRow(); if (row < 0) row -= 0.5; else row += 0.5; count = 1; charge = strips[0].GetCharge(); } colB = col; rowB = row; } else { //std::cout << "Multiple Hits!" << std::endl; for (uint i = 0; i < strips.size(); i++){ //std::cout << "ActCol / Row" << col << " " << row << " added Col/Row " << strips[i].GetCol() << " " << strips[i].GetRow() << std::endl; if (strips[i].GetCharge() > 0){ col += (strips[i].GetCol()*strips[i].GetCharge()); row += (strips[i].GetRow()*strips[i].GetCharge()); charge += strips[i].GetCharge(); colB += (strips[i].GetCol()); rowB += (strips[i].GetRow()); count++; } } if (count > 0){ colB /= count; rowB /= count; if (charge > 0){ col /= charge; row /= charge; } else { col = colB; row = rowB; } if (col < 0){ col -= 0.5; colB -= 0.5; } else { col += 0.5; colB += 0.5; } if (row < 0){ row -= 0.5; rowB -= 0.5; } else { row += 0.5; rowB += 0.5; } } else col = colB = row = rowB = 0; } if (fVerbose > 1){ std::cout << "Col: " << col << " Row: " << row << std::endl; std::cout << "ColB: " << colB << " RowB: " << rowB << std::endl; } double resultD[3]; resultD[0] = col * fly; resultD[1] = row * flx; resultD[2] = 0; double resultFinal[3]; TVector3 offset = GetSensorDimensions(strips[0].GetDetName()); resultD[0] -= offset.x(); resultD[1] -= offset.y(); //resultD[2] -= offset.z(); TGeoHMatrix trans = GetTransformation(strips[0].GetDetName()); //std::cout << "Transformation for: " << strips[0].GetDetName() << std::endl; //trans.Print(""); trans.LocalToMaster(resultD, resultFinal); //result.setXYZ(resultFinal[0], resultFinal[1], resultFinal[2]); TVector3 dpos(flx/sqrt(12.0),fly/sqrt(12.0),0); TVector3 pos(resultFinal[0], resultFinal[1], resultFinal[2]); return (MvdHit(0, 0, strips[0].GetDetName(), pos, dpos, 0, charge, count)); } //CbmGeoTransform MvdStripHitProducer::GetTransformation(std::string detName) TGeoHMatrix MvdStripHitProducer::GetTransformation(std::string detName) { gGeoManager->cd(detName.c_str()); TGeoHMatrix* transMat = gGeoManager->GetCurrentMatrix(); if (fVerbose > 1) transMat->Print(""); return *transMat; } TVector3 MvdStripHitProducer::GetSensorDimensions(std::string detName) { 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(MvdStripHitProducer)