//
// description
//             MVD main program
// author
//             J.H.Park Yonsei Uni.
// revision    
//             S.S.Ryu Yonsei Uni.   29\AUG\99
//                 translated to fit into PHOOL 
//
#include <iostream>
#include <iomanip>

//For histograms
#include "PhHistogramFactory.hh"

//PHOOL interface
#include "PHNode.h"
#include "PHIODataNode.h"
#include "PHCompositeNode.h"
#include "PHNodeIterator.h"
#include "mMvdReco.hh"
#include "PdbMvdMap.hh"
#include "PdbMvdCalib.hh"

//MVD classes
#include "MvdReco.hh"
#include "MvdClumps.hh"
#include "MvdBbcVertex.hh"
#include "MvdTimer.hh"

#include "MvdCalibrator.hh"
#include "MvdGeoRd.hh"

// STAF input tables
#include "dMvdGeoWrapper.h"
#include "dMvdDCMWrapper.h"
#include "dBbcOutWrapper.h"

//STAF output tables
#include "dMvbRawWrapper.h"
#include "dMvcRawWrapper.h"
#include "dMvdClmpWrapper.h"

#include "dMvdVertexOutWrapper.h"
#include "dMvddNdEtaOutWrapper.h"
#include "dMvdMultOutWrapper.h"

#define MAX_SINGLE_EV_OUT  10  // N dN/deta hists for single events 
#define ETA_MIN_HIST -7.       // min eta in histograms   
#define ETA_MAX_HIST  7.       // max eta in histograms     
#define N_ETA_BINS_IN_HIST 56  // number of eta bins in histograms

typedef PHIODataNode<PHTable> TableNode_t;

mMvdReco::mMvdReco(PHBoolean flag=false):enablehist(flag){}

mMvdReco::~mMvdReco() { }

