diff --git a/mkcode.cc b/mkcode.cc index 011aef5d13e2c7f42723ffdb2e8d459265dcce75..9bb452de67b2a0e03fa17e80f139c625d08c4c18 100644 --- a/mkcode.cc +++ b/mkcode.cc @@ -33,6 +33,9 @@ using namespace CS; #include "tracking_new.h" #endif +//MC simpletracking +#include "simu.h" + int procev(RecoEvent &e); void endrun(); @@ -249,7 +252,8 @@ int main(int argc, char *argv[]) return 1; } - const int Nevents = 100000000; + //const int Nevents = 100000000; + const int Nevents = 2000; int NevMCRead = 0; // event loop, read event by event @@ -544,6 +548,25 @@ int procev(RecoEvent &e0) if(mymom < MomentumLowerCut) ipass = 0; if(mymom > MomentumUpperCut) ipass = 0; + } + else { + if(MCFileFormatVersion >= 8) { + if(ievproc == 0) { + Reset_Calib(); + Init_Calib_2022B_MM(1664928923); + LoadGeometryFromCalo(); + } + DoDigitization(e); + + man.Reset(); + man.AddMicromega(e.MM1); + man.AddMicromega(e.MM2); + man.AddMicromega(e.MM3); + man.AddMicromega(e.MM4); + + track = simple_tracking_4MM_2016(e); + if(track) mymom = track.momentum; + } } #endif #ifdef TRACKINGTOOLS @@ -976,7 +999,7 @@ int procev(RecoEvent &e0) double mychi2 = 0.; //mychi2 = calcShowerChi2(e, 2, 2, 3, 3); - mychi2 = calcShowerChi2(e, imaxcellx, imaxcelly, 3, SP_ECAL); + //mychi2 = calcShowerChi2(e, imaxcellx, imaxcelly, 3, SP_ECAL); if(ECALChi2UpperCut < 101. && mychi2 > ECALChi2UpperCut) ipass = 0; if(ECALMaxCell > -1) { // check max cell if(e.mc || ECALMaxCell <= 29) { // set needed cell numbers from cuts @@ -1099,10 +1122,9 @@ int procev(RecoEvent &e0) << " Veto = " << veto << std::endl; } - + momentum->Fill(mymom); //track.momentum - ehplot->Fill(ecal, hcal); srdplot->Fill(srd); srdplot0->Fill(e.SRD[0].energy/MeV); diff --git a/modifyMakefile.sh b/modifyMakefile.sh new file mode 100755 index 0000000000000000000000000000000000000000..66c63c75b8b3faa86a9c985fa465f7ea117958f9 --- /dev/null +++ b/modifyMakefile.sh @@ -0,0 +1,3 @@ +#!/bin/bash +#modify in ../Makefile to allow for simu.h usage +sed -i.bak -r 's#na64\-kirsanov\/mkcode\.exe\: na64\-kirsanov\/\*\.inc#na64\-kirsanov\/mkcode\.exe\: na64\-kirsanov\/\*\.inc na64\-kirsanov\/simu\.h#g' ../Makefile diff --git a/simu.h b/simu.h new file mode 100644 index 0000000000000000000000000000000000000000..e687110f59e34c0f1cd2571c5fd2213e02092fb3 --- /dev/null +++ b/simu.h @@ -0,0 +1,506 @@ +#pragma once + +//root +#include "TVector3.h" +#include "TRandom.h" +#include "TMath.h" + +//na64 +#include "p348reco.h" +#include "mm.h" +#include "conddb.h" + +//GLOBALS +const double MCsignalCalibrationX = 0.05; +const double MCsignalCalibrationY = 0.15; +const int MCsizeX = 5; +const int MCsizeY = 12; +const bool simulatenoise = true; +//method to compute charge spread +enum DigiMethod {STD,GAUS}; +const DigiMethod MM_chargespread = GAUS; +//Maximum threshold for MM +const double MM_threshold = 22e-6; //22 keV, taken from https://www.sciencedirect.com/science/article/pii/S1875389212018019 +const double GEM_threshold = 22e-6; // same, waiting input for Michael Hoesgen +//gem resolution in mm +double GEMyres = 0.060; +double GEMxres = 0.060; +//maximal distance before the merging +double MergingDistance = 1.75; //mm +//decide if to replicate fake hits +const bool ReplicateHits = true; + +//geometry definition +void Init_Geometry_MC2017() +{ + Reset_Calib(); + //imported from standard example1 + ECAL_pos.pos = TVector3( 0, 0, 0); + MAGD_pos.pos = TVector3( 0, 0, -13975.1); + MAGU_pos.pos = TVector3( 0, 0, -17975.1); + MAG_field = 1.7; + MM1_pos.pos = TVector3(315.,0.,-19375); + MM2_pos.pos = TVector3(315.,0.,-18575); + MM3_pos.pos = TVector3(0.,0.,-1840); + MM4_pos.pos = TVector3(0.,0.,-640); + //all same angle for simplicity + MM1_pos.Angle = M_PI/4; + MM2_pos.Angle = M_PI/4; + MM3_pos.Angle = M_PI/4; + MM4_pos.Angle = M_PI/4; + //same pedestal calibration as in 2017 + // timestamp 1506895200 is `2017 Oct 2 0:00 CEST` (last midnight before the end of beam time) + Init_Calib_2017_MM(1506895200); + NominalBeamEnergy = 100.; +} + +void Init_Geometry_MC2018() +{ + Reset_Calib(); + //imported from standard example1 + ECAL_pos.pos = TVector3( 0, 0, 0); + MAGU_pos.pos = TVector3( 0, 0, -19530.1); + //MAGD_pos.pos = TVector3( 0, 0, -15530.1); + //adjusted down position + MAGD_pos.pos = TVector3( 0, 0, -16600.1); + MAG_field = 1.85; + //position of MM + MM3_pos.pos = TVector3(0., 0., -17593.7); + MM4_pos.pos = TVector3(0., 0., -17226.7); + MM5_pos.pos = TVector3(-203.33, 0., 287.35); + MM6_pos.pos = TVector3(-206.667,0., 387.35); + //all same angle for simplicity + MM1_pos.Angle = M_PI/4; + MM2_pos.Angle = M_PI/4; + MM2_pos.Angle = M_PI/4; + MM3_pos.Angle = M_PI/4; + MM5_pos.Angle = M_PI/4; + MM6_pos.Angle = M_PI/4; + //same pedestal calibration as in 2018 + // timestamp 1529445600 is `2018 June 20 0:00 CEST` (last midnight before the end of beam time) + Init_Calib_2018_MM(1529445600); + NominalBeamEnergy = 150.; +} + +//placeholder for dummy geometry +void Init_Geometry_MC2018_invis() +{ + Reset_Calib(); + //imported from standard example1 + ECAL_pos.pos = TVector3( 0, 0, 0); + MAGU_pos.pos = TVector3( 0, 0, -19530.1); + //MAGD_pos.pos = TVector3( 0, 0, -15530.1); + //adjusted down position + MAGD_pos.pos = TVector3( 0, 0, -16600.1); + MAG_field = 1.85; + //position of MM + MM1_pos.pos = TVector3(0., 0., -16930); + MM2_pos.pos = TVector3(0., 0., -16130); + MM3_pos.pos = TVector3(-319., 0., 854); + MM4_pos.pos = TVector3(-329., 0., 1048); + MM5_pos.pos = TVector3(-381, 0., 3309); + MM6_pos.pos = TVector3(-397.667,0., 3521); + //all same angle for simplicity + MM1_pos.Angle = M_PI/4; + MM2_pos.Angle = M_PI/4; + MM2_pos.Angle = M_PI/4; + MM3_pos.Angle = M_PI/4; + MM5_pos.Angle = M_PI/4; + MM6_pos.Angle = M_PI/4; + //same pedestal calibration as in 2018 + // timestamp 1529445600 is `2018 June 20 0:00 CEST` (last midnight before the end of beam time) + Init_Calib_2018_MM(1529445600); + NominalBeamEnergy = 150.; +} + +void Init_Geometry_MC2022B() +{ + //from conddb.h + ECAL_pos.pos = TVector3(-379.9, 8.5, 24887.1); + MAGU_pos.pos = TVector3( 0, 0, 6350.+506.5); + MAGD_pos.pos = TVector3( 0, 0, 6350.+506.5+4000.); + MAG_field = 1.99; + MM1_pos.pos = TVector3( 1.2, -1.3, 4377.4); + MM2_pos.pos = TVector3( -13.0, -2.5, 5279.9); + MM3_pos.pos = TVector3( -313.0, -3.4, 23181.7); + MM4_pos.pos = TVector3( -364.9, -5.4, 24395.6); + MM1_pos.Angle = M_PI/4; + MM2_pos.Angle = M_PI/4; + MM3_pos.Angle = M_PI/4; + MM4_pos.Angle = M_PI/4; + NominalBeamEnergy = 100.; +} + +//collection of function for Micromegas response simulation + +void MCRealHit(vector<MChit>& mchit,const double& EnergyThreshold) +{ + sort(mchit.begin(),mchit.end(),MChitByID); + vector<MChit> new_hit; + for(unsigned int i(0);i<mchit.size();) + { + int id = mchit[i].id; + double energy = mchit[i].energy; + int particle = mchit[i].particle; + double deposit(0),xpos(0),ypos(0),zpos(0); + while(mchit[i].id == id) + { + deposit += mchit[i].deposit; + xpos += mchit[i].deposit*mchit[i].xpos; + ypos += mchit[i].deposit*mchit[i].ypos; + zpos += mchit[i].deposit*mchit[i].zpos; + ++i; + } + //weighted mean of the hit + xpos /= deposit; + ypos /= deposit; + zpos /= deposit; + if(energy > EnergyThreshold) //deposit would be better, but would need Geant4 revision + new_hit.push_back(MChit (xpos,ypos,zpos,deposit,particle,id,energy)); + } + //all computation is over, substitute vector with one ordered by energy + sort(new_hit.begin(),new_hit.end(),MChitByDeposit); + mchit = new_hit; + +}//end MCRealHit + +void SimulateMultiplexing(MicromegaPlane& plane,unsigned int strip,const double& signal) +{ + //const unsigned int detectedsignal=round(gRandom->Poisson(signal)); + const unsigned int detectedsignal=round(signal); + vector<unsigned int>& PlaneStrips = plane.strips; + const int reverse_index = plane.calib->MM_reverse_map.at(strip); + const vector<int> strips = plane.calib->MM_map.at(reverse_index); + for(unsigned int j(0);j<strips.size();++j) + { + //limit ADC count to 4095 + if(PlaneStrips[strips[j]] + detectedsignal < 4095) + { + PlaneStrips[strips[j]] += detectedsignal; + } + else + { + PlaneStrips[strips[j]] = 4095; + } + } +} ///end SimulateMultiplexing + +//method of computation of strip deposit in strip as function of distance from hit point +double ComputeCharge(const double& signal,const double& dist2center,const unsigned int ClusSize) +{ + if( MM_chargespread == STD ) + { + return signal / pow(2,dist2center+1); + } + else if( MM_chargespread == GAUS ) + { + const double chargespread = ClusSize/4.; + return signal*TMath::Gaus(dist2center,0,chargespread); + } + else + { + cout << "WARNING: NO VALID METHOD FOR CHARGE SPREAD OVER STRIPS SELECTED" << endl; + exit(-1); + } +} + + +void SimulateHit(MicromegaPlane& plane,const double center,const unsigned int meansize, const double& signal) +{ + //define size of the cluster + unsigned int size(0); + while (size < 2) + { + size = gRandom->Poisson(meansize); + } + //in case cluster is asymmetric, choose a direction + int direction(1); + if(size%2 == 0) + { + direction = 1-2*gRandom->Integer(2); + } + //substract center + int strip = round(center); + if((strip < 0) || (((unsigned int)strip + 1) > plane.strips.size())) + return; + double dist2center = fabs(center - strip); + double stripsignal = ComputeCharge(signal,dist2center,size); + SimulateMultiplexing(plane,strip,stripsignal); + unsigned int counter = 1; + while(counter < size) + { + int distance = counter/2 +1; + strip = round(center+direction*distance); + dist2center = fabs(center - strip); + stripsignal = ComputeCharge(signal,dist2center,size); + if((strip < 0) || (((unsigned int)strip + 1) > plane.strips.size())){++counter;continue;} + //multiplex the strip + SimulateMultiplexing(plane,strip,stripsignal); + ++counter; + strip = round(center-direction*distance); + dist2center = fabs(center - strip); + stripsignal = ComputeCharge(signal,dist2center,size); + if((strip < 0) || (((unsigned int)strip + 1) > plane.strips.size()) || counter >= (size)){++counter;continue;} + ++counter; + SimulateMultiplexing(plane,strip,stripsignal); + } +} + + + +void MCDoHitAccumulation(Micromega& MM,vector<MChit> MMtruth, const int verbose = 0) +{ + //compute strip hit + const double angle = MM.geom->Angle; + const double Kx = cos(angle); + const double Ky = sin(angle); + const double xcenter = MM.calib->xcalib.MM_Nstrips/2.0 - 0.5; + const double ycenter = MM.calib->ycalib.MM_Nstrips/2.0 - 0.5; + for(unsigned int i(0);i<MMtruth.size();++i) + { + double xpos = MMtruth[i].xpos - MM.geom->pos.X(); + double ypos = MMtruth[i].ypos - MM.geom->pos.Y(); + const double signal = MMtruth[i].deposit; + //add resolution and convert to strip size + static const double StripRes = 1/sqrt(12); + xpos /= MM.calib->xcalib.MM_strip_step; + ypos /= MM.calib->ycalib.MM_strip_step; + //xpos = gRandom->Gaus(xpos,StripRes); + //ypos = gRandom->Gaus(ypos,StripRes); + /*find best strips for the hit and signal + invert the function in MMCalculatePos to account 45 degree conversion*/ + const double xstrip = (-Kx*xpos - Ky*ypos + xcenter); + const double ystrip = (Ky*xpos - Kx*ypos + ycenter); + const double signalX = signal*MCsignalCalibrationX; + const double signalY = signal*MCsignalCalibrationY; + if(verbose > 0) + cout << " hit with xcenter: " << xcenter << " and ycenter: " << ycenter << " and signalX: " << signalX << " and signalY: " << signalY << endl; + if(xstrip >= 0 && xstrip < MM.calib->xcalib.MM_Nstrips) + { + SimulateHit(MM.xplane,xstrip,MCsizeX,signalX); + } + if(ystrip >= 0 && ystrip < MM.calib->ycalib.MM_Nstrips) + { + SimulateHit(MM.yplane,ystrip,MCsizeY,signalY); + } + } +} //end MCDoHitAccumulation + +//function to simulate noise +void SimulatePlaneNoise(MicromegaPlane& plane,const vector<double>& sigma) +{ + //loop over channels + vector<unsigned int>& PlaneStrips = plane.strips; + for(unsigned int i(0);i<plane.calib->MM_Nchannels;++i) + { + const double noise = gRandom->Gaus(0,sigma[i]); + const vector<int> strips = plane.calib->MM_map.at(i); + const double currentsignal = PlaneStrips[strips[0]]; + if(currentsignal + noise < 2*sigma[i]) //sigma used in data + { + continue; + } + else + { + for(unsigned int j(0);j<strips.size();++j) + { + PlaneStrips[strips[j]] += noise; + } + } + } +} + +void SimulateNoise(Micromega& mm) +{ + SimulatePlaneNoise(mm.xplane,mm.calib->xcalib.sigma); + SimulatePlaneNoise(mm.yplane,mm.calib->ycalib.sigma); +} + +void DoSmearing(vector<MChit>& mmhits,const double& xres,const double& yres) +{ + for(unsigned int i = 0; i < mmhits.size(); i++) + { + mmhits[i].xpos = gRandom->Gaus(mmhits[i].xpos,xres); + mmhits[i].ypos = gRandom->Gaus(mmhits[i].ypos,yres); + } +} + +vector<MChit> DoMerging(vector<MChit> mmhits,const double maxD = MergingDistance) +{ + vector<MChit> mmhitsmerged; + for(auto p1 = mmhits.begin(); p1 != mmhits.end(); ++p1) + { + MChit newhit = *p1; + for(auto p2 = p1+1; p2 != mmhits.end(); ) + { + if(abs(newhit.xpos-(*p2).xpos) < maxD && abs(newhit.ypos-(*p2).ypos) < maxD) + { + //DISTANCE IS TOO SHORT, DO MERGING + //TODO: half of distance is taken, use weighted mean instead? + newhit.xpos = 0.5*(newhit.xpos + (*p2).xpos); + newhit.ypos = 0.5*(newhit.ypos + (*p2).ypos); + newhit.zpos = 0.5*(newhit.zpos + (*p2).zpos); + newhit.deposit = newhit.deposit + (*p2).deposit; + //TODO: as convention the one with largest energy dominates, should it be deposit instead? + newhit.energy = newhit.energy > (*p2).energy ? newhit.energy : (*p2).energy; + newhit.id = newhit.energy > (*p2).energy ? newhit.id : (*p2).id; + //TODO: Dark photon always take precedence + if(abs(newhit.particle) == 100) + newhit.particle = newhit.particle; + else if(abs((*p2).particle) == 100) + newhit.particle = (*p2).particle; + else + newhit.particle = newhit.energy > (*p2).energy ? newhit.particle : (*p2).particle; + //erase merged hit + p2 = mmhits.erase(p2); + }// merging condition + else + { + ++p2; + } + } + mmhitsmerged.push_back(newhit); + } + return mmhitsmerged; +} // end do merging function + +//CREATE SAME REDUNDANCY OBSERVED IN THE DATA DUE TO THE FAKE HIT MATCHING BETWEEN PLANES +vector<MChit> CreateFakeHits(const vector<MChit>& mmhits) +{ + vector<MChit> mmhitsfake = mmhits; + for(auto p1 = mmhits.begin(); p1 != mmhits.end(); ++p1) + { + MChit newhit; + for(auto p2 = mmhits.begin(); p2 != mmhits.end(); ++p2) + { + //skip true combination + if(p1 == p2)continue; + //DISTANCE IS TOO SHORT, DO MERGING + newhit.xpos = (*p1).xpos; + newhit.ypos = (*p2).ypos; + newhit.zpos = (*p2).zpos; + newhit.deposit = 0.5*((*p1).deposit + (*p2).deposit); + //TODO: as convention the one with largest energy dominates, should it be deposit instead? + newhit.energy = (*p1).energy > (*p2).energy ? (*p1).energy : (*p2).energy; + //put Geant ID at -1 to see that it is a fake hits + newhit.id = -1; + //TODO: Dark photon always take precedence + if(abs((*p1).particle) == 100) + newhit.particle = (*p1).particle; + else if(abs((*p2).particle) == 100) + newhit.particle = (*p2).particle; + else + newhit.particle = (*p1).energy > (*p2).energy ? (*p1).particle : (*p2).particle; + mmhitsfake.push_back(newhit); + /* cout<< "new hit created" << endl; */ + /* cout << newhit; */ + } + } + return mmhitsfake; +} // end do merging function + + +void DoMCReco (Micromega& mm,vector<MChit>& hits,const MM_calib& calib,const DetGeo& geom) +{ + //eliminate redundant hits + MCRealHit(hits,MM_threshold); + //assign calib and geom + mm.Init(calib, geom); + //simulate multiplexing and real condition + //MCDoHitAccumulation(mm,hits,1); + MCDoHitAccumulation(mm,hits); + //simulate noise + //if(simulatenoise) + //SimulateNoise(mm); + //do the rest of the reconstruction as for data + mm.DoMMReco(); +} + +void DoMCRecoGEM(vector<MChit>& hits) +{ + //take only physical hit + + MCRealHit(hits,GEM_threshold); + + //smear the hits + DoSmearing(hits,GEMxres,GEMyres); + + //MERGED HITS + hits = DoMerging(hits); + + //create fake hits + if(ReplicateHits) + { + hits = CreateFakeHits(hits); + } +} + +//DIGITALIZE MC INFORMATION + +void DoDigitization(RecoEvent& e) +{ + + //do reconstruction in mc for Micromegas + DoMCReco(e.MM1,e.mc->MM1truth,MM1_calib,MM1_pos); + DoMCReco(e.MM2,e.mc->MM2truth,MM2_calib,MM2_pos); + DoMCReco(e.MM3,e.mc->MM3truth,MM3_calib,MM3_pos); + DoMCReco(e.MM4,e.mc->MM4truth,MM4_calib,MM4_pos); + DoMCReco(e.MM5,e.mc->MM5truth,MM1_calib,MM5_pos); + DoMCReco(e.MM6,e.mc->MM6truth,MM2_calib,MM6_pos); + DoMCReco(e.MM7,e.mc->MM7truth,MM3_calib,MM7_pos); + DoMCReco(e.MM8,e.mc->MM8truth,MM4_calib,MM8_pos); + //merge gem hits + DoMCRecoGEM(e.mc->GEM1truth); + DoMCRecoGEM(e.mc->GEM2truth); + DoMCRecoGEM(e.mc->GEM3truth); + DoMCRecoGEM(e.mc->GEM4truth); + +} + +void LoadGeometryFromCalo() +{ + // loop over MC file metadata + for (auto p : MCRunInfo) + { + const string key = p.first; + const string value = p.second; + + if( key == "ECAL:ECALStart") ECAL_pos.pos.SetZ(mcrunoption<double>(key, 0.)); + else if( key == "ECAL:ECALX") ECAL_pos.pos.SetX(mcrunoption<double>(key, 0.)); + else if( key == "General:BeamMomentum") NominalBeamEnergy = mcrunoption<double>(key, 100.); + else if( key == "MM:MM1X") MM1_pos.pos.SetX(mcrunoption<double>(key, 0.)); + else if( key == "MM:MM1Y") MM1_pos.pos.SetY(mcrunoption<double>(key, 0.)); + else if( key == "MM:MM1Z") MM1_pos.pos.SetZ(mcrunoption<double>(key, 0.)); + else if( key == "MM:MM2X") MM2_pos.pos.SetX(mcrunoption<double>(key, 0.)); + else if( key == "MM:MM2Y") MM2_pos.pos.SetY(mcrunoption<double>(key, 0.)); + else if( key == "MM:MM2Z") MM2_pos.pos.SetZ(mcrunoption<double>(key, 0.)); + else if( key == "MM:MM3X") MM3_pos.pos.SetX(mcrunoption<double>(key, 0.)); + else if( key == "MM:MM3Y") MM3_pos.pos.SetY(mcrunoption<double>(key, 0.)); + else if( key == "MM:MM3Z") MM3_pos.pos.SetZ(mcrunoption<double>(key, 0.)); + else if( key == "MM:MM4X") MM4_pos.pos.SetX(mcrunoption<double>(key, 0.)); + else if( key == "MM:MM4Y") MM4_pos.pos.SetY(mcrunoption<double>(key, 0.)); + else if( key == "MM:MM4Z") MM4_pos.pos.SetZ(mcrunoption<double>(key, 0.)); + else if( key == "MM:MM5X") MM5_pos.pos.SetX(mcrunoption<double>(key, 0.)); + else if( key == "MM:MM5Y") MM5_pos.pos.SetY(mcrunoption<double>(key, 0.)); + else if( key == "MM:MM5Z") MM5_pos.pos.SetZ(mcrunoption<double>(key, 0.)); + else if( key == "MM:MM6X") MM6_pos.pos.SetX(mcrunoption<double>(key, 0.)); + else if( key == "MM:MM6Y") MM6_pos.pos.SetY(mcrunoption<double>(key, 0.)); + else if( key == "MM:MM6Z") MM6_pos.pos.SetZ(mcrunoption<double>(key, 0.)); + else if( key == "MM:MM7X") MM7_pos.pos.SetX(mcrunoption<double>(key, 0.)); + else if( key == "MM:MM7Y") MM7_pos.pos.SetY(mcrunoption<double>(key, 0.)); + else if( key == "MM:MM7Z") MM7_pos.pos.SetZ(mcrunoption<double>(key, 0.)); + else if( key == "MM:MM8X") MM8_pos.pos.SetX(mcrunoption<double>(key, 0.)); + else if( key == "MM:MM8Y") MM8_pos.pos.SetY(mcrunoption<double>(key, 0.)); + else if( key == "MM:MM8Z") MM8_pos.pos.SetZ(mcrunoption<double>(key, 0.)); + //else if( key == "Magnet:Membrane0Z") MAGU_pos.pos.SetZ(mcrunoption<double>(key, 0.)); + //else if( key == "Magnet:Membrane1Z") MAGD_pos.pos.SetZ(mcrunoption<double>(key, 0.)); + else if( key == "Magnet:Magnet2End") MAGU_pos.pos.SetZ(mcrunoption<double>(key, 0.)); + else if( key == "Magnet:Magnet1Start") MAGD_pos.pos.SetZ(mcrunoption<double>(key, 0.)-1000.); + else if( key == "Magnet:MagnetMagY") MAG_field = -1.*mcrunoption<double>(key, 0.); + else if( key == "SRD:ZSandwichBegin")SRD_pos.pos.SetZ(mcrunoption<double>(key, 0.)); + + // WCAL_pos is absent on conddb.h for now + //else if( key == "ECAL:ECAL1Start") WCAL_pos.pos.SetZ(mcrunoption<double>(key, 0.)); //ECAL1 is the name of WCAL in the simulation + //else if( key == "ECAL:ECAL1X") WCAL_pos.pos.SetX(mcrunoption<double>(key, 0.)); //ECAL1 is the name of WCAL in the simulation + } +}