// from ROOT #include "TVector3.h" // from hydra #include "hades.h" #include "hmdclayergeompar.h" #include "hmdcsizescells.h" #include "hruntimedb.h" #include "hkaldef.h" #include "hkaldetcradle.h" #include "hkalgeomtools.h" #include "hkalmdcmeaslayer.h" #include "hmaterials.h" #include "hkalmixture.h" ClassImp(HKalDetCradle) //_HADES_CLASS_DESCRIPTION /////////////////////////////////////////////////////////////////////////////// // // A detector cradle is a container that stores the detector's measurement layers. // The constructor HKalDetCradle(Int_t nLayInMdc) will automatically add // Mdc layers to the cradle. // //----------------------------------------------------------------------- /////////////////////////////////////////////////////////////////////////////// // ----------------------------------- // Ctors and Dtor // ----------------------------------- HKalDetCradle::HKalDetCradle() { matsMdc = NULL; matAir = NULL; nLayersInMdc = 1; } HKalDetCradle::HKalDetCradle(Int_t nLayInMdc, const TObjArray *customMats) : TObjArray(kNumSecs*kNumMdcs*nLayInMdc) { // Default contructor that adds Mdc layers to the cradle. // The Mdc layers will be added to the cradle in the following order: // 1. Sort by sector number (0..5) // 2. Sort by module number (0..3) in a sector // 3. Sort by layer number of module in sector // nLayInMdc: the number of mdc layers that should be added to each sector in a module. // customMats: Materials for MDCs and materials between MDCs. If empty // standard materials will be used matsMdc = new HKalMixture*[kNumMdcs]; matAir = NULL; nLayersInMdc = nLayInMdc; if(customMats) { if(customMats->GetEntries() == kNumMdcs+1) { for(Int_t iMod = 0; iMod < kNumMdcs; iMod++) { HKalMixture *mat = dynamic_cast(customMats->At(iMod)); if(mat) { matsMdc[iMod] = mat; } else { Error("HKalDetCradle()", Form("No material found for MDC %i", iMod)); } } HKalMixture *mat = dynamic_cast(customMats->At(kNumMdcs)); if(mat) { matAir = mat; } else { Error("HKalDetCradle()", Form("No object found at index %i", kNumMdcs)); } } else { Error("HKalDetCradle()", Form("Size of materials array is %i. Expected %i.", customMats->GetEntries(), kNumMdcs)); } } else { initMdcMaterials(); } for(Int_t iSec = 0; iSec < kNumSecs; iSec++) { for(Int_t iMod = 0; iMod < kNumMdcs; iMod++) { #if kalDebug > 0 if(!HMdcSizesCells::getObject()) { Warning("HKalDetCradle(Int_t)", "HMdcSizesCells has not been initialized. Measurement layers not added to detector cradle."); break; } #endif addMdcLayersToCradle(HMdcSizesCells::getObject(), iSec, iMod, nLayersInMdc); } } } HKalDetCradle::~HKalDetCradle() { delete [] matsMdc; delete matAir; } // ----------------------------------- // Implementation of public methods // ----------------------------------- void HKalDetCradle::addMdcLayersToCradle(const HMdcSizesCells *fSizesCells, Int_t sector, Int_t module, Int_t nLayers) { // Create a mdc detector with nLayers measurement layers. Each time the constructor // is called, a new set if KalMdcMeasLayer objects will be added to the TObjArray. // Module and sector identify the drift chamber. // // Input: // fSizesCells: object that can do coordinate transformation layer<->sector // sector: the sector where this detector is located (needed for coordinate transformation) // module: the module where this detector is located (needed for coordinate transformation) // nLayers: number of layers in detector HRuntimeDb *rtdb = gHades->getRuntimeDb(); HMdcLayerGeomPar *layerGeomPar = (HMdcLayerGeomPar*)rtdb->getContainer("MdcLayerGeomPar"); for (Int_t layer = 0; layer < nLayers; layer++) { // A measurement vector is modeled as a plane. A plane is defined by // a point on the plane and a normal vector. // The Center of layer is the Null vector and the normal vector is (0,0,1) // in the layer coordinate system. // These vectors need to be transformed into sector coordinates. #if kalDebug > 0 if(!secOk(sector)) { Error("addMdcLayersToCradle()", Form("Invalid sector number: %i. Layers not added to cradle.", sector)); break; } if(!modOk(module)) { Error("addMdcLayersToCradle()", Form("Invalid module number: %i. Layers not added to cradle.", module)); break; } if(!fSizesCells) { Error("addMdcLayersToCradle()", "HMdcSizesCells is not initialized."); } #endif Double_t x = 0.; Double_t y = 0.; Double_t z = 0.; if(nLayers == 1) { HKalGeomTools::TransModuleToSector(x, y, z, sector, module); //((*fSizesCells)[sector][module]).transFrom(x,y,z); // transform to sector coordinate sys } else { HKalGeomTools::TransLayerToSector(x, y, z, sector, module, layer); //((*fSizesCells)[sector][module][layer]).transFrom(x,y,z); // transform to sector coordinate sys } TVector3 vPlaneCenter(x, y, z); x = 0.; y = 0.; z = 1.; if(nLayers == 1) { HKalGeomTools::TransModuleToSector(x, y, z, sector, module); //((*fSizesCells)[sector][module]).transFrom(x,y,z); // transform to sector coordinate system } else { HKalGeomTools::TransLayerToSector(x, y, z, sector, module, layer); //((*fSizesCells)[sector][module][layer]).transFrom(x,y,z); // transform to sector coordinate sys } TVector3 vZ(x, y, z); TVector3 vPlaneNormal = vZ - vPlaneCenter; vPlaneNormal.Unit(); // Choose the normal vector, so that the scalar product of the normal vector // and the point on the plane is positive. if(vPlaneCenter * vPlaneNormal < 0) { vPlaneNormal *= -1.; } HKalMdcMeasLayer *meas; HMdcLayerGeomParLay &layerGeomParLay = (HMdcLayerGeomParLay&)(*layerGeomPar)[sector][module][layer]; Double_t layerThickness = layerGeomParLay.getCatDist(); if(nLayers == 1) { layerThickness *= 8.; // 6 layers + 2 * layer thickness outside } if(nLayers != 1 && layer == 0) { meas = new HKalMdcMeasLayer(vPlaneCenter, vPlaneNormal, matsMdc[module], matAir, module, sector, layer, nLayersInMdc, layerThickness); } else { meas = new HKalMdcMeasLayer(vPlaneCenter, vPlaneNormal, matsMdc[module], matAir, module, sector, layer, nLayersInMdc, layerThickness); } install(meas); } SetOwner(); } const HKalMdcMeasLayer* HKalDetCradle::getMdcLayerAt(Int_t sector, Int_t module, Int_t layer) const { // Returns the mdc layer in a sector // Input: // sector: the sector where this detector is located (needed for coordinate transformation) // module: the module where this detector is located (needed for coordinate transformation) // layer: layer number in the sector of the mdc module. Int_t layIdx = (sector*kNumMdcs + module)*nLayersInMdc; if(layer != 6) { // layer 6 = midplane layIdx += layer; } #ifdef kalDebug if(layIdx < 0 || layIdx > GetEntries() - 1) { Error("getMdcLayerAt()", Form("Mdc layer at sector %i, module %i, layer %i not found in detector cradle. Tried to access element at index %i.", sector, module, layer, layIdx)); exit(1); } if(!At(layIdx)->InheritsFrom("HKalMdcMeasLayer")) { Error("getMdcLayerAt()", Form("Object at index %i in detector cradle is of class %s. Expected class HKalMdcMeasLayer.", layIdx, At(layIdx)->ClassName())); exit(1); } return (HKalMdcMeasLayer*)At(layIdx); #else return (HKalMdcMeasLayer*)At(layIdx); #endif } HKalMixture* HKalDetCradle::getMdcMaterial(mdcComponent comp, Int_t module) const { // Returns the material of an Mdc layer component. // Input: // comp: The Mdc component (Cathode/field/sense wires, gas, window, mylar foil) // module: The module the Mdc layer is located in. if(!modOk(module)) { Error("getMdcMaterial()", Form("%i is not a valid module number", module)); return NULL; } HKalMixture *mdcMix = NULL; Float_t mat[kNumMdcMats]; switch(comp) { // Cathode and field wires consists of the same materials. case kCatWi: case kFieldWi: mdcMix = new HKalMixture("Aluminum", "altu", 1); getProperties(mat, "al"); mdcMix->defineElement(0, mat, 1.F); mdcMix->calcMixture(); break; case kSenseWi: // Material for sense wires. Float_t vol[2], rho[2], w[2]; rho[0] = HMaterials::getDensity("Copper"); rho[1] = HMaterials::getDensity("Beryllium"); // Mdc layers contain a 84:16 volume percent mixture of Argon-Isobutane. vol[0] = .98; vol[1] = .02; // Get mass fractions for gas mixtue. mdcMix = new HKalMixture("Copper/Beryllium 98/2", "co/be", 2); mdcMix->calcMassFracFromVolFrac(w, 2, rho, vol); getProperties(mat, "Copper"); mdcMix->defineElement(0, mat, w[0]); getProperties(mat, "Beryllium"); mdcMix->defineElement(1, mat, w[1]); mdcMix->calcMixture(); break; case kWindow: // Material for mylar foil. mdcMix = new HKalMixture("Mylar", "mylar", 1); getProperties(mat, "mylar"); mdcMix->defineElement(0, mat, 1.F); mdcMix->calcMixture(); break; case kGas: // Gas in Mdc 1 is Argon-Isobutane if(module == 0) { Float_t vol[2], rho[2], w[2]; rho[0] = HMaterials::getDensity("argon"); rho[1] = HMaterials::getDensity("CO2"); // Mdc layers contain a 84:16 volume percent mixture of Argon-Isobutane. vol[0] = .70; vol[1] = .30; // Get mass fractions for gas mixtue. mdcMix = new HKalMixture("Argon-Isobutane 70:30", "arco2_70_30", 2); mdcMix->calcMassFracFromVolFrac(w, 2, rho, vol); getProperties(mat, "argon"); mdcMix->defineElement(0, mat, w[0]); getProperties(mat, "CO2"); mdcMix->defineElement(1, mat, w[1]); mdcMix->calcMixture(); // Gas in the other Mdcs is Helium-Isobutane. } else { Float_t vol[2], rho[2], w[2]; rho[0] = HMaterials::getDensity("helium"); rho[1] = HMaterials::getDensity("butane"); // Mdc layers contain a 60:40 volume percent mixture of Helium-Isobutane. vol[0] = .84; vol[1] = .16; // Get mass fractions. mdcMix = new HKalMixture("Helium-Isobutane 60:40", "hebutane", 2); mdcMix->calcMassFracFromVolFrac(w, 2, rho, vol); getProperties(mat, "helium"); mdcMix->defineElement(0, mat, w[0]); getProperties(mat, "butane"); mdcMix->defineElement(1, mat, w[1]); mdcMix->calcMixture(); } break; } return mdcMix; } void HKalDetCradle::getMdcMaterialVolumeFracs(Float_t volFracs[kNumMdcMats], Int_t sector, Int_t module, Int_t layer) const { // Calculate the volume fractions of the different mdc components. // volFracs: the volume fractions of the components // module: the module number of the mdc layer // sector: the sector number of the mdc layer // layer: the layer number if(!layOk(layer)) { Error("getMdcMaterialVolumeFracs()", Form("%i is not a valid layer number", layer)); return; } if(!modOk(module)) { Error("getMdcMaterialVolumeFracs()", Form("%i is not a valid module number", module)); return; } HRuntimeDb *rtdb = gHades->getRuntimeDb(); HMdcLayerGeomPar *layerGeomPar = (HMdcLayerGeomPar*)rtdb->getContainer("MdcLayerGeomPar"); HMdcLayerGeomParLay &layerGeomParLay = (HMdcLayerGeomParLay&)(*layerGeomPar)[sector][module][layer]; Float_t nLayers = 6; // number of layers in mdc plane // Thickness of cell in mm = distance between cathode layers. Float_t xCell = layerGeomParLay.getCatDist(); Float_t xMylar = .012; // thickness of mylar foils in mm //Float_t xPlane = xCell * (nLayers+2) + 2 * xMylar; // thickness of mdc plane in mm // Distance (pitch) between two wires in mm. Float_t pitchCatWi; if(module == 3) { pitchCatWi = 4.F; } else { pitchCatWi = 2.F; } Float_t pitchFieldWi = xCell; Float_t pitchSenseWi = xCell; Double_t pi = TMath::Pi(); // Radii of wires in mm. Float_t rCatWi = layerGeomParLay.getCathodeWireThickness() / 2; Float_t rFieldWi; if(module == 0 || module == 1) { rFieldWi = 0.04; } else { rFieldWi = 0.05; } Float_t rSenseWi; if(module == 3) { rSenseWi = 0.015; } else { rSenseWi = 0.01; } Float_t thicknessCell[kNumMdcMats]; // Spread out wires with circular cross sectional area over the wires' pitch // to get a wire with a rectangular cross sectional area of equal proportions. // Area of circle: pi * r^2 // Area of rectangle: pitch * x // => Spread out thickness of wire: x = pi * r^2 / pitch thicknessCell[kCatWi] = pi * rCatWi * rCatWi / pitchCatWi; thicknessCell[kFieldWi] = pi * rFieldWi * rFieldWi / pitchFieldWi; thicknessCell[kSenseWi] = pi * rSenseWi * rSenseWi / pitchSenseWi; thicknessCell[kWindow] = xMylar; thicknessCell[kGas] = xCell - (thicknessCell[kCatWi] + thicknessCell[kFieldWi] + thicknessCell[kSenseWi]); // Materials' thickness added over the whole Mdc plane. Float_t thickness[kNumMdcMats]; thickness[kCatWi] = (nLayers+1) * thicknessCell[kCatWi]; thickness[kFieldWi] = nLayers * thicknessCell[kFieldWi]; thickness[kSenseWi] = nLayers * thicknessCell[kSenseWi]; thickness[kWindow] = 2 * thicknessCell[kWindow]; thickness[kGas] = nLayers * thicknessCell[kGas] + xCell * 2; Float_t thicknessTotal = 0.F; for(Int_t i = 0; i < kNumMdcMats; i++) { thicknessTotal += thickness[i]; } for(Int_t i = 0; i < kNumMdcMats; i++) { volFracs[i] = thickness[i] / thicknessTotal; } } void HKalDetCradle::initMdcMaterials() { Float_t mat[5]; // Storage for material properties. matAir = new HKalMixture("Air", "air", 1); getProperties(mat, "air"); matAir->defineElement(0, mat, 1.F); matAir->calcMixture(); // Materials of mdc planes. // Materials of the different components (gas, wires, window foils) of an mdc plane. HKalMixture const *matsMdcComp [kNumMdcMats]; // Density of mdc components. Float_t rhoMatsMdcComp [kNumMdcMats]; // Fractions on total volume of the different mdc components. Float_t volFracsMdcComp [kNumMdcMats]; // Mass fractions of the different mdc components. Float_t massFracsMdcComp[kNumMdcMats]; Int_t sector = 0; Int_t layer = 0; for(Int_t iModule = 0; iModule < kNumMdcs; iModule++) { // Get the materials of the mdc components. for(Int_t iComp = 0; iComp < kNumMdcMats; iComp++) { matsMdcComp [iComp] = getMdcMaterial((mdcComponent)iComp, iModule); rhoMatsMdcComp[iComp] = matsMdcComp[iComp]->GetDensity(); } // Calculate volume fractions of mdc components. getMdcMaterialVolumeFracs(volFracsMdcComp, sector, iModule, layer); // Calculate mass fractions from volume fractions and densities. HKalMixture::calcMassFracFromVolFrac(massFracsMdcComp, kNumMdcMats, rhoMatsMdcComp, volFracsMdcComp); // Make a mixture of the mdc components' materials. matsMdc[iModule] = new HKalMixture(Form("Materials for mdc plane %i.", iModule), "mdc", kNumMdcMats); for(Int_t iComp = 0; iComp < kNumMdcMats; iComp++) { matsMdc[iModule]->defineElement(iComp, matsMdcComp[iComp]->GetA(), matsMdcComp[iComp]->GetZ(), matsMdcComp[iComp]->GetDensity(), matsMdcComp[iComp]->GetExciteEner(), matsMdcComp[iComp]->GetRadLength(), massFracsMdcComp[iComp]); } matsMdc[iModule]->calcMixture(); for(Int_t iComp = 0; iComp < kNumMdcMats; iComp++) { delete matsMdcComp[iComp]; } } } void HKalDetCradle::getProperties(Float_t mat[5], const Char_t *id) { // Returns the atomic mass, atomic number, density, excitation energy and // radiation length in an array and converts the units from cm to mm and eV to MeV. // The array index where a material property is stored is determined by the matIdx // enum in hkaldef.h. // Output: // mat: array with material properties. // Input: // id: material name as in HMaterial mat[Kalman::kMatIdxA] = HMaterials::getA(id); mat[Kalman::kMatIdxZ] = HMaterials::getZ(id); mat[Kalman::kMatIdxDensity] = HMaterials::getDensity(id); mat[Kalman::kMatIdxExEner] = HMaterials::getExciteEner(id); mat[Kalman::kMatIdxRadLength] = HMaterials::getRadLength(id); Float_t corrRho = 1.e-3; // g/cm^3 -> g/mm^3 Float_t corrI = 1.e-6; // eV -> MeV Float_t corrX0 = 10.F; // cm -> mm mat[Kalman::kMatIdxDensity] *= corrRho; mat[Kalman::kMatIdxExEner] *= corrI; mat[Kalman::kMatIdxRadLength] *= corrX0; } void HKalDetCradle::install(HKalMdcMeasLayer *layer) { // Adds a new mdc measurement to the array. Add(layer); }