diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/ICaloGeometry.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/ICaloGeometry.h index ac2693f7c0d473396f76791862ffe3e98da2d2a0..235f8006cd4e993bf9d38d93302f1842de66671e 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/ICaloGeometry.h +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/ICaloGeometry.h @@ -18,7 +18,7 @@ public : virtual const CaloDetDescrElement* getDDE(Identifier identify) = 0; virtual const CaloDetDescrElement* getDDE(int sampling,float eta,float phi,float* distance=0,int* steps=0) = 0; - virtual const CaloDetDescrElement* getFCalDDE(int sampling,float x,float y,float z) = 0; + virtual const CaloDetDescrElement* getFCalDDE(int sampling,float x,float y,float z,float* distance=0,int* steps=0) = 0; virtual double deta(int sample,double eta) const = 0; virtual void minmaxeta(int sample,double eta,double& mineta,double& maxeta) const = 0; diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/CaloGeometry.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/CaloGeometry.h index 68d4d7caa48da6269f584e6bee071e0edad322b9..395ef4b984541a6fbc70967a992739443cd0c032 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/CaloGeometry.h +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/CaloGeometry.h @@ -40,7 +40,8 @@ class CaloGeometry : virtual public ICaloGeometry { virtual const CaloDetDescrElement* getDDE(int sampling, Identifier identify); virtual const CaloDetDescrElement* getDDE(int sampling,float eta,float phi,float* distance=0,int* steps=0); - virtual const CaloDetDescrElement* getFCalDDE(int sampling,float x,float y,float z); + virtual const CaloDetDescrElement* getFCalDDE(int sampling,float x,float y,float z,float* distance=0,int* steps=0); + bool getClosestFCalCellIndex(int sampling,float x,float y,int& ieta, int& iphi,int* steps=0); double deta(int sample,double eta) const; void minmaxeta(int sample,double eta,double& mineta,double& maxeta) const; @@ -96,6 +97,8 @@ class CaloGeometry : virtual public ICaloGeometry { bool m_dographs; std::vector< TGraphErrors* > m_graph_layers; FCAL_ChannelMap m_FCal_ChannelMap; // for hit-to-cell assignment in FCal + std::vector<double> m_FCal_rmin,m_FCal_rmax; + /* double m_min_eta_sample[2][MAX_SAMPLING]; //[side][calosample] diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_ntuple_core.py b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_ntuple_core.py index 38d8cf2817f5db4950c24bacea837b9388ae693a..4d8a10fa16eddbe3db8c3a9fb0d8c2e925458d7b 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_ntuple_core.py +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_ntuple_core.py @@ -38,6 +38,43 @@ ISF_HitAnalysis.GeoFileName = 'ISF_Geometry' ############################## +# The following lines added according to instructions on ATLASSIM-3697 + +include("PixelConditionsServices/PixelDCSSvc_jobOptions.py") + +from SiLorentzAngleSvc.LorentzAngleSvcSetup import lorentzAngleSvc + +from SiPropertiesSvc.SiPropertiesSvcConf import SiPropertiesSvc + +pixelSiPropertiesSvc = SiPropertiesSvc(name = "PixelSiPropertiesSvc",DetectorName="Pixel",SiConditionsServices=lorentzAngleSvc.pixelSiliconConditionsSvc) + +ServiceMgr += pixelSiPropertiesSvc + + + +from IOVDbSvc.CondDB import conddb + +conddb.addFolder("DCS_OFL","/SCT/DCS/CHANSTAT") + +conddb.addFolder("DCS_OFL","/SCT/DCS/MODTEMP") + +conddb.addFolder("DCS_OFL","/SCT/DCS/HV") + + + +from SCT_ConditionsServices.SCT_ConditionsServicesConf import SCT_DCSConditionsSvc + +InDetSCT_DCSConditionsSvc = SCT_DCSConditionsSvc(name="InDetSCT_DCSConditionsSvc") + +ServiceMgr += InDetSCT_DCSConditionsSvc + +sctSiPropertiesSvc = SiPropertiesSvc(name = "SCT_SiPropertiesSvc", DetectorName="SCT", SiConditionsServices=lorentzAngleSvc.sctSiliconConditionsSvc) + +ServiceMgr += sctSiPropertiesSvc + + +########################################################## + ISF_HitAnalysis.CaloBoundaryR = 1148.0 ISF_HitAnalysis.CaloBoundaryZ = 3550.0 #before: 3475.0 ISF_HitAnalysis.CaloMargin=100 #=10cm diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeometry.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeometry.cxx index d82a5f5684d981b65c8c42df7cf484458a6cb260..0c2c808baf6ce7bea2bb2aaf4abc1a02031887d0 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeometry.cxx +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeometry.cxx @@ -549,28 +549,79 @@ const CaloDetDescrElement* CaloGeometry::getDDE(int sampling,float eta,float phi return bestDDE; } -const CaloDetDescrElement* CaloGeometry::getFCalDDE(int sampling,float x,float y,float z){ - int isam = sampling - 20; - int iphi,ieta; - Long64_t mask1[]{0x34,0x34,0x35}; - Long64_t mask2[]{0x36,0x36,0x37}; - if(!m_FCal_ChannelMap.getTileID(isam, x, y, ieta, iphi)) return nullptr; - //cout << ieta << "" - Long64_t id = (ieta << 5) + 2*iphi; - if(isam==2)id+= (8<<8); - - if(z>0) id+=(mask2[isam-1] << 12); - else id+=(mask1[isam-1] << 12); - - id = id << 44; - Identifier identify((unsigned long long)id); - - const CaloDetDescrElement* foundcell=m_cells[identify]; - - return foundcell; +const CaloDetDescrElement* CaloGeometry::getFCalDDE(int sampling,float x,float y,float z,float* distance,int* steps){ + int isam = sampling - 20; + int iphi(-100000),ieta(-100000); + Long64_t mask1[]{0x34,0x34,0x35}; + Long64_t mask2[]{0x36,0x36,0x37}; + bool found = m_FCal_ChannelMap.getTileID(isam, x, y, ieta, iphi); + if(steps && found) *steps=0; + if(!found) { + //cout << "Warning: Hit is not matched with any FCal cell! Looking for the closest cell" << endl; + found = getClosestFCalCellIndex(sampling, x, y, ieta, iphi,steps); + } + if(!found) { + cout << "Error: Unable to find the closest FCal cell!" << endl; + return nullptr; + } + + + //cout << "CaloGeometry::getFCalDDE: x:" << x << " y: " << y << " ieta: " << ieta << " iphi: " << iphi << " cmap->x(): " << m_FCal_ChannelMap.x(isam,ieta,iphi) << " cmap->y(): " << m_FCal_ChannelMap.y(isam,ieta,iphi) << endl; + Long64_t id = (ieta << 5) + 2*iphi; + if(isam==2)id+= (8<<8); + + + + if(z>0) id+=(mask2[isam-1] << 12); + else id+=(mask1[isam-1] << 12); + + id = id << 44; + Identifier identify((unsigned long long)id); + + const CaloDetDescrElement* foundcell=m_cells[identify]; + if(distance) { + *distance=sqrt(pow(foundcell->x() - x,2) + pow(foundcell->y() - y,2) ); + } + + return foundcell; } +bool CaloGeometry::getClosestFCalCellIndex(int sampling,float x,float y,int& ieta, int& iphi,int* steps){ + + double rmin = m_FCal_rmin[sampling-21]; + double rmax = m_FCal_rmax[sampling-21]; + int isam=sampling-20; + double a=1.; + const double b=0.01; + const int nmax=100; + int i=0; + + const double r = sqrt(x*x +y*y); + if(r==0.) return false; + const double r_inverse=1./r; + + if((r/rmax)>(rmin*r_inverse)){ + x=x*rmax*r_inverse; + y=y*rmax*r_inverse; + while((!m_FCal_ChannelMap.getTileID(isam, a*x, a*y, ieta, iphi)) && i<nmax){ + a-=b; + i++; + } + } + else { + x=x*rmin*r_inverse; + y=y*rmin*r_inverse; + while((!m_FCal_ChannelMap.getTileID(isam, a*x, a*y, ieta, iphi)) && i<nmax){ + a+=b; + i++; + } + + } + if(steps)*steps=i+1; + return i<nmax ? true : false; +} + bool CaloGeometry::PostProcessGeometry() { for(int i=0;i<MAX_SAMPLING;++i) { @@ -879,10 +930,8 @@ void CaloGeometry::LoadFCalGeometryFromFiles(TString filename1,TString filename2 for(int imodule=1;imodule<=3;imodule++){ i=0; - //while(i<50){ while(1){ - //cout << electrodes[imodule-1]->eof() << endl; (*electrodes[imodule-1]) >> tubeName; if(electrodes[imodule-1]->eof())break; (*electrodes[imodule-1]) >> thisTubeId; // ????? @@ -905,52 +954,45 @@ void CaloGeometry::LoadFCalGeometryFromFiles(TString filename1,TString filename2 if (tileStream2) tileStream2 >> a2; if (tileStream3) tileStream3 >> a3; - //unsigned int tileName= (a3 << 16) + a2; stringstream s; m_FCal_ChannelMap.add_tube(tubeNamestring, imodule, thisTubeId, thisTubeI,thisTubeJ, thisTubeX, thisTubeY,seventh_column); - - - - //cout << "FCal electrodes: " << tubeName << " " << second_column << " " << thisTubeI << " " << thisTubeJ << " " << thisTubeX << " " << thisTubeY << " " << seventh_column << " " << eight_column << " " << ninth_column << endl; - //cout << tileStream1.str() << " " << tileStream2.str() << " " << tileStream3.str() << endl; - //cout << a1 << " " << a2 << " " << a3 << " " << tileName << endl; + i++; } } m_FCal_ChannelMap.finish(); // Creates maps - + for(int imodule=1;imodule<=3;imodule++) delete electrodes[imodule-1]; electrodes.clear(); unsigned long long phi_index,eta_index; float x,y,dx,dy; long id; - //auto it =m_cells_in_sampling[20].rbegin(); - //it--; - //unsigned long long identify=it->first; - //for(auto it=m_cells_in_sampling[i].begin(); it!=m_cells_in_sampling[i].end();it++){ - // - //} + long mask1[]{0x34,0x34,0x35}; long mask2[]{0x36,0x36,0x37}; - + + m_FCal_rmin.resize(3,FLT_MAX); + m_FCal_rmax.resize(3,0.); + + for(int imap=1;imap<=3;imap++){ - + int sampling = imap+20; - + if((int)m_cells_in_sampling[sampling].size() != 2*std::distance(m_FCal_ChannelMap.begin(imap), m_FCal_ChannelMap.end(imap))){ cout << "Error: Incompatibility between FCalChannel map and GEO file: Different number of cells in m_cells_in_sampling and FCal_ChannelMap" << endl; cout << "m_cells_in_sampling: " << m_cells_in_sampling[sampling].size() << endl; cout << "FCal_ChannelMap: " << 2*std::distance(m_FCal_ChannelMap.begin(imap), m_FCal_ChannelMap.end(imap)) << endl; } - for(auto it=m_FCal_ChannelMap.begin(imap);it!=m_FCal_ChannelMap.end(imap);it++){ + phi_index = it->first & 0xffff; eta_index = it->first >> 16; x=it->second.x(); @@ -958,6 +1000,11 @@ void CaloGeometry::LoadFCalGeometryFromFiles(TString filename1,TString filename2 m_FCal_ChannelMap.tileSize(imap, eta_index, phi_index,dx,dy); + double r=sqrt(x*x+y*y); + + if(r<m_FCal_rmin[imap-1])m_FCal_rmin[imap-1]=r; + if(r>m_FCal_rmax[imap-1])m_FCal_rmax[imap-1]=r; + id=(mask1[imap-1]<<12) + (eta_index << 5) +2*phi_index; if(imap==2) id+= (8<<8); @@ -987,7 +1034,6 @@ void CaloGeometry::LoadFCalGeometryFromFiles(TString filename1,TString filename2 } - } diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/init_geo.C b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/init_geo.C index 0117b9b66b4ac2f674c2004916447ee045b2d22a..21038c071fdcc7742faf7c0839e05ce528e4877c 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/init_geo.C +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/init_geo.C @@ -33,7 +33,7 @@ void init_geo() cout<<"init geometry done"<<endl; cout << "running run_geo.C" << endl; - gROOT->ProcessLine(".x run_geo.C"); + diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/run_geo.C b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/run_geo.C index 2af9a6dd6ec43403f3fa7e694c334a45ad6066b2..dce677cfa8c844a291ea685fe2cda6fdc83c85ad 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/run_geo.C +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/run_geo.C @@ -11,6 +11,7 @@ void run_geo(); void WriteGeneralInfo(TString /*cut_label*/, TString lumi, float size, float x, float y); void ATLASLabel(Double_t x,Double_t y,const char* text, float tsize, Color_t color=1); void WriteInfo(TString info, float size, float x, float y, int color=1); +void plotFCalCell(CaloGeometryFromFile* geo,int sampling, double x, double y); void run_geo() { @@ -22,7 +23,7 @@ void run_geo() // Warning: cell lookup in the FCal is not working yet! CaloGeometryFromFile* geo=new CaloGeometryFromFile(); geo->SetDoGraphs(1); - geo->LoadGeometryFromFile("/afs/cern.ch/atlas/groups/Simulation/FastCaloSim/ATLAS-GEO-20-00-01.root","ATLAS-GEO-20-00-01"); + geo->LoadGeometryFromFile("/afs/cern.ch/atlas/groups/Simulation/FastCaloSimV2/Geometry-ATLAS-R2-2016-01-00-01.root", "ATLAS-R2-2016-01-00-01"); TString path_to_fcal_geo_files = "/afs/cern.ch/atlas/groups/Simulation/FastCaloSimV2/"; geo->LoadFCalGeometryFromFiles(path_to_fcal_geo_files + "FCal1-electrodes.sorted.HV.09Nov2007.dat", path_to_fcal_geo_files + "FCal2-electrodes.sorted.HV.April2011.dat", path_to_fcal_geo_files + "FCal3-electrodes.sorted.HV.09Nov2007.dat"); //CaloGeometry::m_debug_identity=3179554531063103488; @@ -33,7 +34,7 @@ void run_geo() cell=geo->getDDE(2,0.24,0.24); //This is not working yet for the FCal!!! //cout<<"Found cell id="<<cell->identify()<<" sample="<<cell->getSampling()<<" eta="<<cell->eta()<<" phi="<<cell->phi()<<endl; - unsigned long long cellid64(3179554531063103488); + unsigned long long cellid64(3746994889972252672); Identifier cellid(cellid64); cell=geo->getDDE(cellid); //This is working also for the FCal @@ -157,14 +158,22 @@ void run_geo() } cmap->finish(); // Creates maps - cmap->print_tubemap(1); - cmap->print_tubemap(2); - cmap->print_tubemap(3); + //cmap->print_tubemap(1); + //cmap->print_tubemap(2); + //cmap->print_tubemap(3); int eta_index,phi_index; Long64_t eta_index64,phi_index64; - double x=130; - double y=180; + //double x=423.755; + //double y=123.41; + double x=431.821; + double y=116.694; + //double x=21; + //double y=27; + + + //double x=436.892; + //double y=28.0237; const CaloDetDescrElement* mcell=0; @@ -173,22 +182,33 @@ void run_geo() cout << endl; cout << "Looking for tile corresponding to [x,y] = [" << x << "," << y << "]" << endl; + float* distance = new float(-99.); + int* steps = new int(-99); + for(int imap=1;imap<=3;imap++){ cout << "Looking for tile in module " << imap << endl; if(!cmap->getTileID(imap,x,y,eta_index,phi_index)){ - cout << "Not found" << endl;continue; - } - cout << "Tile found" << endl; - - cout << "Tile Id " << (eta_index << 16) + phi_index << endl; - cout << "eta index: " << eta_index << endl; - cout << "phi index: " << phi_index << endl; - cout << "Tile position: [x,y] = [" << cmap->x(imap,eta_index,phi_index) << "," << cmap->y(imap,eta_index,phi_index) << "]" << endl; - mcell=geo->getFCalDDE(imap+20,x,y,-1); - cout << "Tile position for calogeometry: [x,y] = [" << mcell->x() << "," << mcell->y() << "]" << endl; + cout << "Not found" << endl; + } + else{ + cout << "Tile found" << endl; + cout << "Tile Id " << (eta_index << 16) + phi_index << endl; + cout << "eta index: " << eta_index << endl; + cout << "phi index: " << phi_index << endl; + float dx,dy; + cmap->tileSize(imap, eta_index, phi_index, dx, dy); + cout << "Tile position: [x,y] = [" << cmap->x(imap,eta_index,phi_index) << "," << cmap->y(imap,eta_index,phi_index) << "]" << " [dx,dy] = [" << dx << "," << dy << "] " << endl; + } + mcell=geo->getFCalDDE(imap+20,x,y,+1,distance,steps); + cout << "Original hit position: [x,y] = [" << x << "," << y << "]" << endl; + cout << "Tile position from CaloGeometry: [x,y] = [" << mcell->x() << "," << mcell->y() << "]" << " [dx,dy] = [" << mcell->dx() << "," << mcell->dy() << "] " << " Identifier: " << mcell->identify() << endl; + cout << "Distance: " << *distance << endl; + cout << "Steps: " << *steps << endl; - Identifier identifier= mcell->identify(); + Identifier identifier= mcell->identify(); + + /*eta_index64=eta_index; phi_index64=phi_index; if (imap==2) eta_index64+=100; @@ -203,9 +223,120 @@ void run_geo() } - + + delete distance; + delete steps; + + bool makeFCalCellPlots=false; + + if(makeFCalCellPlots){ + + gSystem->mkdir("FCalCellShapes"); + + for(int sampling=21;sampling<=23;sampling++){ + TString samplingSTR=Form("sampling_%i",sampling); + gSystem->mkdir("FCalCellShapes/"+samplingSTR); + for(auto it=cmap->begin(sampling-20);it!=cmap->end(sampling-20);it++){ + + plotFCalCell(geo,sampling, it->second.x(),it->second.y()); + + + } + } + } -} +} + +void plotFCalCell(CaloGeometryFromFile* geo,int sampling, double x_orig, double y_orig){ + const CaloDetDescrElement* mcell=0; + mcell=geo->getFCalDDE(sampling,x_orig,y_orig,+1); + if(!mcell){ + cout << "Error in plotFCalCell: Cell not found" << endl; + return; + } + + double x=mcell->x(); + double y=mcell->y(); + double dx=mcell->dx(); + double dy=mcell->dy(); + + int N=1000; + + + double yy= y+1.5*dy; + double xx=x-dx; + vector<double> vecx,vecy; + + for(int i=0;i<N+1;i++){ + + + + for (int j=0;j<N+1;j++){ + xx+=(2.*dx)/N; + mcell=geo->getFCalDDE(sampling,xx,yy,+1); + if(!mcell || !TMath::AreEqualRel(mcell->x(),x,1.E-6) || !TMath::AreEqualRel(mcell->y(),y,1.E-6)){ + //cout << " "; + } + else { + //cout << "X"; + vecx.push_back(xx-x); + vecy.push_back(yy-y); + } + } + //cout << endl; + + yy -= 1.5*(2.*dy)/N; + xx = x-dx; + + } + + double xmax= *std::max_element(vecx.begin(),vecx.end()); + double xmin= *std::min_element(vecx.begin(),vecx.end()); + double ymax= *std::max_element(vecy.begin(),vecy.end()); + double ymin= *std::min_element(vecy.begin(),vecy.end()); + + xmax=std::max(xmax,dx); + xmin=std::min(xmin,-dx); + ymax=std::max(ymax,dy); + ymin=std::min(ymin,-dy); + + TGraph* gr = new TGraph (vecx.size(),&vecx[0],&vecy[0]); + + TH2D* hdummy = new TH2D("h","",10,xmin,xmax,10,ymin,ymax); + + TCanvas* c = new TCanvas("Cell_shape","Cell_shape",600,600); + c->cd(); + gr->SetLineColor(1); + //gr->SetTitle(name); + gr->SetMarkerStyle(21); + gr->SetMarkerColor(1); + gr->SetMarkerSize(0.5); + //gr->GetXaxis()->SetRangeUser(-dx,dx); + //gr->GetYaxis()->SetRangeUser(ymin,ymax); + hdummy->GetXaxis()->SetTitle("#deltax [mm]"); + hdummy->GetYaxis()->SetTitle("#deltay [mm]"); + + hdummy->Draw(); + gr->Draw("P"); + + + TString samplingSTR=Form("sampling_%i",sampling); + TString s="FCalCellShapes/" + samplingSTR + "/FCalCellShape"; + s+=Form("_sampling%i_x%4.2f_y%4.2f",sampling,x,y); + + cout << s << endl; + + c->SaveAs(s+".png"); + + delete hdummy; + delete gr; + delete c; + + + +} + + void ATLASLabel(Double_t x,Double_t y,const char* text, float tsize, Color_t color){ TLatex l; //l.SetTextAlign(12);