/* * CbmRichGeoManager.cxx * * Created on: Dec 16, 2015 * Author: slebedev */ #include "CbmRichGeoManager.h" #include #include "FairLogger.h" #include "TGeoNode.h" #include "TGeoManager.h" #include "TGeoSphere.h" #include "FairLogger.h" using namespace std; CbmRichGeoManager::CbmRichGeoManager() : fGP(NULL) { InitGeometry(); } void CbmRichGeoManager::InitGeometry() { LOG(INFO) << "CbmRichGeoManager::InitGeometry" << endl; fGP = new CbmRichRecGeoPar(); //TODO: get refractive index from material fGP->fNRefrac = 1.000446242; DetectGeometryType(); if (fGP->fGeometryType == CbmRichGeometryTypeTwoWings) { InitPmt(); } else if (fGP->fGeometryType == CbmRichGeometryTypeCylindrical) { InitPmtCyl(); } else if (fGP->fGeometryType == CbmRichGeometryTypeNotDefined) { // Fatal("CbmRichGeoManager::InitGeometry()", " Geometry type is CbmRichGeometryTypeNotDefined. Geometry could not be defined automatically."); cout<< endl << endl << "CbmRichGeoManager::InitGeometry(): Geometry type is CbmRichGeometryTypeNotDefined. Geometry could not be defined automatically." << endl << endl; } InitMirror(); fGP->Print(); } void CbmRichGeoManager::DetectGeometryType() { TGeoIterator geoIterator(gGeoManager->GetTopNode()->GetVolume()); geoIterator.SetTopName("/cave_1"); TGeoNode* curNode; geoIterator.Reset(); while ((curNode=geoIterator())) { TString nodeName(curNode->GetName()); TString nodePath; if (curNode->GetVolume()->GetName() == TString("pmt_block_strip")) { fGP->fGeometryType = CbmRichGeometryTypeCylindrical; return; } } geoIterator.Reset(); while ((curNode=geoIterator())) { TString nodeName(curNode->GetName()); TString nodePath; if (curNode->GetVolume()->GetName() == TString("camera_quarter")) { fGP->fGeometryType = CbmRichGeometryTypeTwoWings; return; } } fGP->fGeometryType = CbmRichGeometryTypeNotDefined; } void CbmRichGeoManager::InitPmtCyl() { TGeoIterator geoIterator(gGeoManager->GetTopNode()->GetVolume()); geoIterator.SetTopName("/cave_1"); TGeoNode* curNode; vector xCoord; geoIterator.Reset(); while ((curNode=geoIterator())) { TString nodeName(curNode->GetName()); TString nodePath; if (curNode->GetVolume()->GetName() == TString("pmt_block_strip")) { geoIterator.GetPath(nodePath); const TGeoMatrix* curMatrix = geoIterator.GetCurrentMatrix(); const Double_t* curNodeTr = curMatrix->GetTranslation(); const Double_t* curNodeRot = curMatrix->GetRotationMatrix(); double pmtX = curNodeTr[0]; double pmtY = curNodeTr[1]; double pmtZ = curNodeTr[2]; double rotY = TMath::ASin(-curNodeRot[2]); // around Y double rotZ = TMath::ASin(curNodeRot[1]/TMath::Cos(TMath::ASin(-curNodeRot[2]))); // around Z //double rotX = TMath::ASin(curNodeRot[5]/TMath::Cos(TMath::ASin(-curNodeRot[2]))); // around X double rotX = TMath::ACos(curNodeRot[8]/TMath::Cos(TMath::ASin(-curNodeRot[2]))); // around X fGP->fPmtMap[string(nodePath.Data())].fTheta = rotX; fGP->fPmtMap[string(nodePath.Data())].fPhi = rotY; const TGeoBBox* shape = (const TGeoBBox*)(curNode->GetVolume()->GetShape()); fGP->fPmtMap[string(nodePath.Data())].fWidth = 2. * shape->GetDX(); fGP->fPmtMap[string(nodePath.Data())].fHeight = 2. * shape->GetDY(); fGP->fPmtMap[string(nodePath.Data())].fZ = pmtZ; fGP->fPmtMap[string(nodePath.Data())].fX = pmtX; fGP->fPmtMap[string(nodePath.Data())].fY = pmtY; if (pmtX >=0 && pmtY >=0) { xCoord.push_back(pmtX); } } } std::sort(xCoord.begin(), xCoord.end()); for ( map::iterator it = fGP->fPmtMap.begin(); it != fGP->fPmtMap.end(); it++) { Double_t curX = TMath::Abs(it->second.fX); int pos = -1; //int pos = find(xCoord.begin(), xCoord.end(), curX) - xCoord.begin(); for (int i = 0; i < xCoord.size(); i++) { if (TMath::Abs(curX - xCoord[i]) < 0.1) { pos = i; break; } } it->second.fPmtPositionIndexX = pos; } // We also need to find pmt plane center map mapPmtPlaneMinMax; TString filterNamePixel("pmt_pixel"); geoIterator.Reset(); CbmRichPmtPlaneMinMax pmtPlaneMinMax; while ((curNode=geoIterator())) { TString nodeName(curNode->GetName()); TString nodePath; if (curNode->GetVolume()->GetName() == filterNamePixel) { geoIterator.GetPath(nodePath); string nodePathStr = string(nodePath.Data()); size_t foundIndex1 = nodePathStr.find("pmt_block_strip"); if (foundIndex1 == string::npos) continue; size_t foundIndex2 = nodePathStr.find("/", foundIndex1 + 1); if (foundIndex2 == string::npos) continue; string mapKey = nodePathStr.substr(0, foundIndex2); const TGeoMatrix* curMatrix = geoIterator.GetCurrentMatrix(); const Double_t* curNodeTr = curMatrix->GetTranslation(); double pmtX = curNodeTr[0]; double pmtY = curNodeTr[1]; double pmtZ = curNodeTr[2]; mapPmtPlaneMinMax[mapKey].AddPoint(pmtX, pmtY, pmtZ); } } for ( map::iterator it = fGP->fPmtMap.begin(); it != fGP->fPmtMap.end(); it++) { it->second.fPlaneX = mapPmtPlaneMinMax[it->first].GetMeanX(); it->second.fPlaneY = mapPmtPlaneMinMax[it->first].GetMeanY(); it->second.fPlaneZ = mapPmtPlaneMinMax[it->first].GetMeanZ(); // cout << "name:" << it->first << " strip(x,y,z):" <second.fX << "," << it->second.fY << "," << it->second.fZ << // " pmtPlane(z,y,z)" <second.fPlaneX << "," << it->second.fPlaneY << "," << it->second.fPlaneZ << endl; } // Calculate gap between pmt_block_strips geoIterator.Reset(); double master1[3], master2[3]; while ((curNode=geoIterator())) { TString nodeName(curNode->GetName()); TString nodePath; if (curNode->GetVolume()->GetName() == TString("pmt_block_strip")) { geoIterator.GetPath(nodePath); const TGeoMatrix* curMatrix = geoIterator.GetCurrentMatrix(); const Double_t* curNodeTr = curMatrix->GetTranslation(); const Double_t* curNodeRot = curMatrix->GetRotationMatrix(); double pmtX = curNodeTr[0]; double pmtY = curNodeTr[1]; double pmtZ = curNodeTr[2]; if ( pmtX < 0 || pmtY < 0) continue; const TGeoBBox* shape = (const TGeoBBox*)(curNode->GetVolume()->GetShape()); double loc[3]; if (fGP->fPmtMap[string(nodePath.Data())].fPmtPositionIndexX == 1) { loc[0] = shape->GetDX(); loc[1] = shape->GetDY(); loc[2] = shape->GetDZ(); curMatrix->LocalToMaster(loc, master1); } else if (fGP->fPmtMap[string(nodePath.Data())].fPmtPositionIndexX == 2) { loc[0] = -shape->GetDX(); loc[1] = shape->GetDY(); loc[2] = shape->GetDZ(); curMatrix->LocalToMaster(loc, master2); } } } //cout << master1[0] << " "<< master1[1] << " "<< master1[2]<< endl; //cout << master2[0] << " "<< master2[1] << " "<< master2[2]<< endl; double dist = TMath::Sqrt( (master1[0] - master2[0]) * (master1[0] - master2[0]) + (master1[1] - master2[1]) * (master1[1] - master2[1]) + (master1[2] - master2[2]) * (master1[2] - master2[2]) ); fGP->fPmtStripGap = dist; } void CbmRichGeoManager::InitPmt() { TGeoIterator geoIterator(gGeoManager->GetTopNode()->GetVolume()); TGeoNode* curNode; // PMT plane position\rotation TString filterName_pixel("pmt_pixel"); geoIterator.Reset(); CbmRichPmtPlaneMinMax pmtPlaneMinMax; while ((curNode=geoIterator())) { TString nodeName(curNode->GetName()); TString nodePath; if (curNode->GetVolume()->GetName() == filterName_pixel) { geoIterator.GetPath(nodePath); const TGeoMatrix* curMatrix = geoIterator.GetCurrentMatrix(); const Double_t* curNodeTr = curMatrix->GetTranslation(); const Double_t* curNodeRot = curMatrix->GetRotationMatrix(); double pmtX = curNodeTr[0]; double pmtY = curNodeTr[1]; double pmtZ = curNodeTr[2]; if (pmtX > 0. && pmtY > 0) { //printf ("%08f\t%08f\t%08f\t\n", curNodeTranslation[0], curNodeTranslation[1], curNodeTranslation[2]); double rotY = TMath::ASin(-curNodeRot[2]); // around Y double rotZ = TMath::ASin(curNodeRot[1]/TMath::Cos(TMath::ASin(-curNodeRot[2]))); // around Z //double rotX = TMath::ASin(curNodeRot[5]/TMath::Cos(TMath::ASin(-curNodeRot[2]))); // around X double rotX = TMath::ACos(curNodeRot[8]/TMath::Cos(TMath::ASin(-curNodeRot[2]))); // around X fGP->fPmt.fTheta = rotX; fGP->fPmt.fPhi = rotY; pmtPlaneMinMax.AddPoint(pmtX, pmtY, pmtZ); } } } fGP->fPmt.fPlaneX = pmtPlaneMinMax.GetMeanX(); fGP->fPmt.fPlaneY = pmtPlaneMinMax.GetMeanY(); fGP->fPmt.fPlaneZ = pmtPlaneMinMax.GetMeanZ(); geoIterator.Reset(); while ((curNode=geoIterator())) { TString nodeName(curNode->GetName()); TString nodePath; if (curNode->GetVolume()->GetName() == TString("camera_quarter")) { geoIterator.GetPath(nodePath); const TGeoMatrix* curMatrix = geoIterator.GetCurrentMatrix(); const Double_t* curNodeTr = curMatrix->GetTranslation(); //const Double_t* curNodeRot = curMatrix->GetRotationMatrix(); double pmtX = curNodeTr[0]; double pmtY = curNodeTr[1]; double pmtZ = curNodeTr[2]; if (pmtX > 0. && pmtY > 0) { const TGeoBBox* shape = (const TGeoBBox*)(curNode->GetVolume()->GetShape()); fGP->fPmt.fWidth = shape->GetDX(); fGP->fPmt.fHeight = shape->GetDY(); fGP->fPmt.fZ = pmtZ; fGP->fPmt.fX = pmtX; fGP->fPmt.fY = pmtY; } } } } void CbmRichGeoManager::InitMirror() { TGeoIterator geoIterator(gGeoManager->GetTopNode()->GetVolume()); TGeoNode* curNode; geoIterator.Reset(); //mirror position\rotation TString mirrorName0("mirror_tile_type0"); TString mirrorName1("mirror_tile_type1"); TString mirrorName2("mirror_tile_type2"); TString mirrorName3("mirror_tile_type3"); //TString mirrorName4("mirror_tile_type4"); //TString mirrorName5("mirror_tile_type5"); // these names are needed for misaligned geometry TString mirrorMisAlignName0("mirror_tile_0"); TString mirrorMisAlignName1("mirror_tile_1"); TString mirrorMisAlignName2("mirror_tile_2"); TString mirrorMisAlignName3("mirror_tile_3"); geoIterator.Reset(); double maxTheta = 0.; double minTheta = 999999999.; double mirrorX = 0.; double mirrorY = 0.; double mirrorZ = 0.; double mirrorRadius = 0.; while ((curNode=geoIterator())) { TString nodeName(curNode->GetName()); TString nodePath; if (nodeName.Contains(mirrorName0) || nodeName.Contains(mirrorName1) || nodeName.Contains(mirrorName2) || nodeName.Contains(mirrorName3) || nodeName.Contains(mirrorMisAlignName0) || nodeName.Contains(mirrorMisAlignName1) || nodeName.Contains(mirrorMisAlignName2) || nodeName.Contains(mirrorMisAlignName3)) { geoIterator.GetPath(nodePath); const TGeoMatrix* curMatrix = geoIterator.GetCurrentMatrix(); const Double_t* curNodeTr = curMatrix->GetTranslation(); mirrorX = TMath::Abs(curNodeTr[0]); mirrorY = TMath::Abs(curNodeTr[1]); mirrorZ = TMath::Abs(curNodeTr[2]); const TGeoSphere* shape = dynamic_cast (curNode->GetVolume()->GetShape()); //const TGeoSphere* shape = (const TGeoSphere*)(curNode->GetVolume()->GetShape()); if (shape != NULL) { mirrorRadius = shape->GetRmin(); double theta1 = shape->GetTheta1() ; double theta2 = shape->GetTheta2() ; if (maxTheta < theta1 || maxTheta < theta2) { maxTheta = TMath::Max(theta1, theta2); } if (minTheta > theta1 || minTheta > theta2) { minTheta = TMath::Min(theta1, theta2); } } } } fGP->fMirrorTheta = -((maxTheta + minTheta)/2. - 90.) * TMath::DegToRad(); // rotation angle around x-axis fGP->fMirrorX = mirrorX; fGP->fMirrorY = mirrorY; fGP->fMirrorZ = mirrorZ; fGP->fMirrorR = mirrorRadius; } void CbmRichGeoManager::RotatePoint( TVector3 *inPos, TVector3 *outPos, Bool_t noTilting) { if (fGP == nullptr) { LOG(ERROR) << "CbmRichGeoManager::RotatePoint RICH geometry is not initialized. fGP == NULL" << FairLogger::endl; } if (fGP->fGeometryType == CbmRichGeometryTypeTwoWings) { RotatePointTwoWings(inPos, outPos, noTilting); } else if (fGP->fGeometryType == CbmRichGeometryTypeCylindrical) { RotatePointCyl(inPos, outPos, noTilting); } else { outPos->SetXYZ(inPos->X(), inPos->Y(), inPos->Z()); } } void CbmRichGeoManager::RotatePointTwoWings( TVector3 *inPos, TVector3 *outPos, Bool_t noTilting) { if (noTilting == false){ RotatePointImpl(inPos, outPos, fGP->fPmt.fPhi, fGP->fPmt.fTheta, fGP->fPmt.fX, fGP->fPmt.fY, fGP->fPmt.fZ); } else { outPos->SetXYZ(inPos->X(), inPos->Y(), inPos->Z()); } } void CbmRichGeoManager::RotatePointCyl( TVector3 *inPos, TVector3 *outPos, Bool_t noTilting, Bool_t noShift) { if (noTilting == false){ TGeoNode* node = gGeoManager->FindNode(inPos->X(), inPos->Y(), inPos->Z()); string path(gGeoManager->GetPath()); CbmRichRecGeoParPmt pmtPar = fGP->GetGeoRecPmtByBlockPath(path); RotatePointImpl(inPos, outPos, -TMath::Abs(pmtPar.fPhi), TMath::Abs(pmtPar.fTheta), TMath::Abs(pmtPar.fX), TMath::Abs(pmtPar.fY), TMath::Abs(pmtPar.fZ)); // After rotation wee need to shift point if (!noShift) { // All pmt blocks centers will be move to this Y position // TODO: We can also take smallest Y from all pmt blocks Double_t baseLineY = (outPos->Y() >= 0)?160. : -160.; //cm Double_t dY = pmtPar.fY - baseLineY; outPos->SetY(outPos->Y() - dY); // Calculate pmt block center after rotation TVector3 inPosPmt(pmtPar.fX, pmtPar.fY, pmtPar.fZ); TVector3 outPosPmt; RotatePointImpl(&inPosPmt, &outPosPmt, -TMath::Abs(pmtPar.fPhi), TMath::Abs(pmtPar.fTheta), TMath::Abs(pmtPar.fX), TMath::Abs(pmtPar.fY), TMath::Abs(pmtPar.fZ)); //RotatePointCyl(&inPosPmt, &outPosPmt, false, true); // calculate ideal position assuming the same width for all pmt blocks //TODO:Actually we need to implement general solution if pmt-block widths are not the same Double_t gap = fGP->fPmtStripGap; Double_t padding = gap / 2.; Double_t idealX = padding + 0.5 * pmtPar.fWidth + pmtPar.fPmtPositionIndexX * (pmtPar.fWidth + gap); if (outPos->X() < 0) idealX = -idealX; Double_t dX = idealX - outPosPmt.X(); outPos->SetX(outPos->X() + dX); } } else { outPos->SetXYZ(inPos->X(), inPos->Y(), inPos->Z()); } } void CbmRichGeoManager::RotatePointImpl( TVector3 *inPos, TVector3 *outPos, Double_t phi, Double_t theta, Double_t pmtX, Double_t pmtY, Double_t pmtZ) { Double_t xDet = 0.,yDet = 0.,zDet = 0.; Double_t x = inPos->X(); Double_t y = inPos->Y(); Double_t z = inPos->Z(); Double_t sinTheta = TMath::Sin(theta); Double_t cosTheta = TMath::Cos(theta); Double_t sinPhi = TMath::Sin(phi); Double_t cosPhi = TMath::Cos(phi); if (x > 0 && y > 0) { y -= pmtY; x -= pmtX; z -= pmtZ; //xDet = x*cosPhi + z*sinPhi;// - detZOrig*sinPhi; //yDet = -x*sinTheta*sinPhi + y*cosTheta + z*sinTheta*cosPhi; //zDet = -x*cosTheta*sinPhi - y*sinTheta + z*cosTheta*cosPhi; xDet = x * cosPhi - y * sinPhi * sinTheta + z * cosTheta * sinPhi; yDet = y * cosTheta + z * sinTheta; zDet = - x * sinPhi - y * sinTheta * cosPhi + z * cosTheta * cosPhi; yDet += pmtY; xDet += pmtX; zDet += pmtZ; } else if (x > 0 && y < 0) { y += pmtY; x -= pmtX; z -= pmtZ; // xDet = x*cosPhi + z*sinPhi;// - detZOrig*sinPhi; //yDet = x*sinTheta*sinPhi + y*cosTheta - z*sinTheta*cosPhi; //zDet = -x*cosTheta*sinPhi + y*sinTheta + z*cosTheta*cosPhi; xDet = x * cosPhi + y * sinPhi * sinTheta + z * cosTheta * sinPhi; yDet = y * cosTheta - z * sinTheta; zDet = - x * sinPhi + y * sinTheta * cosPhi + z * cosTheta * cosPhi; yDet -= pmtY; xDet += pmtX; zDet += pmtZ; } else if (x < 0 && y < 0) { y += pmtY; x += pmtX; z -= pmtZ; // xDet = x*cosPhi - z*sinPhi;// + detZOrig*sinPhi; //yDet = -x*sinTheta*sinPhi + y*cosTheta - z*sinTheta*cosPhi; //zDet = x*cosTheta*sinPhi + y*sinTheta + z*cosTheta*cosPhi; xDet = x * cosPhi - y * sinPhi * sinTheta - z * cosTheta * sinPhi; yDet = y * cosTheta - z * sinTheta; zDet = x * sinPhi + y * sinTheta * cosPhi + z * cosTheta * cosPhi; yDet -= pmtY; xDet -= pmtX; zDet += pmtZ; } else if (x < 0 && y > 0) { y -= pmtY; x += pmtX; z -= pmtZ; //xDet = x*cosPhi - z*sinPhi;// + detZOrig*sinPhi; //yDet = x*sinTheta*sinPhi + y*cosTheta + z*sinTheta*cosPhi; //zDet = x*cosTheta*sinPhi - y*sinTheta + z*cosTheta*cosPhi; xDet = x * cosPhi + y * sinPhi * sinTheta - z * cosTheta * sinPhi; yDet = y * cosTheta + z * sinTheta; zDet = x * sinPhi - y * sinTheta * cosPhi + z * cosTheta * cosPhi; yDet += pmtY; xDet -= pmtX; zDet += pmtZ; } else { outPos->SetXYZ(x,y,z); } outPos->SetXYZ(xDet,yDet,zDet); } Bool_t CbmRichGeoManager::IsPointInsidePmt( const TVector3* rotatedPoint) { if (fGP->fGeometryType == CbmRichGeometryTypeTwoWings) { CbmRichRecGeoPar* gp = CbmRichGeoManager::GetInstance().fGP; Double_t pmtPlaneX = gp->fPmt.fPlaneX; Double_t pmtPlaneY = gp->fPmt.fPlaneY; Double_t pmtWidth = gp->fPmt.fWidth; Double_t pmtHeight = gp->fPmt.fHeight; Double_t marginX = 2.; // [cm] Double_t marginY = 2.; // [cm] // upper pmt planes Double_t pmtYTop = TMath::Abs(pmtPlaneY) + pmtHeight + marginY; Double_t pmtYBottom = TMath::Abs(pmtPlaneY) - pmtHeight - marginY; Double_t absYDet = TMath::Abs(rotatedPoint->y()); Bool_t isYOk = (absYDet <= pmtYTop && absYDet >= pmtYBottom); Double_t pmtXMin = -TMath::Abs(pmtPlaneX) - pmtWidth - marginX; Double_t pmtXMax = TMath::Abs(pmtPlaneX) + pmtWidth + marginX; Bool_t isXOk = (rotatedPoint->x() >= pmtXMin && rotatedPoint->x() <= pmtXMax); return (isXOk && isYOk); } else if (fGP->fGeometryType == CbmRichGeometryTypeCylindrical) { return true; } else { return false; } }