#include "CARF/G3Event/interface/G3EventProxy.h"

#include <iostream>
#include <fstream>
#include "Tracker/TkLayout/interface/CmsTracker.h"
#include "Tracker/TkLayout/interface/FullTracker.h"
#include "CommonDet/DetLayout/interface/DetLayer.h"
#include "CommonDet/BasicDet/interface/SimDet.h"
#include "CommonDet/BasicDet/interface/DetUnit.h"
#include "CommonDet/BasicDet/interface/DetType.h"
#include "Tracker/SiPixelDet/interface/PixelReadout.h"
#include "Tracker/SiPixelDet/interface/PixelDet.h"
#include "Tracker/SiPixelDet/interface/PixelDetType.h"
#include "Tracker/SiStripDet/interface/SiStripDet.h"
#include "Tracker/SiStripDet/interface/SiStripDetType.h"
#include <Tracker/CommonStripDet/interface/StripReadout.h>


#include "Utilities/Configuration/interface/Architecture.h"

#include "Utilities/Notification/interface/PackageInitializer.h"
#include "Utilities/Notification/interface/Observer.h" 
//#define SIMHITS
#define DIGIS 
// forward declare classes we have only pointers to
class G3EventProxy; 

class CuTkModule  {
      public:
	CuTkModule(int idc, int iring, int ilayer);
	virtual ~CuTkModule();
	float posx, posy, posz;
	float length, width, thickness, widthAtHalfLength;
	int value;
	int count;
	DetUnit* detUnit;
	int getId(){return idModule; }
	bool notInUse(){return notused;}
        void setUsed(){notused=false;}
        int idModule;
        int ring;
        int layer;
        bool notused;	
};
CuTkModule::CuTkModule(int idc, int ring, int layer){
 idModule = idc;
 this->ring=ring;
 this->layer = layer;
 notused=true;
}

CuTkModule::~CuTkModule(){
}
#include <map>

using namespace std;
class CuTkModuleMap {
 public:
  static map<const DetUnit * , CuTkModule *> moduleMap;
 };
map<const DetUnit * , CuTkModule *>
CuTkModuleMap::moduleMap=map<const DetUnit * , CuTkModule *>();

class SvgModuleMap {
 public:
  static map<const int  , CuTkModule *> smoduleMap;
 };
map< const int  , CuTkModule *>

SvgModuleMap::smoduleMap=map<const int  , CuTkModule *>();


class TestMap : private Observer<G3EventProxy *> {
public:
  TestMap();   //!< default constructor
  ~TestMap();  //!< default destructor

  /// this method will do the user analysis 
  void myAnalysis(G3EventProxy * ev);
  void build();
  void print();
  void drawModule(CuTkModule * mod, int key, int layer);

    double phival(double x, double y){
     double phi;
     double phi1=atan(y/x);
     phi = phi1;
     if(y<0. && x>0) phi = phi1+2.*M_PI;
     if(x<0.)phi=phi1+M_PI;
     if(fabs(y)<0.000001 && x>0)phi=0;
     if(fabs(y)<0.000001&&x<0)phi=M_PI;
     if(fabs(x)<0.000001&&y>0)phi=M_PI/2.;
     if(fabs(x)<0.000001&&y<0)phi=3.*M_PI/2.;

     return phi;

    }

