diff --git a/core/detectors/trd/CbmTrdRadiator.cxx b/core/detectors/trd/CbmTrdRadiator.cxx index 0adeb44987e41454f6ac59a08e877f9b705d5c7a..43e2a654087a9ddc74f7a2503b534a6cc55f0583 100644 --- a/core/detectors/trd/CbmTrdRadiator.cxx +++ b/core/detectors/trd/CbmTrdRadiator.cxx @@ -18,34 +18,18 @@ #include <TRandom3.h> // for gRandom #include <TVector3.h> // for TVector3 -#include <iomanip> // for setprecision, __iom_t5 +#include <iomanip> // for setprecision, __iom_t5 + #include <stdio.h> // for printf, sprintf #include <stdlib.h> // for exit using std::setprecision; - -// ----- Default constructor --------------------------------------- -CbmTrdRadiator::CbmTrdRadiator() : CbmTrdRadiator(kTRUE, 130, 0.0013, 0.02) {} -//----------------------------------------------------------------------------- - -// ----- Constructor -------------------------------------------------- -CbmTrdRadiator::CbmTrdRadiator(Bool_t SimpleTR, - Int_t Nfoils, - Float_t FoilThick, - Float_t GapThick) - : CbmTrdRadiator(SimpleTR, Nfoils, FoilThick, GapThick, "polyethylen", "") {} -//----------------------------------------------------------------------------- - -// ----- Constructor -------------------------------------------------- -CbmTrdRadiator::CbmTrdRadiator(Bool_t SimpleTR, - Int_t Nfoils, - Float_t FoilThick, - Float_t GapThick, - TString material, - TString prototype) - : fWindowFoil("") - , fRadType(prototype) +// ----- Default constructor ------------------------------------------------ +CbmTrdRadiator::CbmTrdRadiator(Bool_t SimpleTR, Int_t Nfoils, Float_t FoilThick, Float_t GapThick, TString material, + TString window) // If "Kapton" is set, 0.05 mu aluminum will be added + : fWindowFoil(window) + , fRadType("") , fDetType(-1) , fFirstPass(kTRUE) , fSimpleTR(SimpleTR) @@ -58,13 +42,14 @@ CbmTrdRadiator::CbmTrdRadiator(Bool_t SimpleTR, , fGapDens(-1.) , fFoilOmega(-1.) , fGapOmega(-1.) - , fnPhotonCorr(1.0) + , fnPhotonCorr(0.65) , fFoilThickCorr(1.) , fGapThickCorr(1.) , fGasThickCorr(1.) , fnTRprod(-1) , fWinDens(-1.) , fWinThick(-1.) + , WINDOW() , fCom1(-1.) , fCom2(-1.) , fSpBinWidth((Float_t) fSpRange / (Float_t) fSpNBins) @@ -80,39 +65,39 @@ CbmTrdRadiator::CbmTrdRadiator(Bool_t SimpleTR, , fnTRabs() , fnTRab(-1.) , fELoss(-1.) - , fMom(-1., -1., -1.) { + , fMom(-1., -1., -1.) +{ for (Int_t i = 0; i < fNMom; i++) { fFinal[i] = nullptr; } - //Init(); - CreateHistograms(); } //----------------------------------------------------------------------------- -// ----- Constructor -------------------------------------------------- -CbmTrdRadiator::CbmTrdRadiator(Bool_t SimpleTR, TString prototype) - : fWindowFoil("Kapton") - , //which is the gas window of MS2012 prototypes. Anything else will result in a mylar gas window -> MuBu default - fRadType(prototype) +// ----- Constructor using predefined radiator prototype------------------------------------ +CbmTrdRadiator::CbmTrdRadiator(Bool_t SimpleTR, TString prototype, + TString window) // If "Kapton" is set, 0.05 mu aluminum will be added + : fWindowFoil(window) + , fRadType(prototype) , fDetType(-1) , fFirstPass(kTRUE) , fSimpleTR(SimpleTR) - , fNFoils(130) - , fFoilThick(0.0013) - , fGapThick(0.02) + , fNFoils(337) + , fFoilThick(0.0012) + , fGapThick(0.09) , fFoilMaterial("") , fGasThick(-1) , fFoilDens(-1.) , fGapDens(-1.) , fFoilOmega(-1.) , fGapOmega(-1.) - , fnPhotonCorr(1.0) + , fnPhotonCorr(0.65) , fFoilThickCorr(1.) , fGapThickCorr(1.) , fGasThickCorr(1.) , fnTRprod(-1) , fWinDens(-1.) , fWinThick(-1.) + , WINDOW() , fCom1(-1.) , fCom2(-1.) , fSpBinWidth((Float_t) fSpRange / (Float_t) fSpNBins) @@ -128,19 +113,17 @@ CbmTrdRadiator::CbmTrdRadiator(Bool_t SimpleTR, TString prototype) , fnTRabs() , fnTRab(-1.) , fELoss(-1.) - , fMom(-1., -1., -1.) { + , fMom(-1., -1., -1.) +{ for (Int_t i = 0; i < fNMom; i++) { fFinal[i] = nullptr; } - - //Init(SimpleTR, prototype); - //SetRadPrototype(prototype); - CreateHistograms(); } //----------------------------------------------------------------------------- // ----- Destructor ---------------------------------------------------- -CbmTrdRadiator::~CbmTrdRadiator() { +CbmTrdRadiator::~CbmTrdRadiator() +{ // FairRootManager *fManager = FairRootManager::Instance(); // fManager->Write(); if (fSpectrum) { @@ -157,7 +140,8 @@ CbmTrdRadiator::~CbmTrdRadiator() { } // ----- Create Histogramms -------------------------------------------------- -void CbmTrdRadiator::CreateHistograms() { +void CbmTrdRadiator::CreateHistograms() +{ // Create the needed histograms @@ -172,22 +156,13 @@ void CbmTrdRadiator::CreateHistograms() { fSpectrum = new TH1D("fSpectrum", "TR spectrum", fSpNBins, SpLower, SpUpper); if (fWinSpectrum) delete fWinSpectrum; - fWinSpectrum = new TH1D( - "fWinSpectrum", "TR spectrum in Mylar", fSpNBins, SpLower, SpUpper); + fWinSpectrum = new TH1D("fWinSpectrum", "TR spectrum in Mylar", fSpNBins, SpLower, SpUpper); if (fDetSpectrum) delete fDetSpectrum; - fDetSpectrum = new TH1D("fDetSpectrum", - "TR spectrum escaped from detector", - fSpNBins, - SpLower, - SpUpper); + fDetSpectrum = new TH1D("fDetSpectrum", "TR spectrum escaped from detector", fSpNBins, SpLower, SpUpper); if (fDetSpectrumA) delete fDetSpectrumA; - fDetSpectrumA = new TH1D("fDetSpectrumA", - "TR spectrum absorbed in detector", - fSpNBins, - SpLower, - SpUpper); + fDetSpectrumA = new TH1D("fDetSpectrumA", "TR spectrum absorbed in detector", fSpNBins, SpLower, SpUpper); for (Int_t i = 0; i < fNMom; i++) { sprintf(name, "fFinal%d", i + 1); @@ -197,166 +172,122 @@ void CbmTrdRadiator::CreateHistograms() { } } -// ----- Init function ---------------------------------------------------- -void CbmTrdRadiator::Init(Bool_t SimpleTR, - Int_t Nfoils, - Float_t FoilThick, - Float_t GapThick, - TString material = "polyethylen") { +// ----- Init function ---------------------------------------------------- +void CbmTrdRadiator::Init() +{ LOG(info) << "================CbmTrdRadiator==============="; - // Set initial parameters defining the radiator + CreateHistograms(); - fNFoils = Nfoils; - fGapThick = GapThick; - fFoilThick = FoilThick; - fSimpleTR = SimpleTR; - fnPhotonCorr = 1.0; + //Check if a radiator type has been set + if (fRadType == "") { + //Set fnPhotonCorr also if no prototype has been specified + fnPhotonCorr = 0.65; + //Check specified radiator types + } + else if (fRadType == "tdr18") { //Default TDR18 radiator 30cm Box + 1.5 cm in carbon grid + fFoilMaterial = "pefoam20"; + fNFoils = 337; + fGapThick = 900 * 0.0001; + fFoilThick = 12 * 0.0001; + fnPhotonCorr = 0.65; + } + else if (fRadType == "B++" || fRadType == "K++") { + fFoilMaterial = "pokalon"; + fNFoils = 350; + fGapThick = 0.07; + fFoilThick = 0.0024; + fnPhotonCorr = 0.65; + } + else if (fRadType == "G30") { + fFoilMaterial = "pefiber"; + fNFoils = 1791; + fGapThick = 50 * 0.0001; + fFoilThick = 17 * 0.0001; + fnPhotonCorr = 0.65; + } + else { // If radiator typ is neither empty nor one of the above (because of typo or else) set defaults + + LOG(warn) << "CbmTrdRadiator::Init : *********** Radiator prototype not known ***********"; + LOG(warn) << "CbmTrdRadiator::Init : switch to default parameters "; + fFoilMaterial = "pefoam20"; + fNFoils = 337; + fGapThick = 900 * 0.0001; + fFoilThick = 12 * 0.0001; + fnPhotonCorr = 0.65; + } - CreateHistograms(); + //Material has been specified now, either in the constructor or set by radiator type above. Now, go through materials and set properties accordingly - LOG(info) << "CbmTrdRadiator::Init : Foil material : " << material; - LOG(info) << "CbmTrdRadiator::Init : Nr. of foils : " << Nfoils; - LOG(info) << "CbmTrdRadiator::Init : Foil thickness : " - << setprecision(4) << FoilThick << " cm"; - LOG(info) << "CbmTrdRadiator::Init : Gap thickness : " << GapThick - << " cm"; - LOG(info) << "CbmTrdRadiator::Init : Simple TR production: " - << setprecision(2) << SimpleTR; - - // material dependent parameters for the radiator material - // (polyethylen) - if (material == "" || material == "polyethylen") { - fFoilDens = 0.92; // [g/cm3 ] - fFoilOmega = 20.9; //plasma frequency ( from Anton ) - } else if (material == "pefoam20") { - // (polyethylen) + + // material dependent parameters for the gas between the foils of the radiator + //Properties of the air gaps between foils + fGapDens = 0.001205; // [g/cm3] + fGapOmega = 0.7; //plasma frequency ( from Anton ) + //material = fFoilMaterial; + if (fFoilMaterial == "pefoam20") { + // (polyethylen foam foil) fFoilDens = 0.90; // [g/cm3 ] fFoilOmega = 20.64; //plasma frequency ( from Cyrano ) - } else if (material == "pefiber") { + } + else if (fFoilMaterial == "polyethylen") { + // (polyethylen) + fFoilDens = 0.92; // [g/cm3 ] + fFoilOmega = 20.9; //plasma frequency ( from Anton ) + } + else if (fFoilMaterial == "pefiber") { // (polyethylen) fFoilDens = 0.90; // [g/cm3 ] fFoilOmega = 20.64; //plasma frequency ( from Cyrano ) - } else if (material == "mylar") { + } + else if (fFoilMaterial == "mylar") { // (Mylar) fFoilDens = 1.393; // [g/cm3 ] fFoilOmega = 24.53; //plasma frequency ( from Cyrano ) - } else if (material == "pokalon") { - // (Pokalon) + } + else if (fFoilMaterial == "pokalon") { + // (PokalonN470) fFoilDens = 1.150; // [g/cm3 ] - fFoilOmega = 22.50; //plasma frequency ( from Cyrano ) - } else { - // (polyethylen) - material = "polyethylen"; - fFoilDens = 0.92; // [g/cm3 ] - fFoilOmega = 20.9; //plasma frequency ( from Anton ) + fFoilOmega = 22.43; //plasma frequency ( from Cyrano ) + } + else { // If material is empty or one of the above (because of typo or else) set default + LOG(warn) << "CbmTrdRadiator::Init : *********** Radiator material not known ***********"; + LOG(warn) << "CbmTrdRadiator::Init : switch to default material "; + // (polyethylen foam foil) + fFoilMaterial = "pefoam20"; + fFoilDens = 0.90; // [g/cm3 ] + fFoilOmega = 20.64; //plasma frequency ( from Cyrano ) } - fFoilMaterial = material; - // material dependent parameters for the gas between the foils of the - // radiator - // changed gap density at 16.04.07 FU - // density of air is 0.001205, 0.00186 - //fGapDens = 0.00186; // [g/cm3] - fGapDens = 0.001205; // [g/cm3] - fGapOmega = 0.7; //plasma frequency ( from Anton ) // foil between radiator and gas chamber - // TODO: implement kapton foil also in trd geometry - fWinDens = 1.39; // [g/cm3] - fWinThick = 0.0025; // [cm] - - - // Get all the gas properties from CbmTrdGas - CbmTrdGas* fTrdGas = CbmTrdGas::Instance(); - if (fTrdGas == 0) { - fTrdGas = new CbmTrdGas(); - fTrdGas->Init(); + if (fWindowFoil == "Kapton") { // 0.05 mu Aluminum is added in SigmaWin() for Kapton + fWinDens = 1.42; //[g/cm3] ?? define values ?? default mylar + fWinThick = 0.0025; // [cm] ?? define values ?? default mylar } - - fCom2 = fTrdGas->GetNobleGas(); - fCom1 = fTrdGas->GetCO2(); - fDetType = fTrdGas->GetDetType(); - fGasThick = fTrdGas->GetGasThick(); - - LOG(info) << "CbmTrdRadiator::Init : Detector noble gas : " - << fTrdGas->GetNobleGas(); - LOG(info) << "CbmTrdRadiator::Init : Detector cruncher gas: " - << fTrdGas->GetCO2(); - LOG(info) << "CbmTrdRadiator::Init : Detector type : " - << fTrdGas->GetDetType(); - LOG(info) << "CbmTrdRadiator::Init : Detector gas thick. : " - << fTrdGas->GetGasThick() << " cm"; - - //LOG(info) << "************* End of Trd Radiator Init **************"; - - // If simplified version is used one didn't calculate the TR - // for each CbmTrdPoint in full glory, but create at startup - // some histograms for several momenta which are later used to - // get the TR much faster. - if (fSimpleTR == kTRUE) { - - ProduceSpectra(); - - TFile* oldfile = gFile; - TFile* f1 = new TFile("TRhistos.root", "recreate"); - f1->cd(); - - for (Int_t i = 0; i < fNMom; i++) { - fFinal[i]->Write(); - } - - f1->Close(); - f1->Delete(); - gFile = oldfile; + else if (fWindowFoil != "Mylar") { + if (!WINDOW.Init(fWindowFoil)) fWindowFoil = "Mylar"; + } + if (fWindowFoil == "Mylar") { + fWinDens = 1.39; //[g/cm3] + fWinThick = 0.0100; //[cm] + } + if (!WINDOW.fN) { + LOG(info) << "CbmTrdRadiator::Init : Window type : " << fWindowFoil; + LOG(info) << "CbmTrdRadiator::Init : Window density : " << fWinDens << " g/cm^3"; + LOG(info) << "CbmTrdRadiator::Init : Window thickness : " << fWinThick << " cm"; } -} -//---------------------------------------------------------------------------- - -// ----- Init function ---------------------------------------------------- -void CbmTrdRadiator::Init(Bool_t SimpleTR, - Int_t Nfoils, - Float_t FoilThick, - Float_t GapThick) { - - LOG(info) << "================CbmTrdRadiator==============="; - // Set initial parameters defining the radiator - fNFoils = Nfoils; - fGapThick = GapThick; - fFoilThick = FoilThick; - fSimpleTR = SimpleTR; - fnPhotonCorr = 1.0; - CreateHistograms(); - - //LOG(info) << "CbmTrdRadiator::Init : Foil material : " << material; - LOG(info) << "CbmTrdRadiator::Init : Nr. of foils : " << Nfoils; - LOG(info) << "CbmTrdRadiator::Init : Foil thickness : " - << setprecision(4) << FoilThick << " cm"; - LOG(info) << "CbmTrdRadiator::Init : Gap thickness : " << GapThick - << " cm"; - LOG(info) << "CbmTrdRadiator::Init : Simple TR production: " - << setprecision(2) << SimpleTR; - - - // material dependend parameters for the radiator material - // (polyethylen) - fFoilDens = 0.92; // [g/cm3 ] - fFoilOmega = 20.9; //plasma frequency ( from Anton ) - - // material dependent parameters for the gas between the foils of the - // radiator - // changed gap density at 16.04.07 FU - // density of air is 0.001205, 0.00186 - //fGapDens = 0.00186; // [g/cm3] - fGapDens = 0.001205; // [g/cm3] - fGapOmega = 0.7; //plasma frequency ( from Anton ) - - // foil between radiator and gas chamber - // TODO: implement kapton foil also in trd geometry - fWinDens = 1.39; // [g/cm3] - fWinThick = 0.0025; // [cm] + LOG(info) << "CbmTrdRadiator::Init : Prototype : " << fRadType; + LOG(info) << "CbmTrdRadiator::Init : Scaling factor : " << fnPhotonCorr; + LOG(info) << "CbmTrdRadiator::Init : Foil material : " << fFoilMaterial; + LOG(info) << "CbmTrdRadiator::Init : Nr. of foils : " << fNFoils; + LOG(info) << "CbmTrdRadiator::Init : Foil thickness : " << setprecision(4) << fFoilThick << " cm"; + LOG(info) << "CbmTrdRadiator::Init : Gap thickness : " << fGapThick << " cm"; + LOG(info) << "CbmTrdRadiator::Init : Simple TR production: " << setprecision(2) << fSimpleTR; + LOG(info) << "CbmTrdRadiator::Init : Foil plasm. frequ. : " << fFoilOmega << " keV"; + LOG(info) << "CbmTrdRadiator::Init : Gap plasm. frequ. : " << fGapOmega << " keV"; // Get all the gas properties from CbmTrdGas @@ -370,14 +301,13 @@ void CbmTrdRadiator::Init(Bool_t SimpleTR, fCom1 = fTrdGas->GetCO2(); fDetType = fTrdGas->GetDetType(); fGasThick = fTrdGas->GetGasThick(); - /* - LOG(info) << "Detector noble gas : "<< fTrdGas->GetNobleGas(); - LOG(info) << "Detector cruncher gas: "<< fTrdGas->GetCO2(); - LOG(info) << "Detector type : "<< fTrdGas->GetDetType(); - LOG(info) << "Detector gas thick. : "<< fTrdGas->GetGasThick() <<" cm"; - */ - //LOG(info) << "************* End of Trd Radiator Init **************"; + LOG(info) << "CbmTrdRadiator::Init : Detector noble gas : " << fTrdGas->GetNobleGas(); + LOG(info) << "CbmTrdRadiator::Init : Detector quenching gas: " << fTrdGas->GetCO2(); + LOG(info) << "CbmTrdRadiator::Init : Detector type : " << fTrdGas->GetDetType(); + LOG(info) << "CbmTrdRadiator::Init : Detector gas thick. : " << fTrdGas->GetGasThick() << " cm"; + + //LOG(info) << "************* End of Trd Radiator Init **************"; // If simplified version is used one didn't calculate the TR // for each CbmTrdPoint in full glory, but create at startup // some histograms for several momenta which are later used to @@ -386,236 +316,111 @@ void CbmTrdRadiator::Init(Bool_t SimpleTR, ProduceSpectra(); - TFile* oldfile = gFile; - TFile* f1 = new TFile("TRhistos.root", "recreate"); - f1->cd(); + TFile* f1 = new TFile("TRhistos.root", "recreate"); for (Int_t i = 0; i < fNMom; i++) { fFinal[i]->Write(); } - f1->Close(); f1->Delete(); - gFile = oldfile; } } -//---------------------------------------------------------------------------- -// ----- Init function ---------------------------------------------------- -void CbmTrdRadiator::Init() { - LOG(info) << "================CbmTrdRadiator==============="; - TString material; - CreateHistograms(); - if (fRadType == "") { - fFoilMaterial = "polyethylen"; - fnPhotonCorr = 1.0; - - /* - LOG(info) << "CbmTrdRadiator::Init : Foil material : " << fFoilMaterial; - //LOG(info) << "CbmTrdRadiator::Init : Nr. of foils : " << fNfoils; - LOG(info) << "CbmTrdRadiator::Init : Foil thickness : " << setprecision(4) << fFoilThick <<" cm"; - LOG(info) << "CbmTrdRadiator::Init : Gap thickness : " << fGapThick << " cm"; - LOG(info) << "CbmTrdRadiator::Init : Simple TR production: " << setprecision(2)<< fSimpleTR; - */ - // material dependend parameters for the radiator material - // (polyethylen) - fFoilDens = 0.92; // [g/cm3 ] - fFoilOmega = 20.9; //plasma frequency ( from Anton ) +//---------------------------------------------------------------------------- +CbmTrdRadiator::CbmTrdEntranceWindow::CbmTrdEntranceWindow(TString def) : fN(0) +{ + memset(fThick, 0, NCOMPONENTS * sizeof(Float_t)); + memset(fDens, 0, NCOMPONENTS * sizeof(Float_t)); + memset(GetMu, 0, NCOMPONENTS * sizeof(GetMuPtr)); + Init(def); +} - // material dependent parameters for the gas between the foils of the - // radiator - // changed gap density at 16.04.07 FU - // density of air is 0.001205, 0.00186 - //fGapDens = 0.00186; // [g/cm3] - fGapDens = 0.001205; // [g/cm3] - fGapOmega = 0.7; //plasma frequency ( from Anton ) - - // foil between radiator and gas chamber - // TODO: implement kapton foil also in trd geometry - fWinDens = 1.39; // [g/cm3] - fWinThick = 0.0025; // [cm] - } else { - if (fRadType == "default") { - material = "polyethylen"; - fNFoils = 130; - fGapThick = 0.02; - fFoilThick = 0.0013; - fnPhotonCorr = 1.0; - } else if (fRadType == "A") { - material = "polyethylen"; - fNFoils = 200; - fGapThick = 700 * 0.0001; - fFoilThick = 15 * 0.0001; - fnPhotonCorr = 0.40; - } else if (fRadType == "Bshort" || fRadType == "Kshort") { - material = "pokalon"; - fNFoils = 100; - fGapThick = 0.07; - fFoilThick = 0.0024; - fnPhotonCorr = 0.65; - } else if (fRadType == "B" || fRadType == "K") { - material = "pokalon"; - fNFoils = 250; - fGapThick = 0.07; - fFoilThick = 0.0024; - fnPhotonCorr = 0.65; - } else if (fRadType == "B++" || fRadType == "K++") { - material = "pokalon"; - fNFoils = 350; - fGapThick = 0.07; - fFoilThick = 0.0024; - fnPhotonCorr = 0.65; - } else if (fRadType == "C") { - material = "polyethylen"; - fNFoils = 200; - fGapThick = 700 * 0.0001; - fFoilThick = 15 * 0.0001; - fnPhotonCorr = 0.30; - } else if (fRadType == "D") { - material = "polyethylen"; - fNFoils = 100; - fGapThick = 500 * 0.0001; - fFoilThick = 15 * 0.0001; - fnPhotonCorr = 0.45; - } else if (fRadType == "E") { - material = "polyethylen"; - fNFoils = 120; - fGapThick = 500 * 0.0001; - fFoilThick = 20 * 0.0001; - fnPhotonCorr = 0.60; - } else if (fRadType == "F") { - material = "polyethylen"; - fNFoils = 220; - fGapThick = 250 * 0.0001; - fFoilThick = 20 * 0.0001; - fnPhotonCorr = 0.65; - } else if (fRadType == "G30") { - material = "pefiber"; - fNFoils = 200; - fGapThick = 700 * 0.0001; - fFoilThick = 15 * 0.0001; - fnPhotonCorr = 0.65; - } else if (fRadType == "H") { - material = "pefoam20"; - fNFoils = 200; - fGapThick = 700 * 0.0001; - fFoilThick = 15 * 0.0001; - fnPhotonCorr = 0.65; - } else if (fRadType == "H++") { - material = "pefoam20"; - fNFoils = 200; - fGapThick = 700 * 0.0001; - fFoilThick = 15 * 0.0001; - fnPhotonCorr = 0.65; - } else { - LOG(warn) << "CbmTrdRadiator::Init : *********** Radiator prototype not " - "known ***********"; - LOG(warn) << "CbmTrdRadiator::Init : switch to default " - "parameters "; - material = "polyethylen"; - fNFoils = 130; - fGapThick = 0.02; - fFoilThick = 0.0013; - fnPhotonCorr = 1.0; +//---------------------------------------------------------------------------- +Bool_t CbmTrdRadiator::CbmTrdEntranceWindow::Init(TString def) +{ + TObjArray* mats = def.Tokenize(";"); + for (fN = 0; fN < mats->GetEntriesFast(); fN++) { + TObjString* mat = (TObjString*) (*mats)[fN]; + TString mname = mat->String(); + if (mname.CompareTo("Al") == 0) { + fMat[fN] = "Al"; + fDens[fN] = 2.7; + GetMu[fN] = CbmTrdRadiator::GetMuAl; } - fFoilMaterial = material; - if (material == "" || material == "polyethylen") { - fFoilDens = 0.92; // [g/cm3 ] - fFoilOmega = 20.9; //plasma frequency ( from Anton ) - } else if (material == "pefoam20") { - // (polyethylen) - fFoilDens = 0.90; // [g/cm3 ] - fFoilOmega = 20.64; //plasma frequency ( from Cyrano ) - } else if (material == "pefiber") { - // (polyethylen) - fFoilDens = 0.90; // [g/cm3 ] - fFoilOmega = 20.64; //plasma frequency ( from Cyrano ) - } else if (material == "mylar") { - // (Mylar) - fFoilDens = 1.393; // [g/cm3 ] - fFoilOmega = 24.53; //plasma frequency ( from Cyrano ) - } else if (material == "pokalon") { - // (Pokalon) - fFoilDens = 1.150; // [g/cm3 ] - fFoilOmega = 22.50; //plasma frequency ( from Cyrano ) - } else { - // (polyethylen) - fFoilDens = 0.92; // [g/cm3 ] - fFoilOmega = 20.9; //plasma frequency ( from Anton ) + else if (mname.CompareTo("C") == 0) { + fMat[fN] = "C"; + fDens[fN] = 2.267; + GetMu[fN] = CbmTrdRadiator::GetMuC; } - fGapDens = 0.001205; // [g/cm3] - fGapOmega = 0.7; //plasma frequency ( from Anton ) - // foil between radiator and gas chamber - fWinDens = - 1.42; //[g/cm3] for Kapton -> Alu is added only in case of kapton !!! - fWinThick = 0.0025; // [cm] - - LOG(info) << "CbmTrdRadiator::Init : Prototype : " << fRadType; - LOG(info) << "CbmTrdRadiator::Init : Scaling factor : " - << fnPhotonCorr; - LOG(info) << "CbmTrdRadiator::Init : Foil material : " - << fFoilMaterial; - LOG(info) << "CbmTrdRadiator::Init : Nr. of foils : " << fNFoils; - LOG(info) << "CbmTrdRadiator::Init : Foil thickness : " - << setprecision(4) << fFoilThick << " cm"; - LOG(info) << "CbmTrdRadiator::Init : Gap thickness : " << fGapThick - << " cm"; - LOG(info) << "CbmTrdRadiator::Init : Simple TR production: " - << setprecision(2) << fSimpleTR; - LOG(info) << "CbmTrdRadiator::Init : Foil plasm. frequ. : " << fFoilOmega - << " keV"; - LOG(info) << "CbmTrdRadiator::Init : Gap plasm. frequ. : " << fGapOmega - << " keV"; - LOG(info) << "CbmTrdRadiator::Init : Window density : " << fWinDens - << " g/cm^3"; - LOG(info) << "CbmTrdRadiator::Init : Window thickness : " << fWinThick - << " cm"; - } - - - // Get all the gas properties from CbmTrdGas - CbmTrdGas* fTrdGas = CbmTrdGas::Instance(); - if (fTrdGas == 0) { - fTrdGas = new CbmTrdGas(); - fTrdGas->Init(); - } - - fCom2 = fTrdGas->GetNobleGas(); - fCom1 = fTrdGas->GetCO2(); - fDetType = fTrdGas->GetDetType(); - fGasThick = fTrdGas->GetGasThick(); - - LOG(info) << "CbmTrdRadiator::Init : Detector noble gas : " - << fTrdGas->GetNobleGas(); - LOG(info) << "CbmTrdRadiator::Init : Detector cruncher gas: " - << fTrdGas->GetCO2(); - LOG(info) << "CbmTrdRadiator::Init : Detector type : " - << fTrdGas->GetDetType(); - LOG(info) << "CbmTrdRadiator::Init : Detector gas thick. : " - << fTrdGas->GetGasThick() << " cm"; - - //LOG(info) << "************* End of Trd Radiator Init **************"; - // If simplified version is used one didn't calculate the TR - // for each CbmTrdPoint in full glory, but create at startup - // some histograms for several momenta which are later used to - // get the TR much faster. - if (fSimpleTR == kTRUE) { - - ProduceSpectra(); - - TFile* f1 = new TFile("TRhistos.root", "recreate"); - - for (Int_t i = 0; i < fNMom; i++) { - fFinal[i]->Write(); + else if (mname.CompareTo("HC") == 0) { + fMat[fN] = "HC"; + fDens[fN] = 4.8e-2; // 48 Kg/m^3 + GetMu[fN] = CbmTrdRadiator::GetMuAir; } - f1->Close(); - f1->Delete(); + else if (mname.CompareTo("Po") == 0) { + fMat[fN] = "Po"; + fDens[fN] = 0.92; + GetMu[fN] = CbmTrdRadiator::GetMuPo; + } + else if (mname.CompareTo("Pok") == 0) { + fMat[fN] = "Pok"; + fDens[fN] = 1.150; + GetMu[fN] = CbmTrdRadiator::GetMuPok; + } + else if (mname.CompareTo("Ka") == 0) { + fMat[fN] = "Ka"; + fDens[fN] = 1.39; + GetMu[fN] = CbmTrdRadiator::GetMuKa; + } + else if (mname.CompareTo("Air") == 0) { + fMat[fN] = "Air"; + fDens[fN] = 1.205e-3; + GetMu[fN] = CbmTrdRadiator::GetMuAir; + } + else if (mname.CompareTo("Xe") == 0) { + fMat[fN] = "Xe"; + fDens[fN] = 5.9e-3; + GetMu[fN] = CbmTrdRadiator::GetMuXe; + } + else if (mname.CompareTo("CO2") == 0) { + fMat[fN] = "CO2"; + fDens[fN] = 1.9767e-3; + GetMu[fN] = CbmTrdRadiator::GetMuCO2; + } + else if (mname.CompareTo("My") == 0) { + fMat[fN] = "My"; + fDens[fN] = 1.393; + GetMu[fN] = CbmTrdRadiator::GetMuMy; + } + else { + LOG(warn) << "CbmTrdEntranceWindow::Init : material " << mname.Data() << " not in the DB. Can't create " + << def.Data() << " entrance window structure. Default to Mylar."; + fN = 0; + return kFALSE; + } + LOG(info) << "CbmTrdEntranceWindow::Init : C" << fN << " type : " << fMat[fN]; + LOG(info) << "CbmTrdEntranceWindow::Init : density : " << fDens[fN] << " g/cm^3"; + LOG(info) << "CbmTrdEntranceWindow::Init : thickness : " << fThick[fN] << " cm"; } + mats->Delete(); + delete mats; + return kTRUE; } + //---------------------------------------------------------------------------- +void CbmTrdRadiator::SetEWwidths(Int_t n, Float_t* w) +{ + if (n > NCOMPONENTS) { + LOG(warn) << "CbmTrdEntranceWindow : can support only up to " << NCOMPONENTS + << " plane-parallel elements. Please check your simulations"; + n = NCOMPONENTS; + } + memcpy(WINDOW.fThick, w, n * sizeof(Float_t)); +} //---- Lattice Hit------------------------------------------------------------ -Bool_t CbmTrdRadiator::LatticeHit(const CbmTrdPoint* point) { +Bool_t CbmTrdRadiator::LatticeHit(const CbmTrdPoint* point) +{ //printf("---------------------------------------------------\n"); Double_t point_in[3]; Double_t point_out[3]; @@ -627,28 +432,23 @@ Bool_t CbmTrdRadiator::LatticeHit(const CbmTrdPoint* point) { point_out[0] = point->GetXOut(); point_out[1] = point->GetYOut(); point_out[2] = point->GetZOut(); - Double_t back_direction[3] = {(point_in[0] - point_out[0]), - (point_in[1] - point_out[1]), + Double_t back_direction[3] = {(point_in[0] - point_out[0]), (point_in[1] - point_out[1]), (point_in[2] - point_out[2])}; if (back_direction[2] > 0) { - LOG(debug2) - << "CbmTrdRadiator::LatticeHit: MC-track points towards target!"; + LOG(debug2) << "CbmTrdRadiator::LatticeHit: MC-track points towards target!"; //for(Int_t i = 0; i < 3; i++) //printf("%i: in:%8.2f out:%8.2f direction:%8.2f\n",i,point_in[i],point_out[i],back_direction[i]); return false; } - Double_t trackLength = TMath::Sqrt(back_direction[0] * back_direction[0] - + back_direction[1] * back_direction[1] + Double_t trackLength = TMath::Sqrt(back_direction[0] * back_direction[0] + back_direction[1] * back_direction[1] + back_direction[2] * back_direction[2]); trackLength *= 10.; // cm -> mm to get a step width of 1mm //printf("track length:%7.2fmm\n",trackLength); - gGeoManager->FindNode((point_out[0] + point_in[0]) / 2, - (point_out[1] + point_in[1]) / 2, + gGeoManager->FindNode((point_out[0] + point_in[0]) / 2, (point_out[1] + point_in[1]) / 2, (point_out[2] + point_in[2]) / 2); - Double_t pos[3] = { - point_in[0], point_in[1], point_in[2]}; // start at entrance point + Double_t pos[3] = {point_in[0], point_in[1], point_in[2]}; // start at entrance point for (Int_t i = 0; i < 3; i++) { back_direction[i] /= trackLength; @@ -658,8 +458,7 @@ Bool_t CbmTrdRadiator::LatticeHit(const CbmTrdPoint* point) { Int_t stepCount = 0; while ( /*(!TString(gGeoManager->GetPath()).Contains("lattice") || !TString(gGeoManager->GetPath()).Contains("radiator")) &&*/ - pos[2] >= point_in[2] - 5 - && stepCount < 50) { + pos[2] >= point_in[2] - 5 && stepCount < 50) { stepCount++; //printf("step%i\n",stepCount); for (Int_t i = 0; i < 3; i++) { @@ -671,23 +470,27 @@ Bool_t CbmTrdRadiator::LatticeHit(const CbmTrdPoint* point) { if (TString(gGeoManager->GetPath()).Contains("radiator")) { //printf ("%s false\n",TString(gGeoManager->GetPath()).Data()); return false; - } else if (TString(gGeoManager->GetPath()).Contains("lattice")) { + } + else if (TString(gGeoManager->GetPath()).Contains("lattice")) { //printf ("%s true <----------------------------------------\n",TString(gGeoManager->GetPath()).Data()); return true; - } else if (TString(gGeoManager->GetPath()).Contains("frame")) { + } + else if (TString(gGeoManager->GetPath()).Contains("frame")) { //printf ("%s true <----------------------------------------\n",TString(gGeoManager->GetPath()).Data()); return true; - } else { + } + else { //printf ("%s\n",TString(gGeoManager->GetPath()).Data()); } } - } else { - LOG(error) << "CbmTrdRadiator::LatticeHit: MC-track not in TRD! Node:" - << TString(gGeoManager->GetPath()).Data() + } + else { + LOG(error) << "CbmTrdRadiator::LatticeHit: MC-track not in TRD! Node:" << TString(gGeoManager->GetPath()).Data() << " gGeoManager->MasterToLocal() failed!"; return false; } - } else { + } + else { LOG(error) << "CbmTrdRadiator::LatticeHit: CbmTrdPoint == nullptr!"; return false; } @@ -695,12 +498,12 @@ Bool_t CbmTrdRadiator::LatticeHit(const CbmTrdPoint* point) { } //---------------------------------------------------------------------------- // ---- Spectra Production --------------------------------------------------- -void CbmTrdRadiator::ProduceSpectra() { +void CbmTrdRadiator::ProduceSpectra() +{ fTrackMomentum = new Double_t[fNMom]; - Double_t trackMomentum[fNMom] = { - 0.1, 0.25, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0}; + Double_t trackMomentum[fNMom] = {0.1, 0.25, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0}; for (Int_t imom = 0; imom < fNMom; imom++) { fTrackMomentum[imom] = trackMomentum[imom]; @@ -721,15 +524,15 @@ void CbmTrdRadiator::ProduceSpectra() { tmp = fDetSpectrumA->GetBinContent(j + 1); fFinal[i]->SetBinContent(j + 1, tmp); } - fFinal[i]->SetTitle( - Form("%s-momentum%.2f", fFinal[i]->GetTitle(), fTrackMomentum[i])); + fFinal[i]->SetTitle(Form("%s-momentum%.2f", fFinal[i]->GetTitle(), fTrackMomentum[i])); } } //---------------------------------------------------------------------------- // ---- GetTR --------------------------------------------------------------- -Float_t CbmTrdRadiator::GetTR(TVector3 Mom) { +Float_t CbmTrdRadiator::GetTR(TVector3 Mom) +{ /* * Simplified simulation of the TR @@ -750,7 +553,6 @@ Float_t CbmTrdRadiator::GetTR(TVector3 Mom) { } ELoss(i); - } /* @@ -768,7 +570,8 @@ Float_t CbmTrdRadiator::GetTR(TVector3 Mom) { //---------------------------------------------------------------------------- //----- main TR processing function ------------------------------ -void CbmTrdRadiator::ProcessTR() { +void CbmTrdRadiator::ProcessTR() +{ // Compute the angle normalization factor - // for different angles the detector thicknesses are different @@ -810,7 +613,8 @@ void CbmTrdRadiator::ProcessTR() { //---------------------------------------------------------------------------- // ---- Compute the energy losses -------------------------------------------- -Int_t CbmTrdRadiator::ELoss(Int_t index) { +Int_t CbmTrdRadiator::ELoss(Int_t index) +{ // calculate the energy loss of the TR photons in the // detector gas. fNTRab is the number of TR photons @@ -835,7 +639,8 @@ Int_t CbmTrdRadiator::ELoss(Int_t index) { eTR = fDetSpectrumA->GetRandom(); eLoss += eTR; } - } else { + } + else { if (0 == fFinal[index]) { LOG(error) << " -Err :: No input spectra "; fELoss = -1; @@ -859,7 +664,8 @@ Int_t CbmTrdRadiator::ELoss(Int_t index) { //---------------------------------------------------------------------------- //----- TR spectrum --------------------------------------------------- -Int_t CbmTrdRadiator::TRspectrum() { +Int_t CbmTrdRadiator::TRspectrum() +{ // calculate the number of TR photons in the radiator which are not // absorbed in the radiator. @@ -894,16 +700,13 @@ Int_t CbmTrdRadiator::TRspectrum() { Float_t csFoil = fFoilOmega / energyeV; Float_t csGap = fGapOmega / energyeV; - Float_t rho1 = energyeV * fFoilThickCorr * 1e4 * 2.5 - * (1.0 / (gamma * gamma) + csFoil * csFoil); - Float_t rho2 = energyeV * fFoilThickCorr * 1e4 * 2.5 - * (1.0 / (gamma * gamma) + csGap * csGap); + Float_t rho1 = energyeV * fFoilThickCorr * 1e4 * 2.5 * (1.0 / (gamma * gamma) + csFoil * csFoil); + Float_t rho2 = energyeV * fFoilThickCorr * 1e4 * 2.5 * (1.0 / (gamma * gamma) + csGap * csGap); // Calculate the sum over the production angle tetan Float_t sum = 0; for (Int_t iSum = 0; iSum < kSumMax; iSum++) { - Float_t tetan = (TMath::Pi() * 2.0 * (iSum + 1) - (rho1 + kappa * rho2)) - / (kappa + 1.0); + Float_t tetan = (TMath::Pi() * 2.0 * (iSum + 1) - (rho1 + kappa * rho2)) / (kappa + 1.0); if (tetan < 0.0) { tetan = 0.0; } Float_t aux = 1.0 / (rho1 + tetan) - 1.0 / (rho2 + tetan); sum += tetan * (aux * aux) * (1.0 - TMath::Cos(rho1 + tetan)); @@ -914,8 +717,7 @@ Int_t CbmTrdRadiator::TRspectrum() { Float_t energykeV = energyeV * 0.001; // dN / domega - Float_t wn = - kAlpha * 4.0 / (fSigma[iBin] * (kappa + 1.0)) * conv * sum / energykeV; + Float_t wn = kAlpha * 4.0 / (fSigma[iBin] * (kappa + 1.0)) * conv * sum / energykeV; /* wn = 4.0 * kAlpha / (energykeV * (kappa + 1.0)) @@ -941,7 +743,8 @@ Int_t CbmTrdRadiator::TRspectrum() { //------------------------------------------------------------------ //----- WinTRspectrum ----------------------------------------------- -Int_t CbmTrdRadiator::WinTRspectrum() { +Int_t CbmTrdRadiator::WinTRspectrum() +{ // // Computes the spectrum after absorption in Mylar foil // @@ -954,9 +757,8 @@ Int_t CbmTrdRadiator::WinTRspectrum() { for (Int_t iBin = 0; iBin < fSpNBins; iBin++) { Float_t sp = 0; - if (fFirstPass == true) { - sp = fSpectrum->GetBinContent(iBin + 1); - } else { + if (fFirstPass == true) { sp = fSpectrum->GetBinContent(iBin + 1); } + else { sp = fDetSpectrum->GetBinContent(iBin + 1); } // compute the absorption coefficient @@ -978,7 +780,8 @@ Int_t CbmTrdRadiator::WinTRspectrum() { //---------------------------------------------------------------------------- //----- DetTRspectrum ------------------------------------------------ -Int_t CbmTrdRadiator::DetTRspectrum() { +Int_t CbmTrdRadiator::DetTRspectrum() +{ // // Computes the spectrum absorbed and passed in the gas mixture // @@ -1019,7 +822,8 @@ Int_t CbmTrdRadiator::DetTRspectrum() { //---------------------------------------------------------------------------- //----- Gamma factor-------------------------------------------------- -Float_t CbmTrdRadiator::GammaF() { +Float_t CbmTrdRadiator::GammaF() +{ // put this in the header ? Float_t massE = 5.1099907e-4; // GeV/c2 @@ -1032,7 +836,8 @@ Float_t CbmTrdRadiator::GammaF() { //---------------------------------------------------------------------------- //----- SetSigma--------------------------------------------------------- -void CbmTrdRadiator::SetSigma(Int_t SigmaT) { +void CbmTrdRadiator::SetSigma(Int_t SigmaT) +{ // // Sets the absorbtion crosssection for the energies of the TR spectrum // @@ -1068,7 +873,8 @@ void CbmTrdRadiator::SetSigma(Int_t SigmaT) { //---------------------------------------------------------------------------- //----- Sigma------------------------------------------------------------- -Float_t CbmTrdRadiator::Sigma(Float_t energykeV) { +Float_t CbmTrdRadiator::Sigma(Float_t energykeV) +{ // // Calculates the absorbtion crosssection for a one-foil-one-gap-radiator // @@ -1077,8 +883,7 @@ Float_t CbmTrdRadiator::Sigma(Float_t energykeV) { Float_t energyMeV = energykeV * 0.001; Float_t foil = 0.0; - if (fFoilMaterial == "polyethylen") - foil = GetMuPo(energyMeV) * fFoilDens * fFoilThickCorr; + if (fFoilMaterial == "polyethylen") foil = GetMuPo(energyMeV) * fFoilDens * fFoilThickCorr; else if (fFoilMaterial == "pefoam20") foil = GetMuPo(energyMeV) * fFoilDens * fFoilThickCorr; else if (fFoilMaterial == "pefiber") @@ -1094,40 +899,42 @@ Float_t CbmTrdRadiator::Sigma(Float_t energykeV) { if (energyMeV >= 0.001) { Float_t result = (foil + (GetMuAir(energyMeV) * fGapDens * fGapThickCorr)); return result; - } else { + } + else { return 1e6; } } //---------------------------------------------------------------------------- //----- SigmaWin ------------------------------------------------------- -Float_t CbmTrdRadiator::SigmaWin(Float_t energykeV) { +Float_t CbmTrdRadiator::SigmaWin(Float_t energykeV) +{ // // Calculates the absorbtion crosssection for a one-foil // // keV -> MeV Float_t energyMeV = energykeV * 0.001; - //if (fWindowFoil=="Kapton"){ - //if (energyMeV >= 0.001) { - //return((GetMuKa(energyMeV) * fWinDens * fWinThick) + (GetMuAl(energyMeV) * 2.70/*[g/cm^3]*/ * 5E-6/*[cm]*/)); - //} - // else { - // return 1e6; - // } - // } - // else { if (energyMeV >= 0.001) { - return (GetMuMy(energyMeV) * fWinDens * fWinThick); - } else { - return 1e6; + if (fWindowFoil == "Kapton") // aluminized kapton + return GetMuKa(energyMeV) * fWinDens * fWinThick + GetMuAl(energyMeV) * 2.70 /*[g/cm^3]*/ * 5E-6 /*[cm]*/; + else if (fWindowFoil == "Mylar") // mylar foil + return GetMuMy(energyMeV) * fWinDens * fWinThick; + else { // general sandwich structure + Float_t sigma(0.); + for (Int_t iwin(0); iwin < WINDOW.fN; iwin++) + sigma += WINDOW.GetMu[iwin](energyMeV) * WINDOW.fDens[iwin] * WINDOW.fThick[iwin]; + return sigma; + } } - // } + else + return 1e6; } //---------------------------------------------------------------------------- //----- SigmaDet -------------------------------------------------------- -Float_t CbmTrdRadiator::SigmaDet(Float_t energykeV) { +Float_t CbmTrdRadiator::SigmaDet(Float_t energykeV) +{ // //Calculates the absorbtion crosssection for a choosed gas // @@ -1141,257 +948,106 @@ Float_t CbmTrdRadiator::SigmaDet(Float_t energykeV) { // Values are from http://pdg.lbl.gov/AtomicNuclearProperties/ // return(GetMuCO2(energyMeV) * 0.001806 * fGasThick * fCom1 + // GetMuXe(energyMeV) * 0.00589 * fGasThick * fCom2); - return (GetMuCO2(energyMeV) * 0.00184 * fGasThick * fCom1 - + GetMuXe(energyMeV) * 0.00549 * fGasThick * fCom2); - } else { + return (GetMuCO2(energyMeV) * 0.00184 * fGasThick * fCom1 + GetMuXe(energyMeV) * 0.00549 * fGasThick * fCom2); + } + else { return 1e6; } } //---------------------------------------------------------------------------- //----- GetMuPok -------------------------------------------------------- -Float_t CbmTrdRadiator::GetMuPok(Float_t energyMeV) { +Float_t CbmTrdRadiator::GetMuPok(Float_t energyMeV) +{ // // Returns the photon absorbtion cross section for pokalon N470 // - const Int_t kN = 46; - Float_t en[kN] = { - 1.000E-03, 1.500E-03, 2.000E-03, 3.000E-03, 4.000E-03, 5.000E-03, 6.000E-03, - 8.000E-03, 1.000E-02, 1.500E-02, 2.000E-02, 3.000E-02, 4.000E-02, 5.000E-02, - 6.000E-02, 8.000E-02, 1.000E-01, 1.500E-01, 2.000E-01, 3.000E-01, 4.000E-01, - 5.000E-01, 6.000E-01, 8.000E-01, 1.000E+00, 1.022E+00, 1.250E+00, 1.500E+00, - 2.000E+00, 2.044E+00, 3.000E+00, 4.000E+00, 5.000E+00, 6.000E+00, 7.000E+00, - 8.000E+00, 9.000E+00, 1.000E+01, 1.100E+01, 1.200E+01, 1.300E+01, 1.400E+01, - 1.500E+01, 1.600E+01, 1.800E+01, 2.000E+01}; - Float_t mu[kN] = { - 2.537E+03, 8.218E+02, 3.599E+02, 1.093E+02, 4.616E+01, 2.351E+01, 1.352E+01, - 5.675E+00, 2.938E+00, 9.776E-01, 5.179E-01, 2.847E-01, 2.249E-01, 2.003E-01, - 1.866E-01, 1.705E-01, 1.600E-01, 1.422E-01, 1.297E-01, 1.125E-01, 1.007E-01, - 9.193E-02, 8.501E-02, 7.465E-02, 6.711E-02, 6.642E-02, 6.002E-02, 5.463E-02, - 4.686E-02, 4.631E-02, 3.755E-02, 3.210E-02, 2.851E-02, 2.597E-02, 2.409E-02, - 2.263E-02, 2.149E-02, 2.056E-02, 1.979E-02, 1.916E-02, 1.862E-02, 1.816E-02, - 1.777E-02, 1.743E-02, 1.687E-02, 1.644E-02}; - return Interpolate(energyMeV, en, mu, kN); + return Interpolate(energyMeV, fEn_pok, fMu_pok, fKN_pok); } //---------------------------------------------------------------------------- //----- GetMuKa -------------------------------------------------------- -Float_t CbmTrdRadiator::GetMuKa(Float_t energyMeV) { +Float_t CbmTrdRadiator::GetMuKa(Float_t energyMeV) +{ // // Returns the photon absorbtion cross section for kapton // - const Int_t kN = 46; - Float_t en[kN] = { - 1.000E-03, 1.500E-03, 2.000E-03, 3.000E-03, 4.000E-03, 5.000E-03, 6.000E-03, - 8.000E-03, 1.000E-02, 1.500E-02, 2.000E-02, 3.000E-02, 4.000E-02, 5.000E-02, - 6.000E-02, 8.000E-02, 1.000E-01, 1.500E-01, 2.000E-01, 3.000E-01, 4.000E-01, - 5.000E-01, 6.000E-01, 8.000E-01, 1.000E+00, 1.022E+00, 1.250E+00, 1.500E+00, - 2.000E+00, 2.044E+00, 3.000E+00, 4.000E+00, 5.000E+00, 6.000E+00, 7.000E+00, - 8.000E+00, 9.000E+00, 1.000E+01, 1.100E+01, 1.200E+01, 1.300E+01, 1.400E+01, - 1.500E+01, 1.600E+01, 1.800E+01, 2.000E+01}; - Float_t mu[kN] = { - 2.731E+03, 8.875E+02, 3.895E+02, 1.185E+02, 5.013E+01, 2.555E+01, 1.470E+01, - 6.160E+00, 3.180E+00, 1.043E+00, 5.415E-01, 2.880E-01, 2.235E-01, 1.973E-01, - 1.830E-01, 1.666E-01, 1.560E-01, 1.385E-01, 1.263E-01, 1.095E-01, 9.799E-02, - 8.944E-02, 8.270E-02, 7.262E-02, 6.529E-02, 6.462E-02, 5.839E-02, 5.315E-02, - 4.561E-02, 4.507E-02, 3.659E-02, 3.133E-02, 2.787E-02, 2.543E-02, 2.362E-02, - 2.223E-02, 2.114E-02, 2.025E-02, 1.953E-02, 1.893E-02, 1.842E-02, 1.799E-02, - 1.762E-02, 1.730E-02, 1.678E-02, 1.638E-02}; - return Interpolate(energyMeV, en, mu, kN); + return Interpolate(energyMeV, fEn_ka, fMu_ka, fKN_ka); } //---------------------------------------------------------------------------- //----- GetMuAl -------------------------------------------------------- -Float_t CbmTrdRadiator::GetMuAl(Float_t energyMeV) { +Float_t CbmTrdRadiator::GetMuAl(Float_t energyMeV) +{ // // Returns the photon absorbtion cross section for Al // - - const Int_t kN = 48; - - Float_t en[kN] = { - 1.000E-03, 1.500E-03, 1.560E-03, 1.560E-03, 2.000E-03, 3.000E-03, 4.000E-03, - 5.000E-03, 6.000E-03, 8.000E-03, 1.000E-02, 1.500E-02, 2.000E-02, 3.000E-02, - 4.000E-02, 5.000E-02, 6.000E-02, 8.000E-02, 1.000E-01, 1.500E-01, 2.000E-01, - 3.000E-01, 4.000E-01, 5.000E-01, 6.000E-01, 8.000E-01, 1.000E+00, 1.022E+00, - 1.250E+00, 1.500E+00, 2.000E+00, 2.044E+00, 3.000E+00, 4.000E+00, 5.000E+00, - 6.000E+00, 7.000E+00, 8.000E+00, 9.000E+00, 1.000E+01, 1.100E+01, 1.200E+01, - 1.300E+01, 1.400E+01, 1.500E+01, 1.600E+01, 1.800E+01, 2.000E+01}; - - Float_t mu[kN] = { - 1.185E+03, 4.023E+02, 3.621E+02, 3.957E+03, 2.263E+03, 7.881E+02, 3.605E+02, - 1.934E+02, 1.153E+02, 5.032E+01, 2.621E+01, 7.955E+00, 3.442E+00, 1.128E+00, - 5.684E-01, 3.681E-01, 2.778E-01, 2.018E-01, 1.704E-01, 1.378E-01, 1.223E-01, - 1.042E-01, 9.276E-02, 8.445E-02, 7.802E-02, 6.841E-02, 6.146E-02, 6.080E-02, - 5.496E-02, 5.006E-02, 4.324E-02, 4.277E-02, 3.541E-02, 3.106E-02, 2.836E-02, - 2.655E-02, 2.529E-02, 2.437E-02, 2.369E-02, 2.318E-02, 2.279E-02, 2.249E-02, - 2.226E-02, 2.208E-02, 2.195E-02, 2.185E-02, 2.173E-02, 2.168E-02}; - - return Interpolate(energyMeV, en, mu, kN); + return Interpolate(energyMeV, fEn_al, fMu_al, fKN_al); } //---------------------------------------------------------------------------- +//----- GetMuC -------------------------------------------------------- +Float_t CbmTrdRadiator::GetMuC(Float_t energyMeV) +{ + // + // Returns the photon absorbtion cross section for Carbon + // + return Interpolate(energyMeV, fEn_c, fMu_c, fKN_c); +} +//---------------------------------------------------------------------------- //----- GetMuPo -------------------------------------------------------- -Float_t CbmTrdRadiator::GetMuPo(Float_t energyMeV) { +Float_t CbmTrdRadiator::GetMuPo(Float_t energyMeV) +{ // // Returns the photon absorbtion cross section for polypropylene // - - const Int_t kN = 36; - - Float_t mu[kN] = { - 1.894E+03, 5.999E+02, 2.593E+02, 7.743E+01, 3.242E+01, 1.643E+01, - 9.432E+00, 3.975E+00, 2.088E+00, 7.452E-01, 4.315E-01, 2.706E-01, - 2.275E-01, 2.084E-01, 1.970E-01, 1.823E-01, 1.719E-01, 1.534E-01, - 1.402E-01, 1.217E-01, 1.089E-01, 9.947E-02, 9.198E-02, 8.078E-02, - 7.262E-02, 6.495E-02, 5.910E-02, 5.064E-02, 4.045E-02, 3.444E-02, - 3.045E-02, 2.760E-02, 2.383E-02, 2.145E-02, 1.819E-02, 1.658E-02}; - - Float_t en[kN] = { - 1.000E-03, 1.500E-03, 2.000E-03, 3.000E-03, 4.000E-03, 5.000E-03, - 6.000E-03, 8.000E-03, 1.000E-02, 1.500E-02, 2.000E-02, 3.000E-02, - 4.000E-02, 5.000E-02, 6.000E-02, 8.000E-02, 1.000E-01, 1.500E-01, - 2.000E-01, 3.000E-01, 4.000E-01, 5.000E-01, 6.000E-01, 8.000E-01, - 1.000E+00, 1.250E+00, 1.500E+00, 2.000E+00, 3.000E+00, 4.000E+00, - 5.000E+00, 6.000E+00, 8.000E+00, 1.000E+01, 1.500E+01, 2.000E+01}; - - return Interpolate(energyMeV, en, mu, kN); + return Interpolate(energyMeV, fEn_po, fMu_po, fKN_po); } //---------------------------------------------------------------------------- //----- GetMuAir -------------------------------------------------------- -Float_t CbmTrdRadiator::GetMuAir(Float_t energyMeV) { +Float_t CbmTrdRadiator::GetMuAir(Float_t energyMeV) +{ // // Returns the photon absorbtion cross section for Air // - - const Int_t kN = 38; - - - Float_t mu[kN] = { - 0.35854E+04, 0.11841E+04, 0.52458E+03, 0.16143E+03, 0.15722E+03, - 0.14250E+03, 0.77538E+02, 0.40099E+02, 0.23313E+02, 0.98816E+01, - 0.51000E+01, 0.16079E+01, 0.77536E+00, 0.35282E+00, 0.24790E+00, - 0.20750E+00, 0.18703E+00, 0.16589E+00, 0.15375E+00, 0.13530E+00, - 0.12311E+00, 0.10654E+00, 0.95297E-01, 0.86939E-01, 0.80390E-01, - 0.70596E-01, 0.63452E-01, 0.56754E-01, 0.51644E-01, 0.44382E-01, - 0.35733E-01, 0.30721E-01, 0.27450E-01, 0.25171E-01, 0.22205E-01, - 0.20399E-01, 0.18053E-01, 0.18057E-01}; - - Float_t en[kN] = { - 0.10000E-02, 0.15000E-02, 0.20000E-02, 0.30000E-02, 0.32029E-02, - 0.32029E-02, 0.40000E-02, 0.50000E-02, 0.60000E-02, 0.80000E-02, - 0.10000E-01, 0.15000E-01, 0.20000E-01, 0.30000E-01, 0.40000E-01, - 0.50000E-01, 0.60000E-01, 0.80000E-01, 0.10000E+00, 0.15000E+00, - 0.20000E+00, 0.30000E+00, 0.40000E+00, 0.50000E+00, 0.60000E+00, - 0.80000E+00, 0.10000E+01, 0.12500E+01, 0.15000E+01, 0.20000E+01, - 0.30000E+01, 0.40000E+01, 0.50000E+01, 0.60000E+01, 0.80000E+01, - 0.10000E+02, 0.15000E+02, 0.20000E+02}; - - return Interpolate(energyMeV, en, mu, kN); + return Interpolate(energyMeV, fEn_air, fMu_air, fKN_air); } //---------------------------------------------------------------------------- //----- GetMuXe -------------------------------------------------------- -Float_t CbmTrdRadiator::GetMuXe(Float_t energyMeV) { +Float_t CbmTrdRadiator::GetMuXe(Float_t energyMeV) +{ // // Returns the photon absorbtion cross section for xenon // - - const Int_t kN = 48; - - Float_t mu[kN] = { - 9.413E+03, 8.151E+03, 7.035E+03, 7.338E+03, 4.085E+03, 2.088E+03, 7.780E+02, - 3.787E+02, 2.408E+02, 6.941E+02, 6.392E+02, 6.044E+02, 8.181E+02, 7.579E+02, - 6.991E+02, 8.064E+02, 6.376E+02, 3.032E+02, 1.690E+02, 5.743E+01, 2.652E+01, - 8.930E+00, 6.129E+00, 3.316E+01, 2.270E+01, 1.272E+01, 7.825E+00, 3.633E+00, - 2.011E+00, 7.202E-01, 3.760E-01, 1.797E-01, 1.223E-01, 9.699E-02, 8.281E-02, - 6.696E-02, 5.785E-02, 5.054E-02, 4.594E-02, 4.078E-02, 3.681E-02, 3.577E-02, - 3.583E-02, 3.634E-02, 3.797E-02, 3.987E-02, 4.445E-02, 4.815E-02}; - - Float_t en[kN] = { - 1.00000E-03, 1.07191E-03, 1.14900E-03, 1.14900E-03, 1.50000E-03, - 2.00000E-03, 3.00000E-03, 4.00000E-03, 4.78220E-03, 4.78220E-03, - 5.00000E-03, 5.10370E-03, 5.10370E-03, 5.27536E-03, 5.45280E-03, - 5.45280E-03, 6.00000E-03, 8.00000E-03, 1.00000E-02, 1.50000E-02, - 2.00000E-02, 3.00000E-02, 3.45614E-02, 3.45614E-02, 4.00000E-02, - 5.00000E-02, 6.00000E-02, 8.00000E-02, 1.00000E-01, 1.50000E-01, - 2.00000E-01, 3.00000E-01, 4.00000E-01, 5.00000E-01, 6.00000E-01, - 8.00000E-01, 1.00000E+00, 1.25000E+00, 1.50000E+00, 2.00000E+00, - 3.00000E+00, 4.00000E+00, 5.00000E+00, 6.00000E+00, 8.00000E+00, - 1.00000E+01, 1.50000E+01, 2.00000E+01}; - - return Interpolate(energyMeV, en, mu, kN); + return Interpolate(energyMeV, fEn_xe, fMu_xe, fKN_xe); } //---------------------------------------------------------------------------- //----- GetMuCO2 ------------------------------------------------------ -Float_t CbmTrdRadiator::GetMuCO2(Float_t energyMeV) { +Float_t CbmTrdRadiator::GetMuCO2(Float_t energyMeV) +{ // // Returns the photon absorbtion cross section for CO2 // - - const Int_t kN = 36; - - Float_t mu[kN] = {0.39383E+04, 0.13166E+04, 0.58750E+03, 0.18240E+03, - 0.77996E+02, 0.40024E+02, 0.23116E+02, 0.96997E+01, - 0.49726E+01, 0.15543E+01, 0.74915E+00, 0.34442E+00, - 0.24440E+00, 0.20589E+00, 0.18632E+00, 0.16578E+00, - 0.15394E+00, 0.13558E+00, 0.12336E+00, 0.10678E+00, - 0.95510E-01, 0.87165E-01, 0.80587E-01, 0.70769E-01, - 0.63626E-01, 0.56894E-01, 0.51782E-01, 0.44499E-01, - 0.35839E-01, 0.30825E-01, 0.27555E-01, 0.25269E-01, - 0.22311E-01, 0.20516E-01, 0.18184E-01, 0.17152E-01}; - - Float_t en[kN] = {0.10000E-02, 0.15000E-02, 0.20000E-02, 0.30000E-02, - 0.40000E-02, 0.50000E-02, 0.60000E-02, 0.80000E-02, - 0.10000E-01, 0.15000E-01, 0.20000E-01, 0.30000E-01, - 0.40000E-01, 0.50000E-01, 0.60000E-01, 0.80000E-01, - 0.10000E+00, 0.15000E+00, 0.20000E+00, 0.30000E+00, - 0.40000E+00, 0.50000E+00, 0.60000E+00, 0.80000E+00, - 0.10000E+01, 0.12500E+01, 0.15000E+01, 0.20000E+01, - 0.30000E+01, 0.40000E+01, 0.50000E+01, 0.60000E+01, - 0.80000E+01, 0.10000E+02, 0.15000E+02, 0.20000E+02}; - - return Interpolate(energyMeV, en, mu, kN); + return Interpolate(energyMeV, fEn_co2, fMu_co2, fKN_co2); } //---------------------------------------------------------------------------- //----- GetMuMy ------------------------------------------------------- -Float_t CbmTrdRadiator::GetMuMy(Float_t energyMeV) { +Float_t CbmTrdRadiator::GetMuMy(Float_t energyMeV) +{ // // Returns the photon absorbtion cross section for mylar // - - const Int_t kN = 36; - - Float_t mu[kN] = { - 2.911E+03, 9.536E+02, 4.206E+02, 1.288E+02, 5.466E+01, 2.792E+01, - 1.608E+01, 6.750E+00, 3.481E+00, 1.132E+00, 5.798E-01, 3.009E-01, - 2.304E-01, 2.020E-01, 1.868E-01, 1.695E-01, 1.586E-01, 1.406E-01, - 1.282E-01, 1.111E-01, 9.947E-02, 9.079E-02, 8.395E-02, 7.372E-02, - 6.628E-02, 5.927E-02, 5.395E-02, 4.630E-02, 3.715E-02, 3.181E-02, - 2.829E-02, 2.582E-02, 2.257E-02, 2.057E-02, 1.789E-02, 1.664E-02}; - - Float_t en[kN] = {1.00000E-03, 1.50000E-03, 2.00000E-03, 3.00000E-03, - 4.00000E-03, 5.00000E-03, 6.00000E-03, 8.00000E-03, - 1.00000E-02, 1.50000E-02, 2.00000E-02, 3.00000E-02, - 4.00000E-02, 5.00000E-02, 6.00000E-02, 8.00000E-02, - 1.00000E-01, 1.50000E-01, 2.00000E-01, 3.00000E-01, - 4.00000E-01, 5.00000E-01, 6.00000E-01, 8.00000E-01, - 1.00000E+00, 1.25000E+00, 1.50000E+00, 2.00000E+00, - 3.00000E+00, 4.00000E+00, 5.00000E+00, 6.00000E+00, - 8.00000E+00, 1.00000E+01, 1.50000E+01, 2.00000E+01}; - - return Interpolate(energyMeV, en, mu, kN); + return Interpolate(energyMeV, fEn_my, fMu_my, fKN_my); } //---------------------------------------------------------------------------- //----- Interpolate ------------------------------------------------------ -Float_t CbmTrdRadiator::Interpolate(Float_t energyMeV, - Float_t* en, - Float_t* mu, - Int_t n) { +Float_t CbmTrdRadiator::Interpolate(Float_t energyMeV, const Float_t* en, const Float_t* mu, Int_t n) +{ // // Interpolates the photon absorbtion cross section // for a given energy <energyMeV>. @@ -1401,22 +1057,22 @@ Float_t CbmTrdRadiator::Interpolate(Float_t energyMeV, Int_t index = 0; Int_t istat = Locate(en, n, energyMeV, index, de); if (istat == 0) { + // Float_t result = (mu[index] - de * (mu[index] - mu[index+1]) + // / (en[index+1] - en[index] )); + // return result; Float_t result = - (mu[index] - - de * (mu[index] - mu[index + 1]) / (en[index + 1] - en[index])); - return result; - } else { + (TMath::Log(mu[index]) - de * (TMath::Log(mu[index]) - TMath::Log(mu[index + 1])) / (en[index + 1] - en[index])); + return TMath::Exp(result); + } + else { return 0.0; } } //---------------------------------------------------------------------------- //----- Locate ------------------------------------------------------------ -Int_t CbmTrdRadiator::Locate(Float_t* xv, - Int_t n, - Float_t xval, - Int_t& kl, - Float_t& dx) { +Int_t CbmTrdRadiator::Locate(const Float_t* xv, Int_t n, Float_t xval, Int_t& kl, Float_t& dx) +{ // // Locates a point (xval) in a 1-dim grid (xv(n)) // @@ -1429,18 +1085,12 @@ Int_t CbmTrdRadiator::Locate(Float_t* xv, kl = 0; while (kh - kl > 1) { - if (xval < xv[km = (kl + kh) / 2]) - kh = km; + if (xval < xv[km = (kl + kh) / 2]) kh = km; else kl = km; } if (xval < xv[kl] || xval > xv[kl + 1] || kl >= n - 1) { - printf("Locate failed xv[%d] %f xval %f xv[%d] %f!!!\n", - kl, - xv[kl], - xval, - kl + 1, - xv[kl + 1]); + printf("Locate failed xv[%d] %f xval %f xv[%d] %f!!!\n", kl, xv[kl], xval, kl + 1, xv[kl + 1]); exit(1); } @@ -1449,6 +1099,43 @@ Int_t CbmTrdRadiator::Locate(Float_t* xv, return 0; } //---------------------------------------------------------------------------- +// Photon absorption cross sections + +// Pokalon N470 +constexpr Float_t CbmTrdRadiator::fEn_pok[46]; +constexpr Float_t CbmTrdRadiator::fMu_pok[46]; + +// Kapton +constexpr Float_t CbmTrdRadiator::fEn_ka[46]; +constexpr Float_t CbmTrdRadiator::fMu_ka[46]; + +// Aluminum +constexpr Float_t CbmTrdRadiator::fEn_al[48]; +constexpr Float_t CbmTrdRadiator::fMu_al[48]; + +// Polypropylene/Polyethylene +constexpr Float_t CbmTrdRadiator::fEn_po[36]; +constexpr Float_t CbmTrdRadiator::fMu_po[36]; + +// Carbon +constexpr Float_t CbmTrdRadiator::fEn_c[46]; +constexpr Float_t CbmTrdRadiator::fMu_c[46]; + +// Air +constexpr Float_t CbmTrdRadiator::fEn_air[38]; +constexpr Float_t CbmTrdRadiator::fMu_air[38]; + +// Xenon +constexpr Float_t CbmTrdRadiator::fEn_xe[48]; +constexpr Float_t CbmTrdRadiator::fMu_xe[48]; + +// CO2 +constexpr Float_t CbmTrdRadiator::fEn_co2[36]; +constexpr Float_t CbmTrdRadiator::fMu_co2[36]; + +// Mylar +constexpr Float_t CbmTrdRadiator::fEn_my[36]; +constexpr Float_t CbmTrdRadiator::fMu_my[36]; ClassImp(CbmTrdRadiator) diff --git a/core/detectors/trd/CbmTrdRadiator.h b/core/detectors/trd/CbmTrdRadiator.h index 08d55b6c2470f5a5c59524cda18375cc9b5e8e8b..b2f1f2c4dffa608c3bb3187d7fe46144174b5d37 100644 --- a/core/detectors/trd/CbmTrdRadiator.h +++ b/core/detectors/trd/CbmTrdRadiator.h @@ -12,6 +12,7 @@ #ifndef CBMTRDRADIATOR_H #define CBMTRDRADIATOR_H +#define NCOMPONENTS 10 #include <Rtypes.h> // for THashConsistencyHolder, ClassDef #include <RtypesCore.h> // for Float_t, Int_t, Bool_t, Double_t @@ -20,44 +21,39 @@ class CbmTrdPoint; class TH1D; - +typedef Float_t (*GetMuPtr)(Float_t); class CbmTrdRadiator { public: + struct CbmTrdEntranceWindow { + CbmTrdEntranceWindow(TString def = ""); + Bool_t Init(TString def); + + Int_t fN; // no of components + TString fMat[NCOMPONENTS]; // components' names + Float_t fThick[NCOMPONENTS]; // thickness [cm] of window components + Float_t fDens[NCOMPONENTS]; // density [g/cm^3] of window components + GetMuPtr GetMu[NCOMPONENTS]; // array of mu=f(E) for each component + }; + /** Default constructor **/ - CbmTrdRadiator(); - - /** Constructor **/ - CbmTrdRadiator(Bool_t SimpleTR, - Int_t Nfoils, - Float_t FoilThick, - Float_t GapThick); - - /** Constructor **/ - CbmTrdRadiator(Bool_t SimpleTR, - Int_t Nfoils, - Float_t FoilThick, - Float_t GapThick, - TString material, - TString prototype); - - /** Constructor 20.02.2013 to include external files with tr-spectra**/ - CbmTrdRadiator(Bool_t SimpleTR, TString prototype); + CbmTrdRadiator(Bool_t SimpleTR = true, Int_t Nfoils = 337, Float_t FoilThick = 0.0012, Float_t GapThick = 0.09, + TString material = "pefoam20", TString window = "Kapton"); + + /** Constructor using predefined radiator prototype **/ + CbmTrdRadiator(Bool_t SimpleTR, TString prototype, TString window = "Kapton"); /* // implemented prototypes are: - "A" ALICE like 7xfiber layer + 2x8mm Rohacell layer - "Bshort" POKALON 24µm, 0.7mm gap, 100 foils - "B" POKALON 24µm, 0.7mm gap, 250 foils + "tdr18" 153x2mm PE foam foil layer (146 layers in box, 7 layers in carbon grid) "B++" POKALON 24µm, 0.7mm gap, 350 foils - "C" PE 15µm, 0.7mm gap, 200 foils - "D" PE 15µm, 0.5mm gap, 100 foils - "E" PE 20µm, 0.5mm gap, 120 foils - "F" PE 20µm, 0.25mm gap, 220 foils - "G30" 30xfiber layer - "H" 125x2mm PE foam foil layer - "H++" 177x2mm PE foam foil layer - "Kshort" POKALON 24µm, 0.7mm gap, 100 foils micro-structured - "K" POKALON 24µm, 0.7mm gap, 250 foils micro-structured "K++" POKALON 24µm, 0.7mm gap, 350 foils micro-structured + "G30" 30xfiber layer + + + For the entrance window the general syntax would be c0;c1;c2;... + where ci are planar components. e.g. Al;C;Air;C;Al is the TRD-2D sandwich + There are two shortcuts for the Munster detector namely: + "Mylar" - simple mylar window + "Kapton" - Al(50nm)kapton(25µm ) */ /** Destructor **/ @@ -66,19 +62,9 @@ public: /** Create the needed histograms **/ void CreateHistograms(); - /** Init function **/ - void Init(Bool_t SimpleTR, - Int_t Nfoils, - Float_t FoilThick, - Float_t GapTick, - TString material); - /** Init function **/ - void Init(Bool_t SimpleTR, Int_t Nfoils, Float_t FoilThick, Float_t GapThick); /** Init function **/ void Init(); - // void SetRadPrototype(TString prototype); - /** Spectra production **/ void ProduceSpectra(); @@ -112,23 +98,24 @@ public: TString fWindowFoil; /** Computation of photon absorption cross sections taken from http://physics.nist.gov/PhysRefData/Xcom/html/xcom1.html**/ - Float_t GetMuAl(Float_t energyMeV); - Float_t GetMuPo(Float_t energyMeV); - Float_t GetMuPok(Float_t energyMeV); - Float_t GetMuKa(Float_t energyMeV); - Float_t GetMuAir(Float_t energyMeV); - Float_t GetMuXe(Float_t energyMeV); - Float_t GetMuCO2(Float_t energyMeV); - Float_t GetMuMy(Float_t energyMeV); + static Float_t GetMuAl(Float_t energyMeV); + static Float_t GetMuPo(Float_t energyMeV); + static Float_t GetMuPok(Float_t energyMeV); + static Float_t GetMuKa(Float_t energyMeV); + static Float_t GetMuAir(Float_t energyMeV); + static Float_t GetMuXe(Float_t energyMeV); + static Float_t GetMuCO2(Float_t energyMeV); + static Float_t GetMuMy(Float_t energyMeV); + static Float_t GetMuC(Float_t energyMeV); /** Calculate the TR for a given momentum **/ Float_t GetTR(TVector3 mom); /** Interpolate between given points of table **/ - Float_t Interpolate(Float_t energyMeV, Float_t* en, Float_t* mu, Int_t n); + static Float_t Interpolate(Float_t energyMeV, const Float_t* en, const Float_t* mu, Int_t n); /** Locate a point in 1-dim grid **/ - Int_t Locate(Float_t* xv, Int_t n, Float_t xval, Int_t& kl, Float_t& dx); + static Int_t Locate(const Float_t* xv, Int_t n, Float_t xval, Int_t& kl, Float_t& dx); /** Setters for private data memebers **/ @@ -136,6 +123,8 @@ public: void SetFoilThick(Float_t t) { fFoilThick = t; } void SetGapThick(Float_t t) { fGapThick = t; } void SetSigma(Int_t SigmaT); + /** define plane-parallel Entrance Window section in [cm] **/ + void SetEWwidths(Int_t n, Float_t* w); Bool_t LatticeHit(const CbmTrdPoint* point); @@ -152,7 +141,7 @@ private: Float_t fFoilThick; // Thickness of the foils (cm) Float_t fGapThick; // Thickness of gaps between the foils (cm) TString fFoilMaterial; - Float_t fGasThick; // Thickness of hte active gas volume. + Float_t fGasThick; // Thickness of the active gas volume. /* Parameters fixed in Init() */ @@ -175,6 +164,7 @@ private: Float_t fWinDens; Float_t fWinThick; + CbmTrdEntranceWindow WINDOW; // generic window structure Float_t fCom1; // first component of the gas Float_t fCom2; // second component of the gas @@ -183,10 +173,8 @@ private: static const Int_t fSpRange = 50; // Maximum (keV) of eloss spectrum Float_t fSpBinWidth; // Bin width=fSpNBins/fSpRange - Float_t* - fSigma; //! [fSpNBins] Array of sigma values for the foil of the radiator - Float_t* - fSigmaWin; //! [fSpNBins] Array of sigma values for the entrance window of detector + Float_t* fSigma; //! [fSpNBins] Array of sigma values for the foil of the radiator + Float_t* fSigmaWin; //! [fSpNBins] Array of sigma values for the entrance window of detector Float_t* fSigmaDet; //! [fSpNBins] Array of sigma values for the active gas TH1D* fSpectrum; //! TR photon energy spectrum @@ -206,7 +194,140 @@ private: TVector3 fMom; // momentum of the electron - ClassDef(CbmTrdRadiator, 1) + /* Photon absorption cross section for pokalon N470 */ + static const Int_t fKN_pok = 46; + static constexpr Float_t fEn_pok[46] = { + 1.000E-03, 1.500E-03, 2.000E-03, 3.000E-03, 4.000E-03, 5.000E-03, 6.000E-03, 8.000E-03, 1.000E-02, 1.500E-02, + 2.000E-02, 3.000E-02, 4.000E-02, 5.000E-02, 6.000E-02, 8.000E-02, 1.000E-01, 1.500E-01, 2.000E-01, 3.000E-01, + 4.000E-01, 5.000E-01, 6.000E-01, 8.000E-01, 1.000E+00, 1.022E+00, 1.250E+00, 1.500E+00, 2.000E+00, 2.044E+00, + 3.000E+00, 4.000E+00, 5.000E+00, 6.000E+00, 7.000E+00, 8.000E+00, 9.000E+00, 1.000E+01, 1.100E+01, 1.200E+01, + 1.300E+01, 1.400E+01, 1.500E+01, 1.600E+01, 1.800E+01, 2.000E+01}; + static constexpr Float_t fMu_pok[46] = { + 2.537E+03, 8.218E+02, 3.599E+02, 1.093E+02, 4.616E+01, 2.351E+01, 1.352E+01, 5.675E+00, 2.938E+00, 9.776E-01, + 5.179E-01, 2.847E-01, 2.249E-01, 2.003E-01, 1.866E-01, 1.705E-01, 1.600E-01, 1.422E-01, 1.297E-01, 1.125E-01, + 1.007E-01, 9.193E-02, 8.501E-02, 7.465E-02, 6.711E-02, 6.642E-02, 6.002E-02, 5.463E-02, 4.686E-02, 4.631E-02, + 3.755E-02, 3.210E-02, 2.851E-02, 2.597E-02, 2.409E-02, 2.263E-02, 2.149E-02, 2.056E-02, 1.979E-02, 1.916E-02, + 1.862E-02, 1.816E-02, 1.777E-02, 1.743E-02, 1.687E-02, 1.644E-02}; + + /* Photon absorption cross section for kapton */ + static const Int_t fKN_ka = 46; + static constexpr Float_t fEn_ka[46] = { + 1.000E-03, 1.500E-03, 2.000E-03, 3.000E-03, 4.000E-03, 5.000E-03, 6.000E-03, 8.000E-03, 1.000E-02, 1.500E-02, + 2.000E-02, 3.000E-02, 4.000E-02, 5.000E-02, 6.000E-02, 8.000E-02, 1.000E-01, 1.500E-01, 2.000E-01, 3.000E-01, + 4.000E-01, 5.000E-01, 6.000E-01, 8.000E-01, 1.000E+00, 1.022E+00, 1.250E+00, 1.500E+00, 2.000E+00, 2.044E+00, + 3.000E+00, 4.000E+00, 5.000E+00, 6.000E+00, 7.000E+00, 8.000E+00, 9.000E+00, 1.000E+01, 1.100E+01, 1.200E+01, + 1.300E+01, 1.400E+01, 1.500E+01, 1.600E+01, 1.800E+01, 2.000E+01}; + static constexpr Float_t fMu_ka[46] = { + 2.731E+03, 8.875E+02, 3.895E+02, 1.185E+02, 5.013E+01, 2.555E+01, 1.470E+01, 6.160E+00, 3.180E+00, 1.043E+00, + 5.415E-01, 2.880E-01, 2.235E-01, 1.973E-01, 1.830E-01, 1.666E-01, 1.560E-01, 1.385E-01, 1.263E-01, 1.095E-01, + 9.799E-02, 8.944E-02, 8.270E-02, 7.262E-02, 6.529E-02, 6.462E-02, 5.839E-02, 5.315E-02, 4.561E-02, 4.507E-02, + 3.659E-02, 3.133E-02, 2.787E-02, 2.543E-02, 2.362E-02, 2.223E-02, 2.114E-02, 2.025E-02, 1.953E-02, 1.893E-02, + 1.842E-02, 1.799E-02, 1.762E-02, 1.730E-02, 1.678E-02, 1.638E-02}; + + /* Photon absorption cross section for aluminum */ + static const Int_t fKN_al = 48; + static constexpr Float_t fEn_al[48] = { + 1.000E-03, 1.500E-03, 1.560E-03, 1.560E-03, 2.000E-03, 3.000E-03, 4.000E-03, 5.000E-03, 6.000E-03, 8.000E-03, + 1.000E-02, 1.500E-02, 2.000E-02, 3.000E-02, 4.000E-02, 5.000E-02, 6.000E-02, 8.000E-02, 1.000E-01, 1.500E-01, + 2.000E-01, 3.000E-01, 4.000E-01, 5.000E-01, 6.000E-01, 8.000E-01, 1.000E+00, 1.022E+00, 1.250E+00, 1.500E+00, + 2.000E+00, 2.044E+00, 3.000E+00, 4.000E+00, 5.000E+00, 6.000E+00, 7.000E+00, 8.000E+00, 9.000E+00, 1.000E+01, + 1.100E+01, 1.200E+01, 1.300E+01, 1.400E+01, 1.500E+01, 1.600E+01, 1.800E+01, 2.000E+01}; + static constexpr Float_t fMu_al[48] = { + 1.185E+03, 4.023E+02, 3.621E+02, 3.957E+03, 2.263E+03, 7.881E+02, 3.605E+02, 1.934E+02, 1.153E+02, 5.032E+01, + 2.621E+01, 7.955E+00, 3.442E+00, 1.128E+00, 5.684E-01, 3.681E-01, 2.778E-01, 2.018E-01, 1.704E-01, 1.378E-01, + 1.223E-01, 1.042E-01, 9.276E-02, 8.445E-02, 7.802E-02, 6.841E-02, 6.146E-02, 6.080E-02, 5.496E-02, 5.006E-02, + 4.324E-02, 4.277E-02, 3.541E-02, 3.106E-02, 2.836E-02, 2.655E-02, 2.529E-02, 2.437E-02, 2.369E-02, 2.318E-02, + 2.279E-02, 2.249E-02, 2.226E-02, 2.208E-02, 2.195E-02, 2.185E-02, 2.173E-02, 2.168E-02}; + + /* Photon absorption cross section for polypropylene/polyethylene */ + static const Int_t fKN_po = 36; + static constexpr Float_t fEn_po[36] = { + 1.000E-03, 1.500E-03, 2.000E-03, 3.000E-03, 4.000E-03, 5.000E-03, 6.000E-03, 8.000E-03, 1.000E-02, + 1.500E-02, 2.000E-02, 3.000E-02, 4.000E-02, 5.000E-02, 6.000E-02, 8.000E-02, 1.000E-01, 1.500E-01, + 2.000E-01, 3.000E-01, 4.000E-01, 5.000E-01, 6.000E-01, 8.000E-01, 1.000E+00, 1.250E+00, 1.500E+00, + 2.000E+00, 3.000E+00, 4.000E+00, 5.000E+00, 6.000E+00, 8.000E+00, 1.000E+01, 1.500E+01, 2.000E+01}; + static constexpr Float_t fMu_po[36] = { + 1.894E+03, 5.999E+02, 2.593E+02, 7.743E+01, 3.242E+01, 1.643E+01, 9.432E+00, 3.975E+00, 2.088E+00, + 7.452E-01, 4.315E-01, 2.706E-01, 2.275E-01, 2.084E-01, 1.970E-01, 1.823E-01, 1.719E-01, 1.534E-01, + 1.402E-01, 1.217E-01, 1.089E-01, 9.947E-02, 9.198E-02, 8.078E-02, 7.262E-02, 6.495E-02, 5.910E-02, + 5.064E-02, 4.045E-02, 3.444E-02, 3.045E-02, 2.760E-02, 2.383E-02, 2.145E-02, 1.819E-02, 1.658E-02}; + + /* Photon absorption cross section for carbon */ + static const Int_t fKN_c = 46; + static constexpr Float_t fEn_c[46] = { + 1.000E-03, 1.500E-03, 2.000E-03, 3.000E-03, 4.000E-03, 5.000E-03, 6.000E-03, 8.000E-03, 1.000E-02, 1.500E-02, + 2.000E-02, 3.000E-02, 4.000E-02, 5.000E-02, 6.000E-02, 8.000E-02, 1.000E-01, 1.500E-01, 2.000E-01, 3.000E-01, + 4.000E-01, 5.000E-01, 6.000E-01, 8.000E-01, 1.000E+00, 1.022E+00, 1.250E+00, 1.500E+00, 2.000E+00, 2.044E+00, + 3.000E+00, 4.000E+00, 5.000E+00, 6.000E+00, 7.000E+00, 8.000E+00, 9.000E+00, 1.000E+01, 1.100E+01, 1.200E+01, + 1.300E+01, 1.400E+01, 1.500E+01, 1.600E+01, 1.800E+01, 2.000E+01}; + static constexpr Float_t fMu_c[46] = { + 2.211e+03, 7.004e+02, 3.026e+02, 9.032e+01, 3.778e+01, 1.912e+01, 1.095e+01, 4.576e+00, 2.373e+00, 8.074e-01, + 4.420e-01, 2.562e-01, 2.076e-01, 1.871e-01, 1.753e-01, 1.610e-01, 1.514e-01, 1.347e-01, 1.229e-01, 1.066e-01, + 9.547e-02, 8.715e-02, 8.058e-02, 7.076e-02, 6.362e-02, 6.296e-02, 5.690e-02, 5.179e-02, 4.443e-02, 4.391e-02, + 3.563e-02, 3.047e-02, 2.708e-02, 2.469e-02, 2.291e-02, 2.154e-02, 2.047e-02, 1.959e-02, 1.888e-02, 1.828e-02, + 1.778e-02, 1.735e-02, 1.698e-02, 1.667e-02, 1.615e-02, 1.575e-02}; + + /* Photon absorption cross section for air */ + static const Int_t fKN_air = 38; + static constexpr Float_t fEn_air[38] = { + 0.10000E-02, 0.15000E-02, 0.20000E-02, 0.30000E-02, 0.32029E-02, 0.32029E-02, 0.40000E-02, 0.50000E-02, + 0.60000E-02, 0.80000E-02, 0.10000E-01, 0.15000E-01, 0.20000E-01, 0.30000E-01, 0.40000E-01, 0.50000E-01, + 0.60000E-01, 0.80000E-01, 0.10000E+00, 0.15000E+00, 0.20000E+00, 0.30000E+00, 0.40000E+00, 0.50000E+00, + 0.60000E+00, 0.80000E+00, 0.10000E+01, 0.12500E+01, 0.15000E+01, 0.20000E+01, 0.30000E+01, 0.40000E+01, + 0.50000E+01, 0.60000E+01, 0.80000E+01, 0.10000E+02, 0.15000E+02, 0.20000E+02}; + static constexpr Float_t fMu_air[38] = { + 0.35854E+04, 0.11841E+04, 0.52458E+03, 0.16143E+03, 0.15722E+03, 0.14250E+03, 0.77538E+02, 0.40099E+02, + 0.23313E+02, 0.98816E+01, 0.51000E+01, 0.16079E+01, 0.77536E+00, 0.35282E+00, 0.24790E+00, 0.20750E+00, + 0.18703E+00, 0.16589E+00, 0.15375E+00, 0.13530E+00, 0.12311E+00, 0.10654E+00, 0.95297E-01, 0.86939E-01, + 0.80390E-01, 0.70596E-01, 0.63452E-01, 0.56754E-01, 0.51644E-01, 0.44382E-01, 0.35733E-01, 0.30721E-01, + 0.27450E-01, 0.25171E-01, 0.22205E-01, 0.20399E-01, 0.18053E-01, 0.18057E-01}; + + /* Photon absorption cross section for xenon */ + static const Int_t fKN_xe = 48; + static constexpr Float_t fEn_xe[48] = { + 1.00000E-03, 1.07191E-03, 1.14900E-03, 1.14900E-03, 1.50000E-03, 2.00000E-03, 3.00000E-03, 4.00000E-03, + 4.78220E-03, 4.78220E-03, 5.00000E-03, 5.10370E-03, 5.10370E-03, 5.27536E-03, 5.45280E-03, 5.45280E-03, + 6.00000E-03, 8.00000E-03, 1.00000E-02, 1.50000E-02, 2.00000E-02, 3.00000E-02, 3.45614E-02, 3.45614E-02, + 4.00000E-02, 5.00000E-02, 6.00000E-02, 8.00000E-02, 1.00000E-01, 1.50000E-01, 2.00000E-01, 3.00000E-01, + 4.00000E-01, 5.00000E-01, 6.00000E-01, 8.00000E-01, 1.00000E+00, 1.25000E+00, 1.50000E+00, 2.00000E+00, + 3.00000E+00, 4.00000E+00, 5.00000E+00, 6.00000E+00, 8.00000E+00, 1.00000E+01, 1.50000E+01, 2.00000E+01}; + static constexpr Float_t fMu_xe[48] = { + 9.413E+03, 8.151E+03, 7.035E+03, 7.338E+03, 4.085E+03, 2.088E+03, 7.780E+02, 3.787E+02, 2.408E+02, 6.941E+02, + 6.392E+02, 6.044E+02, 8.181E+02, 7.579E+02, 6.991E+02, 8.064E+02, 6.376E+02, 3.032E+02, 1.690E+02, 5.743E+01, + 2.652E+01, 8.930E+00, 6.129E+00, 3.316E+01, 2.270E+01, 1.272E+01, 7.825E+00, 3.633E+00, 2.011E+00, 7.202E-01, + 3.760E-01, 1.797E-01, 1.223E-01, 9.699E-02, 8.281E-02, 6.696E-02, 5.785E-02, 5.054E-02, 4.594E-02, 4.078E-02, + 3.681E-02, 3.577E-02, 3.583E-02, 3.634E-02, 3.797E-02, 3.987E-02, 4.445E-02, 4.815E-02}; + + /* Photon absorption cross section for CO2 */ + static const Int_t fKN_co2 = 36; + static constexpr Float_t fEn_co2[36] = {0.10000E-02, 0.15000E-02, 0.20000E-02, 0.30000E-02, 0.40000E-02, 0.50000E-02, + 0.60000E-02, 0.80000E-02, 0.10000E-01, 0.15000E-01, 0.20000E-01, 0.30000E-01, + 0.40000E-01, 0.50000E-01, 0.60000E-01, 0.80000E-01, 0.10000E+00, 0.15000E+00, + 0.20000E+00, 0.30000E+00, 0.40000E+00, 0.50000E+00, 0.60000E+00, 0.80000E+00, + 0.10000E+01, 0.12500E+01, 0.15000E+01, 0.20000E+01, 0.30000E+01, 0.40000E+01, + 0.50000E+01, 0.60000E+01, 0.80000E+01, 0.10000E+02, 0.15000E+02, 0.20000E+02}; + static constexpr Float_t fMu_co2[36] = {0.39383E+04, 0.13166E+04, 0.58750E+03, 0.18240E+03, 0.77996E+02, 0.40024E+02, + 0.23116E+02, 0.96997E+01, 0.49726E+01, 0.15543E+01, 0.74915E+00, 0.34442E+00, + 0.24440E+00, 0.20589E+00, 0.18632E+00, 0.16578E+00, 0.15394E+00, 0.13558E+00, + 0.12336E+00, 0.10678E+00, 0.95510E-01, 0.87165E-01, 0.80587E-01, 0.70769E-01, + 0.63626E-01, 0.56894E-01, 0.51782E-01, 0.44499E-01, 0.35839E-01, 0.30825E-01, + 0.27555E-01, 0.25269E-01, 0.22311E-01, 0.20516E-01, 0.18184E-01, 0.17152E-01}; + + /* Photon absorption cross section for mylar */ + static const Int_t fKN_my = 36; + static constexpr Float_t fEn_my[36] = {1.00000E-03, 1.50000E-03, 2.00000E-03, 3.00000E-03, 4.00000E-03, 5.00000E-03, + 6.00000E-03, 8.00000E-03, 1.00000E-02, 1.50000E-02, 2.00000E-02, 3.00000E-02, + 4.00000E-02, 5.00000E-02, 6.00000E-02, 8.00000E-02, 1.00000E-01, 1.50000E-01, + 2.00000E-01, 3.00000E-01, 4.00000E-01, 5.00000E-01, 6.00000E-01, 8.00000E-01, + 1.00000E+00, 1.25000E+00, 1.50000E+00, 2.00000E+00, 3.00000E+00, 4.00000E+00, + 5.00000E+00, 6.00000E+00, 8.00000E+00, 1.00000E+01, 1.50000E+01, 2.00000E+01}; + static constexpr Float_t fMu_my[36] = { + 2.911E+03, 9.536E+02, 4.206E+02, 1.288E+02, 5.466E+01, 2.792E+01, 1.608E+01, 6.750E+00, 3.481E+00, + 1.132E+00, 5.798E-01, 3.009E-01, 2.304E-01, 2.020E-01, 1.868E-01, 1.695E-01, 1.586E-01, 1.406E-01, + 1.282E-01, 1.111E-01, 9.947E-02, 9.079E-02, 8.395E-02, 7.372E-02, 6.628E-02, 5.927E-02, 5.395E-02, + 4.630E-02, 3.715E-02, 3.181E-02, 2.829E-02, 2.582E-02, 2.257E-02, 2.057E-02, 1.789E-02, 1.664E-02}; + + + ClassDef(CbmTrdRadiator, 3) }; - #endif diff --git a/sim/detectors/trd/CbmTrd.cxx b/sim/detectors/trd/CbmTrd.cxx index 1c85f4788ae9f3b59515ccca0bf3fbe3d7f8a7bc..1b00a4dd837872408e879a08f2da91e92da278e7 100644 --- a/sim/detectors/trd/CbmTrd.cxx +++ b/sim/detectors/trd/CbmTrd.cxx @@ -43,8 +43,9 @@ CbmTrd::CbmTrd() , fPosIndex(0) , fTrdPoints(new TClonesArray("CbmTrdPoint")) , fGeoHandler(new CbmTrdGeoHandler()) - , fUseGlobalPhysicsProcesses(kTRUE) - , fCombiTrans(nullptr) { + , fUseGlobalPhysicsProcesses(kFALSE) + , fCombiTrans(nullptr) +{ fVerboseLevel = 1; } // ---------------------------------------------------------------------------- @@ -63,15 +64,17 @@ CbmTrd::CbmTrd(const char* name, Bool_t active) , fPosIndex(0) , fTrdPoints(new TClonesArray("CbmTrdPoint")) , fGeoHandler(new CbmTrdGeoHandler()) - , fUseGlobalPhysicsProcesses(kTRUE) - , fCombiTrans(nullptr) { + , fUseGlobalPhysicsProcesses(kFALSE) + , fCombiTrans(nullptr) +{ fVerboseLevel = 1; } // ---------------------------------------------------------------------------- // ----- Destructor ------------------------------------------------------- -CbmTrd::~CbmTrd() { +CbmTrd::~CbmTrd() +{ if (fTrdPoints) { fTrdPoints->Delete(); delete fTrdPoints; @@ -82,7 +85,8 @@ CbmTrd::~CbmTrd() { // ----- Initialize ------------------------------------------------------- -void CbmTrd::Initialize() { +void CbmTrd::Initialize() +{ FairDetector::Initialize(); // Initialize the CbmTrdGeoHandler helper class from the @@ -94,7 +98,8 @@ void CbmTrd::Initialize() { // ----- SetSpecialPhysicsCuts -------------------------------------------- -void CbmTrd::SetSpecialPhysicsCuts() { +void CbmTrd::SetSpecialPhysicsCuts() +{ FairRun* fRun = FairRun::Instance(); // Setting the properties for the TRDgas is only tested @@ -135,12 +140,12 @@ void CbmTrd::SetSpecialPhysicsCuts() { trdgas->Print(); LOG(fatal) << "Parameters from Virtual Montecarlo:"; LOG(fatal) << "Name " << name << " Aeff=" << a << " Zeff=" << z; - Fatal("CbmTrd::SetSpecialPhysicsCuts", - "Material parameters different between TVirtualMC and TGeomanager"); + Fatal("CbmTrd::SetSpecialPhysicsCuts", "Material parameters different between TVirtualMC and TGeomanager"); } // Set new properties, physics cuts etc. for the TRDgas if (!fUseGlobalPhysicsProcesses) { + LOG(info) << "Using special parameters for TRDgas"; gMC->Gstpar(matIdVMC, "STRA", 1.0); gMC->Gstpar(matIdVMC, "PAIR", 1.0); gMC->Gstpar(matIdVMC, "COMP", 1.0); @@ -157,6 +162,9 @@ void CbmTrd::SetSpecialPhysicsCuts() { gMC->Gstpar(matIdVMC, "DRAY", 1.0); gMC->Gstpar(matIdVMC, "RAYL", 1.0); } + else + LOG(info) << "!! Using global parameters for TRDgas"; + // cut values gMC->Gstpar(matIdVMC, "CUTELE", 10.e-6); gMC->Gstpar(matIdVMC, "CUTGAM", 10.e-6); @@ -174,7 +182,8 @@ void CbmTrd::SetSpecialPhysicsCuts() { // ----- Public method ProcessHits ---------------------------------------- -Bool_t CbmTrd::ProcessHits(FairVolume*) { +Bool_t CbmTrd::ProcessHits(FairVolume*) +{ // Set parameters at entrance of volume. Reset ELoss. if (gMC->IsTrackEntering()) { fELoss = 0.; @@ -188,8 +197,7 @@ Bool_t CbmTrd::ProcessHits(FairVolume*) { fELoss += gMC->Edep(); // Create CbmTrdPoint at exit of active volume - if (gMC->IsTrackExiting() || gMC->IsTrackStop() - || gMC->IsTrackDisappeared()) { + if (gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()) { gMC->TrackPosition(fPosOut); gMC->TrackMomentum(fMomOut); @@ -201,15 +209,9 @@ Bool_t CbmTrd::ProcessHits(FairVolume*) { Int_t size = fTrdPoints->GetEntriesFast(); new ((*fTrdPoints)[size]) - CbmTrdPoint(trackId, - moduleAddress, - TVector3(fPosIn.X(), fPosIn.Y(), fPosIn.Z()), - TVector3(fMomIn.Px(), fMomIn.Py(), fMomIn.Pz()), - TVector3(fPosOut.X(), fPosOut.Y(), fPosOut.Z()), - TVector3(fMomOut.Px(), fMomOut.Py(), fMomOut.Pz()), - fTime, - fLength, - fELoss); + CbmTrdPoint(trackId, moduleAddress, TVector3(fPosIn.X(), fPosIn.Y(), fPosIn.Z()), + TVector3(fMomIn.Px(), fMomIn.Py(), fMomIn.Pz()), TVector3(fPosOut.X(), fPosOut.Y(), fPosOut.Z()), + TVector3(fMomOut.Px(), fMomOut.Py(), fMomOut.Pz()), fTime, fLength, fELoss); // Increment number of trd points in TParticle CbmStack* stack = (CbmStack*) gMC->GetStack(); @@ -224,7 +226,8 @@ Bool_t CbmTrd::ProcessHits(FairVolume*) { // ----- Public method EndOfEvent ----------------------------------------- -void CbmTrd::EndOfEvent() { +void CbmTrd::EndOfEvent() +{ if (fVerboseLevel) Print(); fTrdPoints->Delete(); fPosIndex = 0; @@ -233,16 +236,14 @@ void CbmTrd::EndOfEvent() { // ----- Public method Register ------------------------------------------- -void CbmTrd::Register() { - FairRootManager::Instance()->Register("TrdPoint", "Trd", fTrdPoints, kTRUE); -} +void CbmTrd::Register() { FairRootManager::Instance()->Register("TrdPoint", "Trd", fTrdPoints, kTRUE); } // ---------------------------------------------------------------------------- // ----- Public method GetCollection -------------------------------------- -TClonesArray* CbmTrd::GetCollection(Int_t iColl) const { - if (iColl == 0) - return fTrdPoints; +TClonesArray* CbmTrd::GetCollection(Int_t iColl) const +{ + if (iColl == 0) return fTrdPoints; else return NULL; } @@ -250,7 +251,8 @@ TClonesArray* CbmTrd::GetCollection(Int_t iColl) const { // ----- Public method Print ---------------------------------------------- -void CbmTrd::Print(Option_t*) const { +void CbmTrd::Print(Option_t*) const +{ Int_t nHits = fTrdPoints->GetEntriesFast(); LOG(info) << fName << ": " << nHits << " points registered in this event."; @@ -264,7 +266,8 @@ void CbmTrd::Print(Option_t*) const { // ----- Public method Reset ---------------------------------------------- -void CbmTrd::Reset() { +void CbmTrd::Reset() +{ fTrdPoints->Delete(); ResetParameters(); } @@ -272,7 +275,8 @@ void CbmTrd::Reset() { // ----- Public method CopyClones ----------------------------------------- -void CbmTrd::CopyClones(TClonesArray* cl1, TClonesArray* cl2, Int_t offset) { +void CbmTrd::CopyClones(TClonesArray* cl1, TClonesArray* cl2, Int_t offset) +{ Int_t nEntries = cl1->GetEntriesFast(); LOG(info) << "CbmTrd: " << nEntries << " entries to add."; TClonesArray& clref = *cl2; @@ -289,34 +293,35 @@ void CbmTrd::CopyClones(TClonesArray* cl1, TClonesArray* cl2, Int_t offset) { // ----- ConstructGeometry -------------------------------------------------- -void CbmTrd::ConstructGeometry() { +void CbmTrd::ConstructGeometry() +{ TString fileName = GetGeometryFileName(); - if (fileName.EndsWith(".root")) { - ConstructRootGeometry(); - } else { - LOG(fatal) << "Geometry format of TRD file " << fileName.Data() - << " not supported."; + if (fileName.EndsWith(".root")) { ConstructRootGeometry(); } + else { + LOG(fatal) << "Geometry format of TRD file " << fileName.Data() << " not supported."; } } // ---------------------------------------------------------------------------- //__________________________________________________________________________ -void CbmTrd::ConstructRootGeometry(TGeoMatrix*) { +void CbmTrd::ConstructRootGeometry(TGeoMatrix*) +{ if (Cbm::GeometryUtils::IsNewGeometryFile(fgeoName)) { LOG(info) << "Importing TRD geometry from ROOT file " << fgeoName.Data(); Cbm::GeometryUtils::ImportRootGeometry(fgeoName, this, fCombiTrans); - } else { + } + else { LOG(info) << "Constructing TRD geometry from ROOT file " << fgeoName.Data(); FairModule::ConstructRootGeometry(); } } // ----- CheckIfSensitive ------------------------------------------------- -Bool_t CbmTrd::CheckIfSensitive(std::string name) { +Bool_t CbmTrd::CheckIfSensitive(std::string name) +{ TString tsname = name; if (tsname.EqualTo("gas")) { - LOG(debug) << "CbmTrd::CheckIfSensitive: Register active volume: " - << tsname; + LOG(debug) << "CbmTrd::CheckIfSensitive: Register active volume: " << tsname; return kTRUE; } return kFALSE; diff --git a/sim/detectors/trd/CbmTrdDigitizer.cxx b/sim/detectors/trd/CbmTrdDigitizer.cxx index 0e3379227f7b6eec4ac1e1a55ff12785e59f7ef5..5525fd0efa6c57f82e648a8cb0d06f093e25e36c 100644 --- a/sim/detectors/trd/CbmTrdDigitizer.cxx +++ b/sim/detectors/trd/CbmTrdDigitizer.cxx @@ -73,7 +73,7 @@ CbmTrdDigitizer::CbmTrdDigitizer(CbmTrdRadiator* radiator) // ,fGeoHandler(new CbmTrdGeoHandler()) , fModuleMap() , fDigiMap() { - if (fRadiator == NULL) fRadiator = new CbmTrdRadiator(kTRUE, "K++"); + if (fRadiator == NULL) fRadiator = new CbmTrdRadiator(kTRUE, "tdr18"); }