// Event method.
PHBoolean mMvdReco::event(PHCompositeNode* topNode,
                          MvdDatabaseNew<PdbMvdMap>* mvdmap,
                          MvdDatabaseNew<PdbMvdCalib>* mvdcalib,
                          MvdDatabaseNew<PdbMvdDeadchannel>* mvddeadchannel){

   // STAF table interface
   PHNodeIterator topIter(topNode);
   TableNode_t *dMvdGeoDataNode,
               *dMvdDCMDataNode,
               *dBbcOutDataNode,
               *dMvbRawDataNode,
               *dMvcRawDataNode,
               *dMvdClmpDataNode,
               *dMvdVertexOutDataNode,
               *dMvddNdEtaOutDataNode,
               *dMvdMultOutDataNode;
   TableNode_t *dMvdTrigParDataNode;
   //TableNode_t *headerDataNode;               

   //input tables
   dMvdGeoWrapper *dMvdGeo;
   dMvdDCMWrapper *dMvdDCM;
   dBbcOutWrapper *dBbcOut;
   //headerWrapper *header;
   dMvdTrigParWrapper *dMvdTrigPar;

   //output tables
   dMvbRawWrapper *dMvbRaw;
   dMvcRawWrapper *dMvcRaw;
   dMvdClmpWrapper *dMvdClmp;
   dMvdVertexOutWrapper *dMvdVertexOut;
   dMvddNdEtaOutWrapper *dMvddNdEtaOut;
   dMvdMultOutWrapper *dMvdMultOut;
     
   //keep track of current event 
   static int run=0;
   run++;
   

   cout <<" give idate: ";
   int idate;
//   cin >> idate;
   idate = 20010724;                 // date for geo definition
   int istat = mvdread(idate);       // in file MvdGeoRd.cc
   cout << " mMvdReco::  mvdread returned status "<<istat<<" \n" ;


   //find the input tables; dMvdGeo,dMvdDCM
   dMvdGeoDataNode = static_cast<TableNode_t*>(topIter.findFirst("PHIODataNode",
                                                     "dMvdGeo"));
   if (dMvdGeoDataNode) {
      dMvdGeo = static_cast<dMvdGeoWrapper*>(dMvdGeoDataNode->getData());
   }
   else{
      cout<<"warning; dMvdGeo not found"<<endl;
   }

   dMvdDCMDataNode = static_cast<TableNode_t*>(topIter.findFirst("PHIODataNode",
                                                     "dMvdDCM"));
   if (dMvdDCMDataNode) {
      dMvdDCM = static_cast<dMvdDCMWrapper*>(dMvdDCMDataNode->getData());
      if (dMvdDCM->RowCount()==0){
         cerr<<"dMvdDCM is empty"<<endl;
         return true;
      }
   }
   else{
      cout<<"error; dMvdDCM not found"<<endl;
      return true;
   }

   dBbcOutDataNode = static_cast<TableNode_t*>(topIter.findFirst("PHIODataNode",
                                                     "dBbcOut"));
   if (dBbcOutDataNode) {
      dBbcOut = static_cast<dBbcOutWrapper*>(dBbcOutDataNode->getData());
   }
   else{
      cout<<"warning : dBbcOut not found"<<endl;
   }

   //headerDataNode = static_cast<TableNode_t*>(topIter.findFirst("PHIODataNode",
   //                                                  "header"));
   //if (headerDataNode) {
   //   header = static_cast<headerWrapper*>(headerDataNode->getData());
   //}
   //else{
   //   cout<<"warning : header not found"<<endl;
   //}

   dMvdTrigParDataNode = static_cast<TableNode_t*>(topIter.findFirst("PHIODataNode",
                                                     "dMvdTrigPar"));
   if (dMvdTrigParDataNode) {
      dMvdTrigPar = static_cast<dMvdTrigParWrapper*>(dMvdTrigParDataNode->getData());
   }
   else{
      cout<<"warning : dMvdTrigPar not found"<<endl;
   }


   //find the output tables; dMvbRaw, dMvcRaw, dMvdClmp, 
   //                        dMvdMultOut,dMvdVertexOut,dMvddNdEtaOut

   dMvbRawDataNode = static_cast<TableNode_t*>(topIter.findFirst("PHIODataNode",
                                                     "dMvbRaw"));
   if (dMvbRawDataNode) {
      dMvbRaw = static_cast<dMvbRawWrapper*>(dMvbRawDataNode->getData());
   }
   else{
      cout<<"warning : dMvbRaw not found"<<endl;
   }

   dMvcRawDataNode = static_cast<TableNode_t*>(topIter.findFirst("PHIODataNode",
                                                     "dMvcRaw"));
   if (dMvcRawDataNode) {
      dMvcRaw = static_cast<dMvcRawWrapper*>(dMvcRawDataNode->getData());
   }
   else{
      cout<<"warning : dMvcRaw not found"<<endl;
   }

   dMvdClmpDataNode = static_cast<TableNode_t*>(topIter.findFirst("PHIODataNode",
                                                     "dMvdClmp"));
   if (dMvdClmpDataNode) {
      dMvdClmp = static_cast<dMvdClmpWrapper*>(dMvdClmpDataNode->getData());
   }
   else{
      cout<<"warning : dMvdClmp not found"<<endl;
   }

   dMvdMultOutDataNode = static_cast<TableNode_t*>(topIter.findFirst("PHIODataNode",
                                                     "dMvdMultOut"));
   if (dMvdMultOutDataNode) {
      dMvdMultOut = static_cast<dMvdMultOutWrapper*>(dMvdMultOutDataNode->getData());
   }
   else{
      cout<<"warning : dMvdMultOut not found"<<endl;
   }

   dMvdVertexOutDataNode = static_cast<TableNode_t*>(topIter.findFirst("PHIODataNode",
                                                     "dMvdVertexOut"));
   if (dMvdVertexOutDataNode) {
      dMvdVertexOut = static_cast<dMvdVertexOutWrapper*>(dMvdVertexOutDataNode->getData());
   }
   else{
      cout<<"warning : dMvdVertexOut not found"<<endl;
   }

   dMvddNdEtaOutDataNode = static_cast<TableNode_t*>(topIter.findFirst("PHIODataNode",
                                                     "dMvddNdEtaOut"));
   if (dMvddNdEtaOutDataNode) {
      dMvddNdEtaOut = static_cast<dMvddNdEtaOutWrapper*>(dMvddNdEtaOutDataNode->getData());
   }
   else{
      cout<<"warning : dMvddNdEtaOut not found"<<endl;
   }


   //initialize par
   MvdParameter par;
   par.Init();

   // initialize geometry parameters
   // if dMvdGeo has been filled up correctly, initilize it by copying it,
   // otherwise, use hardwired default values. 
   MvdGeoParameter geopar;
   if (dMvdGeoDataNode && dMvdGeo->RowCount()>0) geopar.Init(dMvdGeo);
   else{
      cout<<"warning; dMvdGeo is empty"<<dMvdGeo->RowCount()<<endl;
      geopar.Init();
   }

   //initialize geometry 
   MvbGeometry bgeometry(par,geopar);
   MvcGeometry cgeometry(par,geopar);
   bgeometry.Init();
   cgeometry.Init();
   MvdBbcVertex bbcvertex;
   MvdClumps mvdclumps(par,geopar);
   MvdReco mvd(par,geopar,mvdmap, mvdcalib,mvddeadchannel);

   if (dBbcOutDataNode && dBbcOut->RowCount()>0) {
       bbcvertex=MvdBbcVertex(0,0,dBbcOut->get_VertexPoint(0),1);
   } else{
       cout<<"warning : no or empty dBbcOut"<<endl;
   }

   mvd.Init(dMvdDCM);
   mvd.Process(mvdclumps,bbcvertex, dMvbRaw, dMvcRaw, dMvdTrigPar);

   //cout<<"dMvbRaw->RowCount()="<<dMvbRaw->RowCount()<<endl;
   //mvd.Show_vertex();
   //cout<<"bbc vertex="<<dBbcOut->get_VertexPoint(0)<<endl;

/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////

   // dump to tables
   mvd.dump2table(dMvdClmp);
   MvdMultiOut  mult_tmp=mvd.BMulti();
   MvdVertex vertex_tmp=mvd.CalcVertex();
   MvdOccupancy occup_tmp=mvd.Occups();
   vector<MvddNdEtaOut>  dndetas_tmp=mvd.DNdEta();

   float totalmult,totalmulterr,softmult;

   float x_vertex,y_vertex,z_vertex;
   float error[3][3];
   float clevel;
   short ecode;
   float softvertex;

   float dndeta,dndetaerr,eta,deta,phi,dphi,softdndeta;

   mult_tmp.Put(totalmult,totalmulterr,softmult);

   for(short k=0; k<1; ++k) {
       (*dMvdMultOut)[k].id=run;
       (*dMvdMultOut)[k].totalmult=totalmult;
       (*dMvdMultOut)[k].totalmulterr=totalmulterr;
       (*dMvdMultOut)[k].softmult=softmult;
   }
   //then don't forget to set the row count!!
   dMvdMultOut->SetRowCount(1);

   x_vertex=vertex_tmp.get_x();
   y_vertex=vertex_tmp.get_y();
   z_vertex=vertex_tmp.get_z();

   for(int i=0; i<3; i++){
      for(int j=0; j<3; j++){
         error[i][j]=vertex_tmp.get_error(i,j);
      }
   }
   clevel=vertex_tmp.get_vertexcl();
   ecode =vertex_tmp.get_errorcode();
   softvertex=vertex_tmp.get_softvertex();

   for(short k=0; k<1; ++k){
      (*dMvdVertexOut)[k].id=run;
      (*dMvdVertexOut)[k].vertexcl=clevel;
      (*dMvdVertexOut)[k].ErrorCode=ecode;
      (*dMvdVertexOut)[k].vertex[0]=x_vertex;
      (*dMvdVertexOut)[k].vertex[1]=y_vertex;
      (*dMvdVertexOut)[k].vertex[2]=z_vertex;
      for (short i=0;i<3;++i){
         for (short j=0;j<3;++j){
            (*dMvdVertexOut)[k].vertexerr[i][j]=error[i][j];
         }
      }
      (*dMvdVertexOut)[k].softvertex=softvertex;
   }
   //don't forget to set the row count correctly!!
   dMvdVertexOut->SetRowCount(1);

   for (unsigned int k=0; k <dndetas_tmp.size(); ++k) {
       dndetas_tmp[k].Put(dndeta,dndetaerr,eta,deta,phi,dphi,softdndeta);
       (*dMvddNdEtaOut)[k].id=run;
       (*dMvddNdEtaOut)[k].dndeta=dndeta;
       (*dMvddNdEtaOut)[k].dndetaerr=dndetaerr;
       (*dMvddNdEtaOut)[k].eta=eta;
       (*dMvddNdEtaOut)[k].deta=deta;
       (*dMvddNdEtaOut)[k].phi=phi;
       (*dMvddNdEtaOut)[k].dphi=dphi;
       (*dMvddNdEtaOut)[k].softdndeta=softdndeta;
   }
   // don't forget to set the row count correctly!!
   dMvddNdEtaOut->SetRowCount(dndetas_tmp.size());


   //histogramming
   if (enablehist){
      //histogram(mvd, header);
      histogram(mvd);
   }

   return true;
}