  int getlayers(int subdet, int partdet){
   int ncomponent=0;
   if(subdet == 1){ //1=pixel
	  if(partdet == 1 || partdet == 3){ //1-3=encap
	    ncomponent = 2;
	  }
	  else { ncomponent = 3; } //barrel
	}
	if(subdet== 2){ //2=inner silicon
	   if(partdet == 1 || partdet == 3){ //1-3=encap
	    ncomponent = 3;
	   }
	   else { ncomponent = 4; } //barrel
	}
	if(subdet== 3){ //3=outer silicon
	   if(partdet == 1 || partdet == 3){ //1-3=encap
	     ncomponent = 9;
           }
	   else { ncomponent = 6; } //barrel
	}
       return(ncomponent);
  }   
double  xdpixel(double x){
  double res;
  res= ((x-xmin)/(xmax-xmin)*xsize)+ix;
    return res;
}
double  ydpixel(double y){
  double res;
  double y1;
  y1 = (y-ymin)/(ymax-ymin);
  if(nlay>28)   res= 2*ysize - (y1*2*ysize)+iy;
  else res= ysize - (y1*ysize)+iy;
  return res;
}

void defwindow(int nlay){
  if(posrel){ // separated modules
           xmin=-2.;ymin=-2.;xmax=2.;ymax=2.;
           if(nlay >12 && nlay < 17){
             xmin=-.40;xmax=.40;ymin=-.40;ymax=.40;
           }
           if(nlay>28){
             xmin=-0.1;xmax=3.;ymin=-0.1;ymax=8.5;
             if(nlay<32){xmin=-0.3;xmax=1.0;}
             if(nlay>31&&nlay<36){xmax=2.0;}
             if(nlay>35){ymax=8.;}//inner
           }
  }else{ //overlayed modules
    xmin=-1.3;ymin=-1.3;xmax=1.3;ymax=1.3;
    if(nlay >12 && nlay < 17){
      xmin=-.20;xmax=.20;ymin=-.20;ymax=.20;
    }
    if(nlay>28){
      xmin=-1.5;xmax=1.5;ymin=-1.;ymax=28.;
      if(nlay<32){xmin=-0.5;xmax=0.5;}
      if(nlay>31&&nlay<36){xmin=-1.;xmax=1.;}
    }

  }
    if(nlay<15){
      ix=0;iy=(14-nlay)*ysize;}
    if(nlay>14&&nlay<29){
      ix=3*xsize;iy=(nlay-15)*ysize;}
    if(nlay>28){
      if(nlay==29){ix=(int)(1.5*xsize);iy=0;}
      if(nlay>29 && nlay%2==0){int il=(nlay-28)/2;ix=xsize;iy=il*2*ysize;}
      if(nlay>29 && nlay%2!=0){int il=(nlay-28)/2;ix=2*xsize;iy=il*2*ysize;}
    }
}


  int getrings(int subdet, int partdet, int layer){
   int ncomponent=0;
   if(subdet== 1){ //1=pixel
	  if(partdet== 1 || partdet== 3){ //end-cap
	    ncomponent = 7;
	  }
	  else{ncomponent = 8;} //barrel
	}	
	if(subdet== 2){ //inner-silicon
	  if(partdet== 1 || partdet== 3){ //end-cap
	     ncomponent = 3;
	  }
	  else{ncomponent = 12;} //barrel
	}
	if(subdet== 3){ //outer-silicon
	  if(partdet== 1){ //end-cap-z
	    if (layer== 1) ncomponent = 4;
	    if (layer== 2 || layer== 3) ncomponent = 5;
	    if (layer== 4 || layer== 5 || layer== 6) ncomponent = 6;
	    if (layer== 7 || layer== 8 || layer== 9) ncomponent = 7;
	  }
	  if(partdet== 3){ //endcap+z
	    if (layer== 9) ncomponent = 4;
	    if (layer== 8 || layer== 7) ncomponent = 5;
	    if (layer== 6 || layer== 5 || layer== 4) ncomponent = 6;
	    if (layer== 3 || layer== 2 || layer== 1) ncomponent = 7;
	  }
	  if(partdet== 2){ //barrel
	    ncomponent = 12;
	  }
	}
    return(ncomponent);
      }
int getmodules(int subdet, int partdet, int layer, int ring){
   int ncomponent=0;
      int spicchif[] ={24,24,40,56,40,56,80};
      int spicchib[] ={18,30,42,30,38,46,56,42,48,54,60,66,74};
      int numero_layer = 0;

        if(partdet == 2){ //barrel
           numero_layer = layer-1;
           if(subdet== 2){ //inner
              numero_layer = numero_layer+3;
           }
           if(subdet == 3){ //outer
              numero_layer = numero_layer+7;
           }
           ncomponent = spicchib[numero_layer];
        }
        if(partdet!= 2){ //endcap
           if(subdet== 1)ncomponent=24;//pixel
               else
                   ncomponent = spicchif[7-ring];
        }
    return(ncomponent);
      }
  bool isPixel(int layer){
    if(layer==13 || layer==14 || layer==15||layer==16||layer==29||layer==30||layer==31) return true; else return false;
 }
  bool isStereo(int key){
   int layer=key/100000;
   int ring = key - layer*100000;
    ring = ring/1000;
    if(layer==32 || layer==33 || layer==36 || layer==37) return true;
    if(layer<13 || (layer>16&&layer<29))
       if(ring==1 || ring==2 || ring==5)return true;
    return false;
   }
int nlayer(int det,int part,int lay){
	           if(det==3 && part==1) return lay;
		  	if(det==2 && part==1) return lay+9;
		       	if(det==1 && part==1) return lay+12;
			if(det==1 && part==3) return lay+14;
			if(det==2 && part==3) return lay+16;
			if(det==3 && part==3) return lay+19;
			if(det==1 && part==2) return lay+28;
			if(det==2 && part==2) return lay+31;
			if(det==3 && part==2) return lay+35;
		    }
private:
  /// don't change the name "upDate" - this method is mandatory. 
  void upDate(G3EventProxy * ev) {if (ev!=0) myAnalysis(ev);}
  
private:
  int eventsAnalysed; //!< just to count events that have been analysed
  int PUeventsUsed;   //!< just to count events that have been analysed
  int runsAnalysed;   //!< count the runs
  int lastrun;        //!< to remember the last run analysed
  float oldz;
  int ntotRing[41];
  int firstRing[41];
  int nlay;
  double xmin,xmax,ymin,ymax;
  int xsize,ysize,ix,iy;
  bool posrel;
  ofstream * svgfile;
  

};



