Added code to create plots of FCal cels shapes

Former-commit-id: 72c9960365fad4640f9e97c6e3532a54e8ec3501

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");
+
  
  
  

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;
@@ -177,18 +186,22 @@ void run_geo()
   
 		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);
+	  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;
 		
-		Identifier identifier= mcell->identify();
+	  Identifier identifier= mcell->identify();
 		/*eta_index64=eta_index;
 		phi_index64=phi_index;
 		if (imap==2) eta_index64+=100;
@@ -203,9 +216,117 @@ void run_geo()
 	
 	
 	}
-	
+
+  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);