//void mMvdReco::histogram(const MvdReco& mvd, headerWrapper* header) {
void mMvdReco::histogram(const MvdReco& mvd) {


   MvdParameter       par=mvd.Par();
   MvdGeoParameter geopar=mvd.GeoPar();
   MvbGeometry  bgeometry=mvd.BGeom();
   MvcGeometry  cgeometry=mvd.CGeom();

   float offset,paneloccup,sumkev,summev,xnptls;
   float innerclmpsum;
   float outerclmpsum;

   float zvertex_fit,xvertex_fit,yvertex_fit;
   short errorcode;
   float softvertex;

   float itrker;
   float xntpl[26];
   float occup[4],mip[4];
   float occupancy;
   float xadc;
   float cosring;
   float bimevt;
   float clsum;       // number of strips in a panel 
   int   ncl;         // number of clumps in a panel
   short EndOfPanel;  // used as a switch 
   float diff;
   float tantheta;
   float nhitfactor;
   float zpanelpos;
   float dndeta,dndetaerr,eta,deta,phi,dphi,softdndeta;

   int nstrip=par.nstrip;;
   int NCHDIM=par.nchdim;
   float ZOFF=par.zoff;
   short NPHUSE=par.nphuse;
   int MvdMaxAdc=par.maxadc_ver;

   float xmin,xmax;
   xmin=0.;
   xmax=float(MvdMaxAdc);

   int i,j,k;

   short shell,row,panel,strip,adc;
   short shell_next,row_next,panel_next,strip_next;
   short end,wedge,column;

   PhHistogramFactory *hf = PhHistogramFactory::instance();

   // The 2d declarations here are all really profile histograms originally.
   // A 2d histogram is not really a substitute for the way the profile
   // histograms were being used. I have not even setup the arguements
   // to be correct for declaring a 2d histograms -- they are just left
   // as they were for the profile histograms

   static PhH1 *TrueMult   = hf->CreateH1F("TrueMult",   "True total Mult.",100,  0.,20000.);
   //static PhH1 *TrueZvertex= hf->CreateH1F("TrueZvertex","true z vertex",    50,-50.,50.);
   static PhH1 *dZvertex_pseudo   = hf->CreateH1F("dZvertex_pseudo","dZ vertex pseudo (found-true)",
                                    500,-5.,5.);
   static PhH1 *dZvertex_correlation= hf->CreateH1F("dZvertex_correlation","dZ vertex correlation(found-true)",
                                   500,-5.,5.);

   static PhH1 *herrorcode = hf->CreateH1F("herrorcode","vertexing error code", 7, -5.5, 1.5); 
   static PhH1 *NclumpVsZ1 = hf->CreateH1F("NclumpVsZ1","Inner Shell : NClump vs. z",
                                   NCHDIM,-ZOFF,ZOFF);
   static PhH1 *NclumpVsZ2 = hf->CreateH1F("NclumpVsZ2","Outer Shell : NClump vs. z",
                                   NCHDIM,-ZOFF,ZOFF);
   static PhH1 *Nclumps1 = hf->CreateH1F("Nclumps1","Nclumps tot R1", 50,  0.,10000.);
   static PhH1 *Nclumps2 = hf->CreateH1F("Nclumps2","Nclumps tot R2", 50,  0.,10000.);
   static PhH1 *StripAdc = hf->CreateH1F("StripAdc","strip Adc value",1000, 0., 1000.);
   static PhH1 *ClumpAdc = hf->CreateH1F("ClumpAdc","ClumpAdc",1000, 0., 1000.);
   static PhH1 *PadAdc1  = hf->CreateH1F("PadAdc1","pad Adc value", 1000, 0., 1000.);

   static PhH1 *PadAdcocclt1  = hf->CreateH1F("PadAdcocclt1","pad Adc occ.lt.0.01",
                                   300, 0., 300.);
   static PhH1 *PadAdcocc1to2 =hf->CreateH1F("PadAdcocc1to2","pad Adc 0.01.gt.occ.lt.0.02",
                                   300, 0., 300.);
   static PhH1 *PadAdcocc2to5 =hf->CreateH1F("PadAdcocc2to5","pad Adc 0.02.gt.occ.lt.0.05",
                                   300, 0., 300.);
   static PhH1 *PadAdcocc5to10=hf->CreateH1F("PadAdcocc5to10","pad Adc 0.05.gt.occ.lt.0.10",
                                   300, 0., 300.);
   static PhH1 *PadAdcoccgt10 =hf->CreateH1F("PadAdcoccgt10","pad Adc occ.gt.0.1",
                                   300, 0., 300.);
   static PhH1 *MeVinSi = hf->CreateH1F("MeVinSi","MeV in Si", 100, 0., 2000.);

   static PhH1 *AdcSumR1= hf->CreateH1F("AdcSumR1","Adc sum inner barrel",
                                   100,0.,20000.);
   static PhH1 *AdcSumR2= hf->CreateH1F("AdcSumR2","Adc sum outer barrel",
                                   100,0.,20000.);
   static PhH1 *AdcSumZp= hf->CreateH1F("AdcSumZp","Adc sum +Z pads",100,0.,20000.);
   static PhH1 *AdcSumZm= hf->CreateH1F("AdcSumZm","Adc sum -Z pads",100,0.,20000.);

   static PhH1 *OccR1 = hf->CreateH1F("OccR1","Avg occ. inner barrel",50,0.,1.);
   static PhH1 *OccR2 = hf->CreateH1F("OccR2","Avg occ. outer barrel",50,0.,1.);
   static PhH1 *OccZp = hf->CreateH1F("OccZp","Avg occ. +Z pads",     50,0.,1.);
   static PhH1 *OccZm = hf->CreateH1F("OccZm","Avg occ. -Z pads",     50,0.,1.);
   static PhH1 *bfm   = hf->CreateH1F("bfm","b(fm) all events",10, 0., 12.);

   static PhProfile *ClumpVsZ1= hf->CreateProfile("ClumpVsZ1",
                                       "Inner Shell:Clumpsize vs z-zvtx",
                                        2*NCHDIM, -64.,64.,kErrorOnMean );
   static PhProfile *ClumpVsZ2= hf->CreateProfile("ClumpVsZ2",
                                       "Outer Shell : Clumpsize vs z-zvtx",
                                        2*NCHDIM, -64.,64.,kErrorOnMean);
   static PhProfile *OccVsZR1 = hf->CreateProfile("OccVsZR1","occ. R1",
                                        NCHDIM,-ZOFF,ZOFF,kErrorOnMean);
   static PhProfile *OccVsZR2 = hf->CreateProfile("OccVsZR2","occ. R2",
                                        NCHDIM,-ZOFF,ZOFF,kErrorOnMean);
   static PhProfile *MeVvsZvert= hf->CreateProfile("MeVvsZvert","MeV in Si vs Zvtx",
                                        2*NCHDIM, -64.,64.,kErrorOnMean);
   static PhProfile *dNdetaAvg = hf->CreateProfile("dNdetaAvg","dN/d[c] (event averaged)",
                                         N_ETA_BINS_IN_HIST,ETA_MIN_HIST,
                                         ETA_MAX_HIST,kErrorOnMean);
   static PhProfile *dNdetaRow0= hf->CreateProfile("dNdetaRow0",
                                        "dN/d[c]d[f] (event averaged) row 0",
                                         N_ETA_BINS_IN_HIST,ETA_MIN_HIST,
                                         ETA_MAX_HIST,kErrorOnMean);
   static PhProfile *dNdetaRow1= hf->CreateProfile("dNdetaRow1",
                                        "dN/d[c]d[f] (event averaged) row 1",
                                         N_ETA_BINS_IN_HIST,ETA_MIN_HIST,
                                         ETA_MAX_HIST,kErrorOnMean);
   static PhProfile *dNdetaRow2= hf->CreateProfile("dNdetaRow2",
                                        "dN/d[c]d[f] (event averaged) row 2",
                                         N_ETA_BINS_IN_HIST,ETA_MIN_HIST,
                                         ETA_MAX_HIST,kErrorOnMean);
   static PhProfile *dNdetaRow3= hf->CreateProfile("dNdetaRow3",
                                        "dN/d[c]d[f] (event averaged) row 3",
                                         N_ETA_BINS_IN_HIST,ETA_MIN_HIST,
                                         ETA_MAX_HIST,kErrorOnMean);
   static PhProfile *dNdetaRow4= hf->CreateProfile("dNdetaRow4",
                                        "dN/d[c]d[f] (event averaged) row 4",
                                         N_ETA_BINS_IN_HIST,ETA_MIN_HIST,
                                         ETA_MAX_HIST,kErrorOnMean);
   static PhProfile *dNdetaRow5= hf->CreateProfile("dNdetaRow5",
                                        "dN/d[c]d[f] (event averaged) row 5",
                                         N_ETA_BINS_IN_HIST,ETA_MIN_HIST,
                                         ETA_MAX_HIST,kErrorOnMean);
   static PhProfile *dNdetaPadCount= hf->CreateProfile("dNdetaPadCount",
                                        "dN/d[c] from pads,event avg,pad count",
                                         N_ETA_BINS_IN_HIST,ETA_MIN_HIST,
                                         ETA_MAX_HIST,kErrorOnMean);
   static PhProfile *dNdetaPadDeconv= hf->CreateProfile("dNdetaPadDeconv",
                                        "dN/d[c] from pads,event avg,deconv",
                                         N_ETA_BINS_IN_HIST,ETA_MIN_HIST,
                                         ETA_MAX_HIST,kErrorOnMean);
   static PhProfile *dNdetaEv =       hf->CreateProfile("dNdetaEv", "dN/d[c] (single event)",
                                         N_ETA_BINS_IN_HIST,ETA_MIN_HIST,
                                         ETA_MAX_HIST,kErrorOnMean);

//   static PhNtuple *MVDNtuple= hf->CreateNtuple("MVDNtuple",
//                                    "results",
//                                    "Zvertex:MeV_SiR1:NPTLS:",
//                         "strhit_1:strhit_2:bfm:ITRKER:Ztrkdif:occR1:occR2:",
//                         "occpad1:occpad2:dsumR1:dsumR2:dsumpad1:dsumpad2:",
//                         "NtrueR1:NtrueR2:Ntruepd1:Ntruepd2:paddcon1:",
//                         "paddecon2:Xtrue:Ytrue:Xcalc:Ycalc");
       
   dNdetaEv->Reset();   // this histograms is for a single event 

   MvdMultiOut  mult_tmp=mvd.BMulti();
   MvdVertex vertex_tmp=mvd.CalcVertex();
   MvdOccupancy occup_tmp=mvd.Occups();
   vector<MvddNdEtaOut>  dndetas_tmp=mvd.DNdEta();
   vector<MvdClump> clumps_tmp=mvd.Clumps();
   vector<MvdNClump> nclmps_tmp=mvd.Nclmps();
   vector<MvdStrip> strips_tmp=mvd.Strips();
   vector<MvdPad> pads_tmp=mvd.Pads();

   // using header table 
   //bimevt=header->get_b(0);             // impact parameter
   //float true_zvertex=header->get_vertex(2,0);      // z of vertex 
   //xnptls=header->get_multiplicity(0);  // true total multiplicity in Geant

//   TrueMult->Fill(xnptls,1.);
//   TrueZvertex->Fill(vertex,1.);
//   bfm->Fill(bimevt,1.);
   
   xvertex_fit=vertex_tmp.get_x();
   yvertex_fit=vertex_tmp.get_y();
   zvertex_fit=vertex_tmp.get_z();
     errorcode=vertex_tmp.get_errorcode();
    softvertex=vertex_tmp.get_softvertex();

   herrorcode->Fill(errorcode);

   //if ( errorcode==1){
   //   if(softvertex>=1.0 && softvertex<2.0){
   //      cout<<"filling dZvertex_pseudo"<<endl;
   //      dZvertex_pseudo->Fill(zvertex_fit-true_zvertex);
   //  }
   //   else if(softvertex>=2.0){
   //      dZvertex_correlation->Fill(zvertex_fit-true_zvertex);
   //   }
   //   else{
   //      cerr<<"unknown vertexing algorithm :"<<softvertex<<endl;
   //   }
   //}

   // Histograms of the panel occupancy vs z position along vertex barrel.
   clsum=0;
   ncl=0;
   vector<const MvdStrip *>::iterator sp;
   vector<const MvdStrip *>::iterator sp_next;

 
   for( i=0;i<clumps_tmp.size();++i){ 
      sp=clumps_tmp[i].Begin();
      (*sp)->Id(shell,row,panel,strip);
      adc=clumps_tmp[i].Adc(); 
//      if(shell==0){
//         Dy=Geo->vislr+Geo->dvisl[0];
//         wid=Geo->dvisl[1];
//      }
//      if (shell==1){
//         Dy=Geo->voslr+Geo->dvosl[0];
//          wid=Geo->dvosl[1];
//      }


      ClumpAdc->Fill((float)adc, 1.);

      ClumpAdc->print();
      zpanelpos=bgeometry.Z(panel);
      offset  =  zpanelpos - zvertex_fit;
      if (shell == 0 )
         ClumpVsZ1->Fill(offset, (float)clumps_tmp[i].Size());
      else
         ClumpVsZ2->Fill(offset, (float)clumps_tmp[i].Size());

      // the end of the list is a special case, it is always
      //   the last clump of the panel it is in
      EndOfPanel = 0;
      if ((i+1)==clumps_tmp.size()) EndOfPanel = 1;
      else{ 
         sp_next=clumps_tmp[i+1].Begin(); 
         (*sp_next)->Id(shell_next,row_next,panel_next,strip_next);

         if (shell!=shell_next || row!=row_next || panel!=panel_next) 
            EndOfPanel=1;
      }
      clsum+=(float)clumps_tmp[i].Size(); // number of strips "on"
      ncl++;                        // number of clumps      

      if ( EndOfPanel == 1 ) {
         paneloccup=clsum/nstrip;

         if ( shell == 0) {
            OccVsZR1->Fill(zpanelpos,paneloccup);
            NclumpVsZ1->Fill(zpanelpos,(float)ncl);
         }
         if ( shell ==1) {
            OccVsZR2->Fill(zpanelpos,paneloccup);
            NclumpVsZ2->Fill(zpanelpos,(float)ncl);
         }
         clsum=0;
         ncl=0;
      }
   }

   // histogram of the total number of clumps found in an event
   // in each of the two barrels

   innerclmpsum=0;
   outerclmpsum=0;
   short number;

   for ( i=0;i< nclmps_tmp.size() ;++i) {
      nclmps_tmp[i].Id(shell,row,panel);
      number=nclmps_tmp[i].Num();
      if ( shell == 0 )
         innerclmpsum+=(float)number;
      if ( shell == 1)
         outerclmpsum+=(float)number;
   }

   Nclumps1->Fill(innerclmpsum,1.);
   Nclumps2->Fill(outerclmpsum,1.);

   // convert to Mev average E loss and histogram
   sumkev=0;
   for ( i=0; i< strips_tmp.size();++i){
      strips_tmp[i].Id(shell,row,panel,strip);
      adc=strips_tmp[i].Adc();
      if ( shell==0 )
         sumkev+=((float)adc)*par.barrel_gain;
   }
   summev = sumkev/1000.;
   MeVvsZvert->Fill(zvertex_fit,summev); 
   MeVinSi->Fill(summev,1.);

   // Histogram occupancy of barrels and endcaps 
   occup_tmp=mvd.Occups();

   OccR1->Fill(occup_tmp.Inner(),1.);
   OccR2->Fill(occup_tmp.Outer(),1.);
   OccZp->Fill(occup_tmp.South(),1.);
   OccZm->Fill(occup_tmp.North(),1.);

   // get the sum of all digitized channels,
   // this is what the trigger will be based on
   float factor;
   float inneradcsum,outeradcsum;
   float northpadsum,southpadsum;

   inneradcsum=outeradcsum=northpadsum=southpadsum=0.;
   factor=par.smax_mip/par.full_scale;

   for (i=0; i< strips_tmp.size();++i) {
      strips_tmp[i].Id(shell,row,panel,strip);
      adc=strips_tmp[i].Adc();
      if (shell == 0)
          inneradcsum+=adc;
      if (shell == 1)
          outeradcsum+=adc;
   }
   mip[0]=factor*inneradcsum;
   mip[1]=factor*outeradcsum;

   AdcSumR1->Fill(mip[0],1.);
   AdcSumR2->Fill(mip[1],1.);

   for (int i=0; i < pads_tmp.size(); ++i) {
      pads_tmp[i].Id(end,wedge,column,row);
      adc=pads_tmp[i].Adc();
      if (end ==0)
         southpadsum+=adc;
      if (end == 1)
         northpadsum+=adc;
   }

   mip[2]=factor*southpadsum;
   mip[3]=factor*northpadsum;

   AdcSumZp->Fill(mip[2],1.);
   AdcSumZm->Fill(mip[3],1.);

   // Adc value histogram 
   // for barrel
   for (i=0;i < strips_tmp.size();++i){
      adc=strips_tmp[i].Adc();
      if ( adc > 0) StripAdc->Fill((float)adc,1.);
   }

   // for pad 
   for (i=0; i <pads_tmp.size(); ++i) {
      pads_tmp[i].Id(end,wedge,column,row);
      adc=pads_tmp[i].Adc();
      diff=zvertex_fit-cgeometry.Z(end);  // using calculated vertex
      diff=diff >0 ? diff : -diff;
      if ( diff>par.tooclose ) {
         tantheta=cgeometry.Radius(row)/diff;
         cosring=1./(float)sqrt((double)(1.0+tantheta*tantheta));
         if (adc >0 ) {
//             xadc=adc*cosring;
             xadc=adc;
             PadAdc1->Fill(xadc,1.);


             // histogram Adc distributions for different occupancies in pads
             if (end==0) occupancy=occup_tmp.Inner();
             else occupancy=occup_tmp.Outer(); 
             if ( occupancy <= 0.01 )
                PadAdcocclt1->Fill(xadc,1.);
             else if (occupancy <= .02 )
                PadAdcocc1to2->Fill(xadc,1.);
             else if (occupancy <= .05 )
                PadAdcocc2to5->Fill(xadc,1.);
             else if ( occupancy<= .10 )
                PadAdcocc5to10->Fill(xadc,1.);
             else
                PadAdcoccgt10->Fill(xadc,1.);
         }
      }
   }


   // make some various dN/deta histograms 

   // now fill them: 
   dndetas_tmp=mvd.DNdEta();

   //for (k=0;k<dndetas_tmp.size();++k){
   //   dndetas_tmp[k].Show();
   //   cout<<endl;
   //}
 
   if ( dndetas_tmp.size() > 0 ) {
      for (k=0;k<=dndetas_tmp.size();++k) {
         dndetas_tmp[k].Put(dndeta,dndetaerr,eta,deta,phi,dphi,softdndeta);

         if ( softdndeta>1.0 && softdndeta<2.0 ) {
           // this section is for dN/deta from the barrel with all rows
           // combined into a single eta bin                           
           dNdetaAvg->Fill(eta,dndeta);
//           if ( Head->event<=MAX_SINGLE_EV_OUT )
//              dNdetaEv->Fill( eta, dndeta);
         }
         else if (softdndeta>2.0 && softdndeta<3.0 ) {
           // this section is for dN/deta/dphi from individual rows 
           //in the barrel

           if ( phi< 0.5236 && phi>-0.52366 )
              dNdetaRow0->Fill(eta,dndeta);
           else if ( phi < 1.5708 )
              dNdetaRow1->Fill(eta,dndeta);
           else if ( phi < 2.6180 )
              dNdetaRow2->Fill(eta,dndeta);
           else if ( phi < 3.6652 )
              dNdetaRow3->Fill(eta,dndeta);
           else if ( phi < 4.7124 )
              dNdetaRow4->Fill(eta,dndeta);
           else if ( phi < 5.7596 )
              dNdetaRow5->Fill(eta,dndeta);
         }
         else if (softdndeta>=3.0 && softdndeta<4.0 ) {
           // this is the section for dN/deta from the average dE/dx
           //   calculation in the pad detectors
           dNdetaPadCount->Fill(eta,dndeta);
         }
         else if (softdndeta>=4.0 && softdndeta<5.0){
            // this is the section for dN/deta from the deconvolution method
            // in the pad detectors
            dNdetaPadDeconv->Fill(eta,dndeta);
         }
      }
   }
}