// -----------------------------------------------------------------------------
// In the constructor we can initialize. It is automatically called when an
// instance of the class is created.
// -----------------------------------------------------------------------------
TestMap::TestMap() {
  cout << "===========================================================" << endl;  
  cout << "=== Start create new TestMap                       ===" << endl;
  // Initialise observer
  init();
  // Initialise counters
  eventsAnalysed = 0;
  PUeventsUsed = 0;
  runsAnalysed = 0;
  lastrun = 0;
  int ntotmod=0;
  xsize=340;ysize=200;
  ix=0;iy=0;

  for (int subdet=1; subdet < 4; subdet++){
    for (int detpart=1; detpart < 4; detpart++){
      // compute number of layers
        int nlayers = getlayers(subdet,detpart); 
         for(int layer=1; layer < nlayers+1; layer++){
           int nrings = getrings(subdet,detpart,layer);
           ntotRing[nlayer(subdet,detpart,layer)-1]=nrings;
           firstRing[nlayer(subdet,detpart,layer)-1]=1;
            if(subdet==3 && detpart!=2)  firstRing[nlayer(subdet,detpart,layer)-1]= 8-nrings;
         for(int ring=1; ring < nrings+1; ring++){
            int idc = ring; 
            if(subdet==2 && detpart!=2)idc = (nrings+1)-ring;
            if(subdet==3 && detpart!=2)idc = 8-ring;
           int nmodules = getmodules(subdet,detpart,layer,ring);
         for(int module=1; module < nmodules+1; module++){
         CuTkModule * smodule = new CuTkModule(module,idc,nlayer(subdet,detpart,layer));
         int key=nlayer(subdet,detpart,layer)*100000+idc*1000+module;
         SvgModuleMap::smoduleMap[key]=smodule;
         ntotmod++;
      }
         if((detpart==2 && subdet !=1 && layer<3) || (detpart!=2 && subdet !=1 && (idc==1||idc==2||idc==5))) for(int module=1; module < nmodules+1; module++){
         CuTkModule * smodule = new CuTkModule(module+100,idc,nlayer(subdet,detpart,layer));
         int key=nlayer(subdet,detpart,layer)*100000+idc*1000+module+100;
         SvgModuleMap::smoduleMap[key]=smodule;
         ntotmod++;
      }
      }
      }
      }
      }
     
  cout << "=== Done create new TestMap with " << ntotmod <<" modules ===" << endl;
  cout << "===========================================================" << endl;
}
 void TestMap::print(){
   svgfile = new ofstream("svg.dat",ios::out);
   for (int layer=1; layer < 42; layer++){
     nlay=layer;
     defwindow(nlay);
     for (int ring=firstRing[layer-1]; ring < ntotRing[layer-1]+firstRing[layer-1];ring++){
   for (int module=1;module<200;module++) { 
         int key=layer*100000+ring*1000+module;
         CuTkModule * mod = SvgModuleMap::smoduleMap[key];
         if(mod !=0 && !mod->notInUse()){
          drawModule(mod,key,layer);

}
    }  
    }  
    }  
}
TestMap::~TestMap() {
   print();
  cout << "===========================================================" << endl;  
  cout << "=== Start delete TestMap                           ===" << endl;
  cout << " Number of events analysed:    " << eventsAnalysed << endl;
  cout << " Number of pileup events used: " << PUeventsUsed << endl;
  cout << " Number of runs analysed:      " << runsAnalysed << endl;
  cout << "=== Done delete TestMap                            ===" << endl;
  cout << "===========================================================" << endl;  
}
 void TestMap::drawModule(CuTkModule * mod, int key,int nlay){
  int x,y;
  double phi,r,dx,dy;
  double xp[4],yp[4],xp1,yp1;
  double vhbot,vhtop,vhapo;
  double rmedio[]={0.041,0.0701,0.0988,0.255,0.340,0.430,0.520,0.610,0.696,0.782,0.868,0.965,1.080};
  double xt1,yt1,xs1=0.,ys1=0.,xt2,yt2,xs2,ys2,pv1,pv2;
  int green = 0;
  posrel=true;
  double xd[4],yd[4];
  int np = 4;
  int numrec=0;
  int numod=0;
   phi = phival(mod->posx,mod->posy);
       r = sqrt(mod->posx*mod->posx+mod->posy*mod->posy);
      vhbot = mod->width;
      vhtop=mod->width;
      vhapo=mod->length;
   if(nlay < 29){ //endcap
        vhbot = mod->widthAtHalfLength/2.-(mod->width/2.-mod->widthAtHalfLength/2.);
        vhtop=mod->width/2.;
        vhapo=mod->length/2.;
        if(nlay >12 && nlay <17){
          if(posrel)r = r+r;
          xp[0]=r-vhtop;yp[0]=-vhapo;
          xp[1]=r+vhtop;yp[1]=-vhapo;
          xp[2]=r+vhtop;yp[2]=vhapo;
          xp[3]=r-vhtop;yp[3]=vhapo;
        }else{
          if(posrel)r = r + r/3.;
          xp[0]=r-vhapo;yp[0]=-vhbot;
          xp[1]=r+vhapo;yp[1]=-vhtop;
          xp[2]=r+vhapo;yp[2]=vhtop;
          xp[3]=r-vhapo;yp[3]=vhbot;
        }
        for(int j=0;j<4;j++){
          xp1 = xp[j]*cos(phi)-yp[j]*sin(phi);
          yp1 = xp[j]*sin(phi)+yp[j]*cos(phi);
          xp[j] = xp1;yp[j]=yp1;
     }
      } else { //barrel
        numod=mod->idModule;if(numod>100)numod=numod-100;
        int vane = mod->ring;
        if(posrel){
          dx = vhapo;
          phi=M_PI;
          xt1=rmedio[nlay-29]; yt1=-vhtop/2.;
          xs1 = xt1*cos(phi)-yt1*sin(phi);
          ys1 = xt1*sin(phi)+yt1*cos(phi);
          xt2=rmedio[nlay-29]; yt2=vhtop/2.;
          xs2 = xt2*cos(phi)-yt2*sin(phi);
          ys2 = xt2*sin(phi)+yt2*cos(phi);
          dy=phival(xs2,ys2)-phival(xs1,ys1);
          xp[0]=vane*(dx+dx/8.);yp[0]=numod*(dy);
          xp[1]=vane*(dx+dx/8.)+dx;yp[1]=numod*(dy);
          xp[2]=vane*(dx+dx/8.)+dx;yp[2]=numod*(dy)+dy;
          xp[3]=vane*(dx+dx/8.);yp[3]=numod*(dy)+dy;
        }else{
          xt1=r; yt1=-vhtop/2.;
          xs1 = xt1*cos(phi)-yt1*sin(phi);
          ys1 = xt1*sin(phi)+yt1*cos(phi);
          xt2=r; yt2=vhtop/2.;
          xs2 = xt2*cos(phi)-yt2*sin(phi);
          ys2 = xt2*sin(phi)+yt2*cos(phi);
          pv1=phival(xs1,ys1);
          pv2=phival(xs2,ys2);
          if(fabs(pv1-pv2)>M_PI && numod==1)pv1=pv1-2.*M_PI;
          if(fabs(pv1-pv2)>M_PI && numod!=1)pv2=pv2+2.*M_PI;
          xp[0]=mod->posz-vhapo/2.;yp[0]=4.2*pv1;
          xp[1]=mod->posz+vhapo/2.;yp[1]=4.2*pv1;
          xp[2]=mod->posz+vhapo/2. ;yp[2]=4.2*pv2;
          xp[3]=mod->posz-vhapo/2.;yp[3]=4.2*pv2;
        }
        }
 if(isStereo(key))
        {
      np = 3;
      if(mod->idModule>100 ){for(int j=0;j<3;j++){
        xd[j]=xdpixel(xp[j]);yd[j]=ydpixel(yp[j]);
      }
      }else {
        xd[0]=xdpixel(xp[2]);yd[0]=ydpixel(yp[2]);
        xd[1]=xdpixel(xp[3]);yd[1]=ydpixel(yp[3]);
        xd[2]=xdpixel(xp[0]);yd[2]=ydpixel(yp[0]);
      }
        } else {
          for(int j=0;j<4;j++){
        xd[j]=xdpixel(xp[j]);yd[j]=ydpixel(yp[j]);
          }
        }
 int average=0;
 green = 0;
 if(mod->count>0){ average= mod->value/mod->count;
 green = mod->value/mod->count; if (green > 255) green=255;}
if(mod->count > 0 && average>0)
 *svgfile <<"<polygon value=\""<<average<<"\" id=\""<<key<<"\" onmouseover=\"showData(evt);\" onmouseout=\"showData(evt);\" POS=\"Module "<<mod->idModule<<" of ring "<<mod->ring<<" in  Layer " << mod->layer<<" \" fill=\"rgb(255,"<<255-green<<",0)\" points=\"";
else
 *svgfile <<"<polygon value=\""<<mod->value<<"\" id=\""<<key<<"\" onmouseover=\"showData(evt);\" onmouseout=\"showData(evt);\" POS=\"Module "<<mod->idModule<<" of ring "<<mod->ring<<" in Layer "<< mod->layer <<" \" fill=\"white\" points=\"";
 for(int k=0;k<np;k++){
        *svgfile << xd[k] << "," << yd[k] << " " ;
      }
   *svgfile <<"\" />" <<endl;

}
void TestMap::build(){
int spicchif[] ={24,24,40,56,40,56,80};
 int spicchib[] ={18,30,42,30,38,46,56,42,48,54,60,66,74};

    const CmsTracker::LayerContainer& theLayerContainer = FullTracker::instance()->forwardLayers();

    CmsTracker::layer_p_iter iLay;
    int det_type=0;
    int nlay=0;
    int nmod1=0;
    int anello = 0;
    int nmod = 0;
    int ntotmod = 0;
    float r;
    int bar_fow = 1;
    float phi,phi1;    
    //float hbotedge,  htopedge, hapothem, hthickness;
    float rmedioS[]={0.27665, 0.3671, 0.4474, 0.5617, 0.6768, 0.8189, 0.9907};    
    float rmedioP[]={0.0643018, 0.0763302, 0.0896792, 0.106088, 0.118808, 0.135627, 0.143702};
for (iLay = theLayerContainer.begin(); iLay != theLayerContainer.end(); iLay++) {
        nlay++;
        Module layer_type = (*iLay)->module();
        const CmsTracker::DetContainer theDetContainer = (*iLay)->detUnits();
        CmsTracker::det_p_iter iDet;
        for (iDet = theDetContainer.begin() ; iDet != theDetContainer.end() ; iDet++) {     
            float posx = (*iDet)->surface().position().x()/ 100.0;  // cm -> m
            float posy = (*iDet)->surface().position().y()/ 100.0;  // cm -> m
            float posz = (*iDet)->surface().position().z()/ 100.0;  // cm -> m
            r=pow(((double)(posx*posx) + posy*posy),0.5);
            phi1=atan(posy/posx);
           phi=phi1;
           if(posy < 0.&&posx>0)phi=phi1+2.*3.1415;
           if(posx < 0.)phi=phi1+3.1415;
           if(fabs(posy)<0.000001&&posx>0)phi=0;
           if(fabs(posy)<0.000001&&posx<0)phi=3.1415;
           if(fabs(posx)<0.000001&&posy>0)phi=3.1415/2.;
           if(fabs(posx)<0.000001&&posy<0)phi=3.1415 + 3.1415/2.;
           det_type=1;
           if (layer_type==silicon) {
	       det_type=2;
	       for (int i=0;i<7; i++){
	            if (fabs(r-rmedioS[i])<0.015){
	                           anello =i+1;
	                           break;
	            }
	       }
	       nmod=(int)((phi/(2.*3.1415))*spicchif[anello-1]+.1)+1;
	   }
           else{
                for (int i=0; i<7; i++){  
                   if(fabs(r -rmedioP[i])<0.0001){
                     anello=i+1;
                     break;
                   } 
                }
                nmod=(int)((phi/(2.*3.1415))*24+.1)+1;    
           }
           if (layer_type==silicon && (*iDet)->type().isStereo()) nmod=nmod+100;
	   ntotmod++;
		bar_fow=1;
		if(nlay>14)bar_fow=3;  
float length = (*iDet)->type ().bounds ().length () / 100.0; // cm -> m
float width = (*iDet)->type ().bounds ().width () / 100.0;   // cm -> m
float thickness = (*iDet)->type ().bounds ().thickness () / 100.0;  // cm -> m
float widthAtHalfLength = (*iDet)->type ().bounds ().widthAtHalfLength () / 100.0;
         int key=nlay*100000+anello*1000+nmod;
	 CuTkModule * mod = SvgModuleMap::smoduleMap[key];
         CuTkModuleMap::moduleMap[(*iDet)]=mod;
       	
if(mod==0)
  cout << "error in module "<<key <<endl; else
  {
   mod->posx = posx;
   mod->posy = posy;
   mod->setUsed();
   mod->value=0;
   mod->count=0;
     mod->posz = posz;
    mod->length = length;
    mod->width = width;
    mod->thickness = thickness;
    mod->widthAtHalfLength = widthAtHalfLength;
   }
	}
    }
const CmsTracker::LayerContainer& theLayerContainer1 = FullTracker::instance()->barrelLayers();
    CmsTracker::layer_p_iter iLay1;
    nmod1=0;
    nmod = 0;
    //float oldz;
    for (iLay1 = theLayerContainer1.begin(); iLay1 != theLayerContainer1.end(); iLay1++) {
	    nlay++;
            anello = -1;
            Module layer_type1 = (*iLay1)->module();
            const CmsTracker::DetContainer theDetContainer1 = (*iLay1)->detUnits();
            CmsTracker::det_p_iter iDet1;
	    for (iDet1 = theDetContainer1.begin() ; iDet1 != theDetContainer1.end() ; iDet1++) {
         	 float posx1 = (*iDet1)->surface().position().x()/ 100.0;  // cm -> m
                 float posy1 = (*iDet1)->surface().position().y()/ 100.0;  // cm -> m
                 float posz1= (*iDet1)->surface().position().z()/ 100.0;  // cm -> m
		 if (anello == -1){
		     oldz = posz1;
		     anello = 1;
                 }
                 else{
	             if((nlay<36 && fabs(oldz - posz1) > 0.001) ||( nlay>35 && fabs(oldz - posz1) > 0.05 )){
	               oldz = posz1;
                       anello++;
                     }
		 }
	         r=pow(((double)(posx1*posx1) + posy1*posy1),0.5);
	         phi1=atan(posy1/posx1);
	         phi=phi1;
	         if(posy1 < 0.&&posx1>0)phi=phi1+2.*3.1415;
	         if(posx1 < 0.)phi=phi1+3.1415;
                 if(fabs(posy1)<0.000001&&posx1>0)phi=0;
                 if(fabs(posy1)<0.000001&&posx1<0)phi=3.1415;
                 if(fabs(posx1)<0.000001&&posy1>0)phi=3.1415/2.;
                 if(fabs(posx1)<0.000001&&posy1<0)phi=3.1415 + 3.1415/2.;
                 nmod=(int)((phi/(2.*3.1415))*spicchib[nlay-29]+.1)+1;
                 det_type=1;
                 if (layer_type1==silicon)  det_type=2;
                 if (layer_type1==silicon && (*iDet1)->type().isStereo()) nmod=nmod+100; 
                 //const vector<float> par = (*iDet1)->type().shape().parameters();	     
                 ntotmod++;
		 bar_fow=2;

float length = (*iDet1)->type ().bounds ().length () / 100.0; // cm -> m
float width = (*iDet1)->type ().bounds ().width () / 100.0;   // cm -> m
float thickness = (*iDet1)->type ().bounds ().thickness () / 100.0;  // cm -> m
float widthAtHalfLength = (*iDet1)->type ().bounds ().widthAtHalfLength () / 100.0;
         int key=nlay*100000+anello*1000+nmod;
	 CuTkModule * mod = SvgModuleMap::smoduleMap[key];
         CuTkModuleMap::moduleMap[(*iDet1)]=mod;
       	
if(mod==0)
  cout << "**************error********************* in module "<<key <<endl; else
  {
   mod->posx = posx1;
   mod->posy = posy1;
   mod->setUsed();
   mod->value=0;
   mod->count=0;
     mod->posz = posz1;
    mod->length = length;
    mod->width = width;
    mod->thickness = thickness;
    mod->widthAtHalfLength = widthAtHalfLength;
   }
}
    }	
}


void TestMap::myAnalysis(G3EventProxy * ev) { 
      if(eventsAnalysed==0)build();
      int simHitsEv = 0;
      int digiEv = 0;
  cout << "===========================================================" << endl;
  cout << "=== Private analysis of " << eventsAnalysed << " event #"<< ev->simSignal()->id().eventInRun() 
       << " in run #" << ev->simSignal()->id().runNumber() ;
  eventsAnalysed++; // some statistics: count events and runs processed
 CmsTracker::LayerIterator ilayer;
  const CmsTracker::LayerContainer& bl =
    FullTracker::instance()->allLayers();
 for ( ilayer = bl.begin(); ilayer != bl.end(); ilayer++) { //loop over layers

    Module theModule = ilayer->module();
    Det::DetUnitContainer theDets = ilayer->detUnits();
    Det::DetUnitContainer::iterator idet;
    for ( idet = theDets.begin(); idet != theDets.end(); idet++) { //loop on detunits
#ifdef SIMHITS
      int simHitsSize = 0;
 const SimDet* simDetP = (*idet)->simDet();  // Get simhits
      if(simDetP!=0  &&  !simDetP->simHits().empty() ){
        simHitsSize = simDetP->simHits().size();
        CuTkModule * mod = CuTkModuleMap::moduleMap[(*idet)];
        if(mod!=0)mod->value=mod->value+simHitsSize;else cout << "**************************error in CuTkModuleMap **************";
        simHitsEv = simHitsEv + simHitsSize;}
#endif
#ifdef DIGIS 
      CuTkModule * mod = CuTkModuleMap::moduleMap[(*idet)];
     int digisSize = 0;
      if(isPixel(mod->layer)){
      PixelDet* pixDet = (PixelDet *) *idet; //OK
      PixelReadout & rd = (PixelReadout &) (pixDet)->specificReadout();
      if( pixDet !=0) {
        digisSize = rd.ndigis();
      }

      // Get the detector module size
      int numColumns = pixDet->specificTopology().ncolumns();
      int numRows = pixDet->specificTopology().nrows();
      if(digisSize !=0){

        PixelReadout::Range digis = rd.digis();
        for(PixelReadout::DigiContainer::const_iterator di=digis.first;
            di!=digis.second;di++) {  // Loop over digis
          digiEv++;
          GlobalPoint p = (*idet)->toGlobal(LocalPoint(0,0,0));

          float posX =p.x();
          float posY =p.y();
          float posZ =p.z();
          int adc = di->adc();
        if(mod!=0){mod->value=mod->value+adc;mod->count=mod->count+1;}
            else cout << "**************************error in CuTkModuleMap **************";

        }//loop over digis
      }//if digis


    } else {
      SiStripDet* stripDet = (SiStripDet *) *idet; //OK
      StripReadout & rd = (StripReadout &) (stripDet)->specificReadout();

      if( stripDet !=0) {
        digisSize = rd.ndigis();
      }

      if(digisSize !=0){

        StripReadout::Range digis = rd.digis();
        for(StripReadout::DigiContainer::const_iterator di=digis.first;
            di!=digis.second;di++) {  // Loop over digis
          digiEv++;
          GlobalPoint p = (*idet)->toGlobal(LocalPoint(0,0,0));

          float posX =p.x();
          float posY =p.y();
          float posZ =p.z();
          int adc = di->adc();
        if(mod!=0){mod->value=mod->value+adc;mod->count=mod->count+1;}
            else cout << "**************************error in CuTkModuleMap **************";

        }//loop over digis
      }//if digis
      }//if isPixel 
#endif


}// end loop on detunits 

}// end loop on layers
       cout << " nDigi " << digiEv << endl;


  if (ev->simSignal()->id().runNumber() != lastrun) {
    lastrun = ev->simSignal()->id().runNumber();
    runsAnalysed++;
  }

  cout << "===========================================================" << endl;
};

static PKBuilder<TestMap> mybuilder("TestMapBuilder");