/[PAMELA software]/PamCut/CaloAxis2.h
ViewVC logotype

Diff of /PamCut/CaloAxis2.h

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1 by pam-fi, Wed May 27 13:30:01 2009 UTC revision 1.2 by pam-fi, Wed Aug 5 14:07:04 2009 UTC
# Line 3  Line 3 
3    
4  #if !defined(__CINT__) || defined(__MAKECINT__)  #if !defined(__CINT__) || defined(__MAKECINT__)
5    
   
6  #include <TString.h>  #include <TString.h>
7  #include <TH1F.h>  #include <TH1F.h>
8  #include <TH2F.h>  #include <TH2F.h>
# Line 27  using namespace std; Line 26  using namespace std;
26  //  //
27  //  //
28  //===============================================================================  //===============================================================================
29  class CaloEvent{  class CaloEvent {
30    
31  public:  public:
32  //  calo event    //  calo event
33  //  ------------------------------    //  ------------------------------
34  //  hit strips    //  hit strips
35  //  ------------------------------    //  ------------------------------
36      vector<int> ids;  //strip 0-95    vector<int> ids; //strip 0-95
37      vector<int> idp;  //plane 0-21    vector<int> idp; //plane 0-21
38      vector<float> q; // signal    vector<float> q; // signal
39      vector<float> zcoord;    vector<float> zcoord;
40      vector<float> xycoord;    vector<float> xycoord;
41    
42      //  ------------------------------
43  //  ------------------------------    //  methods
44  //  methods    //  ------------------------------
45  //  ------------------------------    void Add(int iis, int iip, float fq, float fz, float fxy) {
46      void Add(int iis, int iip, float fq, float fz, float fxy){      ids.push_back(iis);
47          ids.push_back(iis);      idp.push_back(iip);
48          idp.push_back(iip);      q.push_back(fq);
49          q.push_back(fq);      zcoord.push_back(fz);
50          zcoord.push_back(fz);      xycoord.push_back(fxy);
51          xycoord.push_back(fxy);    }
52      };    ;
53    
54      int GetN() {
55        return (int) (q.size());
56      int  GetN(){return (int)(q.size());};    }
57      ;
58      void Reset(){  
59          ids.clear();    void Reset() {
60          idp.clear();      ids.clear();
61          q.clear();      idp.clear();
62          zcoord.clear();      q.clear();
63          xycoord.clear();      zcoord.clear();
64      };      xycoord.clear();
65      }
66      void Delete(){ Reset(); };    ;
67    
68      ~CaloEvent(){ Delete(); };    void Delete() {
69        Reset();
70      CaloEvent(){ };    }
71      ;
72      CaloEvent( CaloLevel1* calo, int view, Bool_t usemechanicalalignement){ Set(calo,view,usemechanicalalignement); };  
73      ~CaloEvent() {
74      CaloEvent( CaloLevel1* calo, int view){ Set(calo,view); };      Delete();
75      }
76      void Set(CaloLevel1* calo, int view, Bool_t usemechanicalalignement){    ;
77    
78          Reset();    CaloEvent() {
79  //      cout << " CaloEvent::Set()"<<endl;    }
80          CaloStrip cstrip;    ;
81          cstrip = CaloStrip(calo,usemechanicalalignement);  
82      CaloEvent(CaloLevel1* calo, int view, Bool_t usemechanicalalignement) {
83          for(int ih=0;ih<calo->istrip;ih++){      Set(calo, view, usemechanicalalignement);
84              int iv=-1;    }
85              int ip=-1;    ;
86              int is=-1;  
87              calo->DecodeEstrip(ih,iv,ip,is);    CaloEvent(CaloLevel1* calo, int view) {
88              if(iv==view){      Set(calo, view);
89                  cstrip.Set(view,ip,is);    }
90                  float fxy = 0.;    ;
91                  if( !view ) fxy = cstrip.GetX();  
92                  else        fxy = cstrip.GetY();    void Set(CaloLevel1* calo, int view, Bool_t usemechanicalalignement) {
93                  float fz = cstrip.GetZ();  
94                  float fq = cstrip.GetE();      Reset();
95                  if(fq>0)Add(is,ip,fq,fz,fxy);      //  cout << " CaloEvent::Set()"<<endl;
96              }      CaloStrip cstrip;
97          }      cstrip = CaloStrip(calo, usemechanicalalignement);
98    
99      };      for (int ih = 0; ih < calo->istrip; ih++) {
100      void Set(CaloLevel1* calo, int view){ Set(calo,view,0); };        int iv = -1;
101          int ip = -1;
102          int is = -1;
103          calo->DecodeEstrip(ih, iv, ip, is);
104          if (iv == view) {
105            cstrip.Set(view, ip, is);
106            float fxy = 0.;
107            if (!view)
108              fxy = cstrip.GetX();
109            else
110              fxy = cstrip.GetY();
111            float fz = cstrip.GetZ();
112            float fq = cstrip.GetE();
113            if (fq > 0)
114              Add(is, ip, fq, fz, fxy);
115          }
116        }
117    
118      }
119      ;
120      void Set(CaloLevel1* calo, int view) {
121        Set(calo, view, 0);
122      }
123      ;
124    
125  };  };
126  //===============================================================================  //===============================================================================
# Line 112  class CaloAxis { Line 133  class CaloAxis {
133    
134  public:  public:
135    
136      CaloEvent *cevent;    CaloEvent *cevent;
   
 //  ------------------------------  
 //  selected points along the axis  
 //  ------------------------------  
 //  fitted points  
     vector<float> x;// independent coordinates (z)  
     vector<float> y;// dependente coordinate (x/y)  
     vector<float> w;// fit weight  
     vector<float> q;// signal  
   
 //  ------------------------------  
 //  axis parameters  
 //  ------------------------------  
     Float_t par[2];  
     Float_t covpar[2][2];  
   
 //  ------------------------------  
 //  general variables  
 //  ------------------------------  
     float cog[22];// baricenter  
     float qpl[22]; // charge in each plane  
   
 //  ------------------------------  
 //  objects to draw the axis  
 //  ------------------------------  
     TPolyMarker *sel;  
     TLine       *axis;  
   
   
 //  ------------------------------  
 //  methods  
 //  ------------------------------  
     void Init(){  
         for(Int_t i=0; i<2;i++){  
             par[i]=0;  
             for(Int_t j=0; j<2;j++)covpar[i][j]=0;  
         };  
         sel  = 0;  
         axis = 0;  
         for(Int_t i=0; i<22; i++){  
             cog[i]=0;  
             qpl[i]=0;  
         }  
 //      if(cevent)cevent->Reset();  
     }  
   
     void Set(CaloLevel1* calo, Int_t view){ cevent->Set(calo,view,0); };  
     void Set(CaloLevel1* calo, Int_t view, Bool_t usemechanicalalignement){ cevent->Set(calo,view,usemechanicalalignement); };  
   
     CaloAxis(){  
         cevent = new CaloEvent();  
         Init();  
     };  
   
     CaloAxis(CaloLevel1* calo, Int_t view){  
         cevent = new CaloEvent();  
         Init();  
         Set(calo,view);  
     };  
     CaloAxis(CaloLevel1* calo, Int_t view, Bool_t usemechanicalalignement){  
         cevent = new CaloEvent();  
         Init();  
         Set(calo,view,usemechanicalalignement);  
     };  
   
   
     void Reset(){  
 //      cout << " CaloAxis::Reset()"<<endl;  
         x.clear();  
         y.clear();  
         w.clear();  
         q.clear();  
 //      cout << sel << endl;  
         if(sel)sel->Delete();  
 //      cout << axis << endl;  
         if(axis)axis->Delete();  
 //      if(cevent)cevent->Reset();  
         Init();  
     }  
   
     void Delete(){ Reset(); delete cevent; };  
   
     ~CaloAxis(){ Delete(); };  
   
   
     void Add(float xin, float yin){  
         x.push_back(xin);  
         y.push_back(yin);  
         w.push_back(1.);  
         q.push_back(0.);  
 /*      Int_t nin = x.size(); */  
 /*      if(nin>2)Fit(); */  
     };  
     void Add(float xin, float yin, float qin){  
         x.push_back(xin);  
         y.push_back(yin);  
         q.push_back(qin);  
         w.push_back(1.);  
 /*      Int_t nin = x.size(); */  
 /*      if(nin>2)Fit(); */  
     };  
     void Add(float xin, float yin, float qin, float win){  
         x.push_back(xin);  
         y.push_back(yin);  
         w.push_back(win);  
         q.push_back(qin);  
 /*      Int_t nin = x.size(); */  
 /*      if(nin>2)Fit(); */  
     };  
   
     int Fit(){  
   
 //      cout << " CaloAxis::Fit()"<<endl;  
         Float_t SSS = 0;  
         Float_t SSX = 0;  
         Float_t SXX = 0;  
   
         Float_t SSY = 0;  
         Float_t SXY = 0;  
   
         Int_t np=0;  
         for(Int_t i=0; i<(int)(x.size()); i++){  
             SSS += w[i]*w[i];  
             SSX += x[i]*w[i]*w[i];  
             SXX += x[i]*x[i]*w[i]*w[i];  
             SSY += y[i]*w[i]*w[i];  
             SXY += x[i]*y[i]*w[i]*w[i];  
             if(w[i])np++;  
         }  
         Float_t det = SSS*SXX-SSX*SSX;  
         if(det==0)return 0;  
         if(np<3)return 0;  
 //      cout << np << " points fitted -- ok "<<endl;  
         par[0] = (SSY*SXX-SXY*SSX)/det;  
         par[1] = (SXY*SSS-SSY*SSX)/det;  
         return 1;  
     };  
   
     void Print(){  
 // useful for debug     for(Int_t i=0; i<(int)(x.size()); i++)if(w[i])cout << x[i] << " - "<<y[i]<<endl;  
     }  
   
     float GetPar(Int_t i){return par[i];};  
   
     int   GetN(){return (int)(x.size());};  
   
     float* GetX(){  
         float *xout = new float[(int)(x.size())];  
         for(Int_t i=0; i<(int)(x.size()); i++)xout[i]=x[i];  
         return xout;  
     }  
   
     float* GetY(){  
         float *yout = new float[(int)(y.size())];  
         for(Int_t i=0; i<(int)(y.size()); i++) yout[i]=y[i];  
         return yout;  
     }  
   
     float* GetQ(){  
         float *qout = new float[(int)(q.size())];  
         for(int i=0; i<(int)(q.size()); i++) qout[i]=q[i];  
         return qout;  
     }  
   
     float GetChi2(){  
         float chi2=0;  
         int   nchi=0;  
         for(Int_t i=0; i<(int)(y.size()); i++){  
             chi2 += w[i]*w[i]*(y[i]-par[0]-par[1]*x[i])*(y[i]-par[0]-par[1]*x[i]);  
             nchi += (int)w[i];  
         };  
   
         if(nchi<3) return -999; // original  
 //      if(nchi<3 || chi2==0)return -999; // modified by Sergio  
   
         chi2=chi2/(nchi-2);  
         return chi2;  
     }  
   
     float GetQaxis(){  
         float qaxis=0;  
         for(Int_t i=0; i<(int)(q.size()); i++) qaxis +=q[i];  
 //      for(Int_t i=0; i<(int)(q.size()); i++) cout<<q[i]<<endl;  
         return qaxis;  
     }  
   
     float GetQaxis(float toll){  
         float qaxis=0;  
         for(int ih=0;ih<cevent->GetN(); ih++){  
             float x = cevent->xycoord[ih];  
             float z = cevent->zcoord[ih];  
             float q = cevent->q[ih];  
             //int ip = cevent->idp[ih];  
             float d = fabs(x-par[0]-par[1]*z)/sqrt(par[1]*par[1]+1);  
             if( d < toll )qaxis+=q;  
         }  
         return qaxis;  
     }  
   
     float GetCOG(Int_t ip){  
         if(ip>=0 && ip<22)return cog[ip];  
         else return 0;  
     }  
   
     float GetQ(Int_t ip){  
         if(ip>=0 && ip<22)return qpl[ip];  
         else return 0;  
     }  
   
     float GetQstrip(){  
   
         float qaxis=0;  
         for(Int_t i=0; i<(int)(q.size()); i++) qaxis +=q[i];  
         if(q.size()==0) return 0.;  
         return qaxis/q.size();  
     }  
   
     float GetYfit(float x){  
         return par[0]+x*par[1];  
     }  
   
     float GetXmin(){  
         float zmin=9999;  
         for(Int_t i=0; i<(int)(x.size()); i++) if(x[i]<zmin) zmin = x[i];  
         return zmin;  
     }  
     float GetXmax(){  
         float zmax=-9999;  
         for(Int_t i=0; i<(int)(x.size()); i++) if(x[i]>zmax) zmax = x[i];  
         return zmax;  
   
     }  
   
   
     void Draw(){  
   
 //      cout << " CaloAxis::Draw()"<<endl;  
         if(GetN()<3)return;  
   
         sel = new TPolyMarker(GetN(),GetY(),GetX());  
         sel->SetMarkerSize(0.3);  
         sel->SetMarkerColor(5);  
   
         axis = new TLine( par[0]+GetXmin()*par[1], GetXmin(), par[0]+GetXmax()*par[1], GetXmax() );  
         axis->SetLineWidth(1);  
         axis->SetLineColor(5);  
   
 //      cout << sel << endl;  
 //      cout << axis << endl;  
   
         sel->Draw();  
         axis->Draw();  
     };  
   
   
     //-----------------------------------------------------------------  
     // fit the axis (optimized for non-interacting patterns)  
     //-----------------------------------------------------------------  
     int FitAxis( CaloLevel1* calo , Int_t view , Float_t toll , Bool_t usemechanicalalignement ){  
   
 //      cout << "CaloAxis::FitAxis(...)"<<endl;  
   
         Set(calo,view,usemechanicalalignement);  
   
         // ------------------------------  
         // first : all hits included, w=1  
         // ------------------------------  
         for(int ih=0;ih<cevent->GetN(); ih++){  
             float x = cevent->xycoord[ih];  
             float z = cevent->zcoord[ih];  
             float q = cevent->q[ih];  
             Add(z,x,q);  
         }  
   
         if(GetN()<3)return 0;  
         Fit();  
 // useful for debug     cout << "(1)"<<endl;  
 //      cout << " n. hits :"<<GetN()<<endl;  
 //      cout << " P0 P1   :"<<par[0]<< " "<<par[1]<<endl;  
 //      cout << " Chi2    :"<<GetChi2()<<endl;  
 //      cout << " Q axis  :"<<GetQaxis()<<endl;  
 //      cout << " Q strip :"<<GetQstrip()<<endl;  
   
         float P0;  
         float P1;  
         int   dof;  
   
         // ---------------------------  
         // second : iteration  
         // ---------------------------  
         int niter=0;                      // number of iterations  
         int pok[22];  
   
         // prova  
         // pesi lineari  
         float W22=1.;  
         float W1=10.;  
         float b=(W22-W1)/(22.-1.);// (w22 - w1)/(22-1)  
         float a=W1-b;  
   
         for(int i=0; i<22; i++)pok[i]=0;  
         do{  
             niter++;  
             float ttoll = toll;           //tolerance (cm)  
             if(niter==1)ttoll = 10.*toll;  
             dof = GetN();           //previous fit  
             P0  = par[0];  
             P1  = par[1];  
             Reset();                //fit reset  
             for(int ih=0;ih<cevent->GetN(); ih++){//loop over selected hits  
                 float x = cevent->xycoord[ih];  
                 float z = cevent->zcoord[ih];  
                 float q = cevent->q[ih];  
                 int ip = cevent->idp[ih];  
                 float d = fabs(x-P0-P1*z)/sqrt(P1*P1+1);//distance to axis  
 //              cout<<d<<endl;  
                 if( d < ttoll && (niter==1 || (niter>1 && pok[ip]==1)) ){  
 //              if( d < ttoll ){  
                     float W=a+b*(ip+1);  
 //                  cout << ip << " "<<W<<endl;  
                     Add(z,x,q,W);  
                     pok[ip]=1;  
                 }  
             }  
 // break the track if more than 3 planes are missing  
             int   nmissing = 0;  
             for(int ip=0; ip<22; ip++){  
                 if(pok[ip]==0){  
                     nmissing++;  
                 }else{  
                     nmissing = 0;  
                 }  
                 if(nmissing==3){  
                     for(int ipp=ip+1; ipp<22; ipp++)pok[ipp]=0;  
                     break;  
                 }  
             }  
   
   
             if(niter==100)break;  
             if(GetN()<3)break;  
             Fit();  
         }while(niter==1 || GetN() != dof);  
 //      cout << " >>> "<<GetN()<<endl;  
         if(GetN()<3)return 0;  
         Fit();  
 //      cout << "(2)"<<endl;  
 //      cout << " n. hits :"<<GetN()<<endl;  
 //      cout << " P0 P1   :"<<par[0]<< " "<<par[1]<<endl;  
 //      cout << " Chi2    :"<<GetChi2()<<endl;  
 //      cout << " Q axis  :"<<GetQaxis()<<endl;  
 //      cout << " Q strip :"<<GetQstrip()<<endl;  
   
         // ---------------------------------------------  
         // third : selecting only closest strip-clusters  
         // and fit baricenters  
         // ---------------------------------------------  
         P0  = par[0];  
         P1  = par[1];  
         Reset();  
   
         float dmin[22];  
         int closest[22];  
         for(int ip=0; ip<22; ip++){  
             dmin[ip]=999;  
             closest[ip]=-1;  
         }  
         for(int ih=0;ih<cevent->GetN(); ih++){  
             float x = cevent->xycoord[ih];  
             float z = cevent->zcoord[ih];  
             //float q = cevent->q[ih];  
             int ip = cevent->idp[ih];  
             //int is = cevent->ids[ih];  
             float d = fabs(x-P0-P1*z)/sqrt(P1*P1+1);  
             if( d < toll && d<dmin[ip] && pok[ip]==1 ){  
 //              Add(z,x,q,q);  
                 closest[ip]=ih;  
 //              cog[ip] += x*q;  
 //              qpl[ip] += q;  
             }  
         }  
         for(Int_t ip=0; ip<22; ip++){  
             if(closest[ip]!=-1){  
                 float x = cevent->xycoord[closest[ip]];  
                 float z = cevent->zcoord[closest[ip]];  
                 float q = cevent->q[closest[ip]];  
                 //int is = cevent->ids[closest[ip]];  
                 Add(z,x,q,q);  
                 cog[ip] += x*q;  
                 qpl[ip] += q;  
 //              cout << ip << " -o- "<<is<<endl;  
             }  
         }  
         // -----------------  
         // add +/- one strip  
         // -----------------  
         for(int ih=0;ih<cevent->GetN(); ih++){  
             float x = cevent->xycoord[ih];  
             float z = cevent->zcoord[ih];  
             float q = cevent->q[ih];  
             int ip = cevent->idp[ih];  
             int is = cevent->ids[ih];  
             if( closest[ip]!=-1 ){  
                 int isc = cevent->ids[closest[ip]];  
                 if( is == isc+1 || is == isc-1 ){  
                     Add(z,x,q,q);  
                     cog[ip] += x*q;  
                     qpl[ip] += q;  
 //                  cout << ip << " -+- "<<is<<endl;  
                 }  
             }  
         }  
         Fit();  
   
         if(GetN()<3)return 0;  
         for(int ip=0; ip<22; ip++){  
             if(qpl[ip]!=0){  
                 cog[ip]=cog[ip]/qpl[ip];  
 //          Add(z,cog[ip],cog[ip],0.7);  
             }  
         }  
   
 //      cout << "(3)"<<endl;  
 //      cout << " n. hits :"<<GetN()<<endl;  
 //      cout << " P0 P1   :"<<par[0]<< " "<<par[1]<<endl;  
 //      cout << " Chi2    :"<<GetChi2()<<endl;  
 //      cout << " Q axis  :"<<GetQaxis()<<endl;  
 //      cout << " Q strip :"<<GetQstrip()<<endl;  
 //      cout << " ---------------------------"<<endl;  
         return 1;  
     };  
   
   
     int FitAxis( CaloLevel1* calo , Int_t view , Float_t toll ){ return FitAxis(calo,view,toll,0); };  
   
   
   
     //-----------------------------------------------------------------  
     // fit the axis (optimized for interacting patterns?? ...forse)  
     //-----------------------------------------------------------------  
     int FitShower( CaloLevel1* calo , Int_t view , Float_t toll , Bool_t usemechanicalalignement ){  
   
 //      cout << "CaloAxis::FitAxis(...)"<<endl;  
   
         Set(calo,view,usemechanicalalignement);  
   
         // ------------------------------  
         // first :  
         // ------------------------------  
         for(int ih=0;ih<cevent->GetN(); ih++){  
             float x = cevent->xycoord[ih];  
             //float z = cevent->zcoord[ih];  
             float q = cevent->q[ih];  
             int ip  = cevent->idp[ih];  
             cog[ip] += x*q;  
             qpl[ip] += q;  
 //          Add(z,x,q);  
         }  
         for(int ip=0; ip<22; ip++){  
             if(qpl[ip]!=0)cog[ip]=cog[ip]/qpl[ip];  
 //          cout << ip << " > "<<qpl[ip]<<" "<<cog[ip]<<endl;  
         }  
         // ----------------  
         // plane of maximum  
         // ----------------  
         float qplmax=1;  
         int ipmax = -1;  
         float dmin[22];  
         int closest[22];  
         for(int ip=0; ip<22; ip++)if(qpl[ip]>qplmax){  
             ipmax=ip;  
             qplmax=qpl[ip];  
             dmin[ip]=1000.;//init  
             closest[ip]=-1;//init  
         }  
         for(int ih=0;ih<cevent->GetN(); ih++){  
             float x = cevent->xycoord[ih];  
             float z = cevent->zcoord[ih];  
             float q = cevent->q[ih];  
             int ip  = cevent->idp[ih];  
 //          if( (ipmax>3 && ip<=ipmax) || ipmax<=3 )Add(z,x,q);  
             if( ip<=ipmax || ip<=3 ){  
                 Add(z,x,q);  
             }  
         }  
         if(GetN()<3)return 0;  
         Fit();  
 //      cout << " n. hits :"<<GetN()<<endl;  
 //      cout << " P0 P1   :"<<par[0]<< " "<<par[1]<<endl;  
 //      cout << " Chi2    :"<<GetChi2()<<endl;  
 //      cout << " Q axis  :"<<GetQaxis()<<endl;  
 //      cout << " Q strip :"<<GetQstrip()<<endl;  
   
         float P0;  
         float P1;  
         int   dof;  
   
         // ---------------------------  
         // second : iteration  
         // ---------------------------  
         int niter=0;                      // number of iterations  
         int pok[22];  
         // prova  
         // pesi lineari  
         float W22=1.;  
         float W1=10.;  
         float b=(W22-W1)/(22.-1.);// (w22 - w1)/(22-1)  
         float a=W1-b;  
   
         for(int i=0; i<22; i++){  
             pok[i]=0;  
         }  
         do{  
             niter++;  
             float ttoll = toll;           //tolerance (cm)  
             if(niter==1)ttoll = 10.*toll;  
             dof = GetN();           //previous fit  
             P0  = par[0];  
             P1  = par[1];  
             Reset();                //fit reset  
             for(int i=0; i<22; i++){  
                 cog[i]=0;  
                 qpl[i]=0;  
             }  
             for(int ih=0;ih<cevent->GetN(); ih++){  
                 float x = cevent->xycoord[ih];  
                 float z = cevent->zcoord[ih];  
                 float q = cevent->q[ih];  
                 int ip = cevent->idp[ih];  
                 float d = fabs(x-P0-P1*z)/sqrt(P1*P1+1);  
 //              cout<<d<<endl;  
                 if(  
                     (niter==1 || (niter>1 && pok[ip]==1)) &&  
                     (ip<=ipmax || ip<=3) &&  
 //                  ((ipmax>3 && ip<=ipmax) || ipmax<=3) &&  
                     true){  
   
                     if(niter==1 && d<dmin[ip]){  
                         dmin[ip]=d;  
                         closest[ip]=ih;  
                     }  
   
                     float W=a+b*(ip+1);  
 //                  cout << ip << " "<<W<<endl;  
                     if( d < ttoll || ( niter==2 && ih==closest[ip] )){  
 /*                      if(ip==ipmax && ipmax>1)Add(z,x,q,W*q);  */  
 /*                      else Add(z,x,q,W);  */  
                         cog[ip] += x*q;  
                         qpl[ip] += q;  
                         Add(z,x,q,W);  
                         pok[ip]=1;  
                         cout << ip << " "<<qpl[ip]<<endl;  
                     }  
                 }  
             }  
 // break the track if more than 3 planes are missing  
             int   nmissing = 0;  
             for(int ip=0; ip<22; ip++){  
                 if(pok[ip]==0){  
                     nmissing++;  
                 }else{  
                     nmissing = 0;  
                 }  
                 if(nmissing==6){  
                     for(int ipp=ip+1; ipp<22; ipp++)pok[ipp]=0;  
                     break;  
                 }  
             }  
   
             for(int ip=0; ip<22; ip++)cout << qpl[ip] << " ";  
             cout << endl;  
   
   
             for(int ip=0; ip<22; ip++)if(qpl[ip]!=0)cog[ip]=cog[ip]/qpl[ip];  
   
   
             if(niter==100)break;  
             if(GetN()<3)break;  
   
             Fit();  
   
         }while(niter==1 || GetN() != dof);  
 //      cout << " >>> "<<GetN()<<endl;  
         if(GetN()<3)return 0;  
         Fit();  
 /* //   cout << " n. hits :"<<GetN()<<endl; */  
 /* //   cout << " P0 P1   :"<<par[0]<< " "<<par[1]<<endl; */  
 /* //   cout << " Chi2    :"<<GetChi2()<<endl; */  
 /* //   cout << " Q axis  :"<<GetQaxis()<<endl; */  
 /* //   cout << " Q strip :"<<GetQstrip()<<endl; */  
   
 /*      // --------------------------------------------- */  
 /*      // third : selecting only closest strip-clusters */  
 /*      // and fit baricenters */  
 /*      // --------------------------------------------- */  
 /*      P0  = par[0]; */  
 /*      P1  = par[1]; */  
 /*      Reset(); */  
   
 /*      float dmin[22]; */  
 /*      int closest[22]; */  
 /*      for(int ip=0; ip<22; ip++){ */  
 /*          dmin[ip]=999; */  
 /*          closest[ip]=-1; */  
 /*      } */  
 /*      for(int ih=0;ih<cevent->GetN(); ih++){ */  
 /*          float x = cevent->xycoord[ih]; */  
 /*          float z = cevent->zcoord[ih]; */  
 /*          float q = cevent->q[ih]; */  
 /*          int ip = cevent->idp[ih]; */  
 /*          int is = cevent->ids[ih]; */  
 /*          float d = fabs(x-P0-P1*z)/sqrt(P1*P1+1); */  
 /*          if( d < toll && d<dmin[ip] && pok[ip]==1 ){ */  
 /* //           Add(z,x,q,q); */  
 /*              closest[ip]=ih; */  
 /* //           cog[ip] += x*q; */  
 /* //           qpl[ip] += q; */  
 /*          } */  
 /*      } */  
 /*      for(Int_t ip=0; ip<22; ip++){ */  
 /*          if(closest[ip]!=-1){ */  
 /*              float x = cevent->xycoord[closest[ip]]; */  
 /*              float z = cevent->zcoord[closest[ip]]; */  
 /*              float q = cevent->q[closest[ip]]; */  
 /*              int is = cevent->ids[closest[ip]]; */  
 /*              Add(z,x,q,q); */  
 /*              cog[ip] += x*q; */  
 /*              qpl[ip] += q; */  
 /* //           cout << ip << " -o- "<<is<<endl; */  
 /*          } */  
 /*      } */  
 /*      // add +/- one strip */  
 /*      for(int ih=0;ih<cevent->GetN(); ih++){ */  
 /*          float x = cevent->xycoord[ih]; */  
 /*          float z = cevent->zcoord[ih]; */  
 /*          float q = cevent->q[ih]; */  
 /*          int ip = cevent->idp[ih]; */  
 /*          int is = cevent->ids[ih]; */  
 /*          if( closest[ip]!=-1 ){ */  
 /*              int isc = cevent->ids[closest[ip]]; */  
 /*              if( is == isc+1 || is == isc-1 ){ */  
 /*                  Add(z,x,q,q); */  
 /*                  cog[ip] += x*q; */  
 /*                  qpl[ip] += q; */  
 /* //               cout << ip << " -+- "<<is<<endl; */  
 /*              } */  
 /*          } */  
 /*      } */  
 /*      Fit(); */  
   
 /*      if(GetN()<3)return 0; */  
 /*      for(int ip=0; ip<22; ip++){ */  
 /*          if(qpl[ip]!=0){ */  
 /*              cog[ip]=cog[ip]/qpl[ip]; */  
 /* //       Add(z,cog[ip],cog[ip],0.7); */  
 /*          } */  
 /*      } */  
   
 //      cout << " n. hits :"<<GetN()<<endl;  
 //      cout << " P0 P1   :"<<par[0]<< " "<<par[1]<<endl;  
 //      cout << " Chi2    :"<<GetChi2()<<endl;  
 //      cout << " Q axis  :"<<GetQaxis()<<endl;  
 //      cout << " Q strip :"<<GetQstrip()<<endl;  
 //      cout << " ---------------------------"<<endl;  
         return 1;  
     };  
   
   
     int FitShower( CaloLevel1* calo , Int_t view , Float_t toll ){ return FitAxis(calo,view,toll,0); };  
137    
138      //  ------------------------------
139      //  selected points along the axis
140      //  ------------------------------
141      //  fitted points
142      vector<float> x;// independent coordinates (z)
143      vector<float> y;// dependente coordinate (x/y)
144      vector<float> w;// fit weight
145      vector<float> q;// signal
146    
147      //  ------------------------------
148      //  axis parameters
149      //  ------------------------------
150      Float_t par[2];
151      Float_t covpar[2][2];
152    
153      //  ------------------------------
154      //  general variables
155      //  ------------------------------
156      float cog[22];// baricenter
157      float qpl[22]; // charge in each plane
158    
159      //  ------------------------------
160      //  objects to draw the axis
161      //  ------------------------------
162      TPolyMarker *sel;
163      TLine *axis;
164    
165      //  ------------------------------
166      //  methods
167      //  ------------------------------
168      void Init() {
169        for (Int_t i = 0; i < 2; i++) {
170          par[i] = 0;
171          for (Int_t j = 0; j < 2; j++)
172            covpar[i][j] = 0;
173        };
174        sel = 0;
175        axis = 0;
176        for (Int_t i = 0; i < 22; i++) {
177          cog[i] = 0;
178          qpl[i] = 0;
179        }
180        //  if(cevent)cevent->Reset();
181      }
182    
183      void Set(CaloLevel1* calo, Int_t view) {
184        cevent->Set(calo, view, 0);
185      }
186      ;
187      void Set(CaloLevel1* calo, Int_t view, Bool_t usemechanicalalignement) {
188        cevent->Set(calo, view, usemechanicalalignement);
189      }
190      ;
191    
192      CaloAxis() {
193        cevent = new CaloEvent();
194        Init();
195      }
196      ;
197    
198      CaloAxis(CaloLevel1* calo, Int_t view) {
199        cevent = new CaloEvent();
200        Init();
201        Set(calo, view);
202      }
203      ;
204      CaloAxis(CaloLevel1* calo, Int_t view, Bool_t usemechanicalalignement) {
205        cevent = new CaloEvent();
206        Init();
207        Set(calo, view, usemechanicalalignement);
208      }
209      ;
210    
211      void Reset() {
212        //  cout << " CaloAxis::Reset()"<<endl;
213        x.clear();
214        y.clear();
215        w.clear();
216        q.clear();
217        //  cout << sel << endl;
218        if (sel)
219          sel->Delete();
220        //  cout << axis << endl;
221        if (axis)
222          axis->Delete();
223        //  if(cevent)cevent->Reset();
224        Init();
225      }
226    
227      void Delete() {
228        Reset();
229        delete cevent;
230      }
231      ;
232    
233      ~CaloAxis() {
234        Delete();
235      }
236      ;
237    
238      void Add(float xin, float yin) {
239        x.push_back(xin);
240        y.push_back(yin);
241        w.push_back(1.);
242        q.push_back(0.);
243        /*  Int_t nin = x.size(); */
244        /*  if(nin>2)Fit(); */
245      }
246      ;
247      void Add(float xin, float yin, float qin) {
248        x.push_back(xin);
249        y.push_back(yin);
250        q.push_back(qin);
251        w.push_back(1.);
252        /*  Int_t nin = x.size(); */
253        /*  if(nin>2)Fit(); */
254      }
255      ;
256      void Add(float xin, float yin, float qin, float win) {
257        x.push_back(xin);
258        y.push_back(yin);
259        w.push_back(win);
260        q.push_back(qin);
261        /*  Int_t nin = x.size(); */
262        /*  if(nin>2)Fit(); */
263      }
264      ;
265    
266      int Fit() {
267    
268        //  cout << " CaloAxis::Fit()"<<endl;
269        Float_t SSS = 0;
270        Float_t SSX = 0;
271        Float_t SXX = 0;
272    
273        Float_t SSY = 0;
274        Float_t SXY = 0;
275    
276        Int_t np = 0;
277        for (Int_t i = 0; i < (int) (x.size()); i++) {
278          SSS += w[i] * w[i];
279          SSX += x[i] * w[i] * w[i];
280          SXX += x[i] * x[i] * w[i] * w[i];
281          SSY += y[i] * w[i] * w[i];
282          SXY += x[i] * y[i] * w[i] * w[i];
283          if (w[i])
284            np++;
285        }
286        Float_t det = SSS * SXX - SSX * SSX;
287        if (det == 0)
288          return 0;
289        if (np < 3)
290          return 0;
291        //  cout << np << " points fitted -- ok "<<endl;
292        par[0] = (SSY * SXX - SXY * SSX) / det;
293        par[1] = (SXY * SSS - SSY * SSX) / det;
294        return 1;
295      }
296      ;
297    
298      void Print() {
299        // useful for debug for(Int_t i=0; i<(int)(x.size()); i++)if(w[i])cout << x[i] << " - "<<y[i]<<endl;
300      }
301    
302      float GetPar(Int_t i) {
303        return par[i];
304      }
305      ;
306    
307      int GetN() {
308        return (int) (x.size());
309      }
310      ;
311    
312      float* GetX() {
313        float *xout = new float[(int) (x.size())];
314        for (Int_t i = 0; i < (int) (x.size()); i++)
315          xout[i] = x[i];
316        return xout;
317      }
318    
319      float* GetY() {
320        float *yout = new float[(int) (y.size())];
321        for (Int_t i = 0; i < (int) (y.size()); i++)
322          yout[i] = y[i];
323        return yout;
324      }
325    
326      float* GetQ() {
327        float *qout = new float[(int) (q.size())];
328        for (int i = 0; i < (int) (q.size()); i++)
329          qout[i] = q[i];
330        return qout;
331      }
332    
333      float GetChi2() {
334        float chi2 = 0;
335        int nchi = 0;
336        for (Int_t i = 0; i < (int) (y.size()); i++) {
337          chi2 += w[i] * w[i] * (y[i] - par[0] - par[1] * x[i]) * (y[i] - par[0] - par[1] * x[i]);
338          nchi += (int) w[i];
339        };
340    
341        if (nchi < 3)
342          return -999; // original
343          //        if(nchi<3 || chi2==0)return -999; // modified by Sergio
344    
345        chi2 = chi2 / (nchi - 2);
346        return chi2;
347      }
348    
349      float GetQaxis() {
350        float qaxis = 0;
351        for (Int_t i = 0; i < (int) (q.size()); i++)
352          qaxis += q[i];
353        //  for(Int_t i=0; i<(int)(q.size()); i++) cout<<q[i]<<endl;
354        return qaxis;
355      }
356    
357      float GetQaxis(float toll) {
358        float qaxis = 0;
359        for (int ih = 0; ih < cevent->GetN(); ih++) {
360          float x = cevent->xycoord[ih];
361          float z = cevent->zcoord[ih];
362          float q = cevent->q[ih];
363          //int ip = cevent->idp[ih];
364          float d = fabs(x - par[0] - par[1] * z) / sqrt(par[1] * par[1] + 1);
365          if (d < toll)
366            qaxis += q;
367        }
368        return qaxis;
369      }
370    
371      float GetCOG(Int_t ip) {
372        if (ip >= 0 && ip < 22)
373          return cog[ip];
374        else
375          return 0;
376      }
377    
378      float GetQ(Int_t ip) {
379        if (ip >= 0 && ip < 22)
380          return qpl[ip];
381        else
382          return 0;
383      }
384    
385      float GetQstrip() {
386    
387        float qaxis = 0;
388        for (Int_t i = 0; i < (int) (q.size()); i++)
389          qaxis += q[i];
390        if (q.size() == 0)
391          return 0.;
392        return qaxis / q.size();
393      }
394    
395      float GetYfit(float x) {
396        return par[0] + x * par[1];
397      }
398    
399      float GetXmin() {
400        float zmin = 9999;
401        for (Int_t i = 0; i < (int) (x.size()); i++)
402          if (x[i] < zmin)
403            zmin = x[i];
404        return zmin;
405      }
406      float GetXmax() {
407        float zmax = -9999;
408        for (Int_t i = 0; i < (int) (x.size()); i++)
409          if (x[i] > zmax)
410            zmax = x[i];
411        return zmax;
412    
413      }
414    
415      void Draw() {
416    
417        //  cout << " CaloAxis::Draw()"<<endl;
418        if (GetN() < 3)
419          return;
420    
421        sel = new TPolyMarker(GetN(), GetY(), GetX());
422        sel->SetMarkerSize(0.3);
423        sel->SetMarkerColor(5);
424    
425        axis = new TLine(par[0] + GetXmin() * par[1], GetXmin(), par[0] + GetXmax() * par[1], GetXmax());
426        axis->SetLineWidth(1);
427        axis->SetLineColor(5);
428    
429        //  cout << sel << endl;
430        //  cout << axis << endl;
431    
432        sel->Draw();
433        axis->Draw();
434      }
435      ;
436    
437      //-----------------------------------------------------------------
438      // fit the axis (optimized for non-interacting patterns)
439      //-----------------------------------------------------------------
440      int FitAxis(CaloLevel1* calo, Int_t view, Float_t toll, Bool_t usemechanicalalignement) {
441    
442        //  cout << "CaloAxis::FitAxis(...)"<<endl;
443    
444        Set(calo, view, usemechanicalalignement);
445    
446        // ------------------------------
447        // first : all hits included, w=1
448        // ------------------------------
449        for (int ih = 0; ih < cevent->GetN(); ih++) {
450          float x = cevent->xycoord[ih];
451          float z = cevent->zcoord[ih];
452          float q = cevent->q[ih];
453          Add(z, x, q);
454        }
455    
456        if (GetN() < 3)
457          return 0;
458        Fit();
459        // useful for debug cout << "(1)"<<endl;
460        //  cout << " n. hits :"<<GetN()<<endl;
461        //  cout << " P0 P1   :"<<par[0]<< " "<<par[1]<<endl;
462        //  cout << " Chi2    :"<<GetChi2()<<endl;
463        //  cout << " Q axis  :"<<GetQaxis()<<endl;
464        //  cout << " Q strip :"<<GetQstrip()<<endl;
465    
466        float P0;
467        float P1;
468        int dof;
469    
470        // ---------------------------
471        // second : iteration
472        // ---------------------------
473        int niter = 0; // number of iterations
474        int pok[22];
475    
476        // prova
477        // pesi lineari
478        float W22 = 1.;
479        float W1 = 10.;
480        float b = (W22 - W1) / (22. - 1.);// (w22 - w1)/(22-1)
481        float a = W1 - b;
482    
483        for (int i = 0; i < 22; i++)
484          pok[i] = 0;
485        do {
486          niter++;
487          float ttoll = toll; //tolerance (cm)
488          if (niter == 1)
489            ttoll = 10. * toll;
490          dof = GetN(); //previous fit
491          P0 = par[0];
492          P1 = par[1];
493          Reset(); //fit reset
494          for (int ih = 0; ih < cevent->GetN(); ih++) {//loop over selected hits
495            float x = cevent->xycoord[ih];
496            float z = cevent->zcoord[ih];
497            float q = cevent->q[ih];
498            int ip = cevent->idp[ih];
499            float d = fabs(x - P0 - P1 * z) / sqrt(P1 * P1 + 1);//distance to axis
500            //              cout<<d<<endl;
501            if (d < ttoll && (niter == 1 || (niter > 1 && pok[ip] == 1))) {
502              //            if( d < ttoll ){
503              float W = a + b * (ip + 1);
504              //                cout << ip << " "<<W<<endl;
505              Add(z, x, q, W);
506              pok[ip] = 1;
507            }
508          }
509          // break the track if more than 3 planes are missing
510          int nmissing = 0;
511          for (int ip = 0; ip < 22; ip++) {
512            if (pok[ip] == 0) {
513              nmissing++;
514            }
515            else {
516              nmissing = 0;
517            }
518            if (nmissing == 3) {
519              for (int ipp = ip + 1; ipp < 22; ipp++)
520                pok[ipp] = 0;
521              break;
522            }
523          }
524    
525          if (niter == 100)
526            break;
527          if (GetN() < 3)
528            break;
529          Fit();
530        } while (niter == 1 || GetN() != dof);
531        //  cout << " >>> "<<GetN()<<endl;
532        if (GetN() < 3)
533          return 0;
534        Fit();
535        //  cout << "(2)"<<endl;
536        //  cout << " n. hits :"<<GetN()<<endl;
537        //  cout << " P0 P1   :"<<par[0]<< " "<<par[1]<<endl;
538        //  cout << " Chi2    :"<<GetChi2()<<endl;
539        //  cout << " Q axis  :"<<GetQaxis()<<endl;
540        //  cout << " Q strip :"<<GetQstrip()<<endl;
541    
542        // ---------------------------------------------
543        // third : selecting only closest strip-clusters
544        // and fit baricenters
545        // ---------------------------------------------
546        P0 = par[0];
547        P1 = par[1];
548        Reset();
549    
550        float dmin[22];
551        int closest[22];
552        for (int ip = 0; ip < 22; ip++) {
553          dmin[ip] = 999;
554          closest[ip] = -1;
555        }
556        for (int ih = 0; ih < cevent->GetN(); ih++) {
557          float x = cevent->xycoord[ih];
558          float z = cevent->zcoord[ih];
559          //float q = cevent->q[ih];
560          int ip = cevent->idp[ih];
561          //int is = cevent->ids[ih];
562          float d = fabs(x - P0 - P1 * z) / sqrt(P1 * P1 + 1);
563          if (d < toll && d < dmin[ip] && pok[ip] == 1) {
564            //              Add(z,x,q,q);
565            closest[ip] = ih;
566            //              cog[ip] += x*q;
567            //              qpl[ip] += q;
568          }
569        }
570        for (Int_t ip = 0; ip < 22; ip++) {
571          if (closest[ip] != -1) {
572            float x = cevent->xycoord[closest[ip]];
573            float z = cevent->zcoord[closest[ip]];
574            float q = cevent->q[closest[ip]];
575            //int is = cevent->ids[closest[ip]];
576            Add(z, x, q, q);
577            cog[ip] += x * q;
578            qpl[ip] += q;
579            //              cout << ip << " -o- "<<is<<endl;
580          }
581        }
582        // -----------------
583        // add +/- one strip
584        // -----------------
585        for (int ih = 0; ih < cevent->GetN(); ih++) {
586          float x = cevent->xycoord[ih];
587          float z = cevent->zcoord[ih];
588          float q = cevent->q[ih];
589          int ip = cevent->idp[ih];
590          int is = cevent->ids[ih];
591          if (closest[ip] != -1) {
592            int isc = cevent->ids[closest[ip]];
593            if (is == isc + 1 || is == isc - 1) {
594              Add(z, x, q, q);
595              cog[ip] += x * q;
596              qpl[ip] += q;
597              //                cout << ip << " -+- "<<is<<endl;
598            }
599          }
600        }
601        Fit();
602    
603        if (GetN() < 3)
604          return 0;
605        for (int ip = 0; ip < 22; ip++) {
606          if (qpl[ip] != 0) {
607            cog[ip] = cog[ip] / qpl[ip];
608            //          Add(z,cog[ip],cog[ip],0.7);
609          }
610        }
611    
612        //  cout << "(3)"<<endl;
613        //          cout << " n. hits :"<<GetN()<<endl;
614        //  cout << " P0 P1   :"<<par[0]<< " "<<par[1]<<endl;
615        //          cout << " Chi2    :"<<GetChi2()<<endl;
616        //  cout << " Q axis  :"<<GetQaxis()<<endl;
617        //          cout << " Q strip :"<<GetQstrip()<<endl;
618        //  cout << " ---------------------------"<<endl;
619        return 1;
620      }
621      ;
622    
623      int FitAxis(CaloLevel1* calo, Int_t view, Float_t toll) {
624        return FitAxis(calo, view, toll, 0);
625      }
626      ;
627    
628      //-----------------------------------------------------------------
629      // fit the axis (optimized for interacting patterns?? ...forse)
630      //-----------------------------------------------------------------
631      int FitShower(CaloLevel1* calo, Int_t view, Float_t toll, Bool_t usemechanicalalignement) {
632    
633        //  cout << "CaloAxis::FitAxis(...)"<<endl;
634    
635        Set(calo, view, usemechanicalalignement);
636    
637        // ------------------------------
638        // first :
639        // ------------------------------
640        for (int ih = 0; ih < cevent->GetN(); ih++) {
641          float x = cevent->xycoord[ih];
642          //float z = cevent->zcoord[ih];
643          float q = cevent->q[ih];
644          int ip = cevent->idp[ih];
645          cog[ip] += x * q;
646          qpl[ip] += q;
647          //            Add(z,x,q);
648        }
649        for (int ip = 0; ip < 22; ip++) {
650          if (qpl[ip] != 0)
651            cog[ip] = cog[ip] / qpl[ip];
652          //            cout << ip << " > "<<qpl[ip]<<" "<<cog[ip]<<endl;
653        }
654        // ----------------
655        // plane of maximum
656        // ----------------
657        float qplmax = 1;
658        int ipmax = -1;
659        float dmin[22];
660        int closest[22];
661        for (int ip = 0; ip < 22; ip++)
662          if (qpl[ip] > qplmax) {
663            ipmax = ip;
664            qplmax = qpl[ip];
665            dmin[ip] = 1000.;//init
666            closest[ip] = -1;//init
667          }
668        for (int ih = 0; ih < cevent->GetN(); ih++) {
669          float x = cevent->xycoord[ih];
670          float z = cevent->zcoord[ih];
671          float q = cevent->q[ih];
672          int ip = cevent->idp[ih];
673          //            if( (ipmax>3 && ip<=ipmax) || ipmax<=3 )Add(z,x,q);
674          if (ip <= ipmax || ip <= 3) {
675            Add(z, x, q);
676          }
677        }
678        if (GetN() < 3)
679          return 0;
680        Fit();
681        //  cout << " n. hits :"<<GetN()<<endl;
682        //  cout << " P0 P1   :"<<par[0]<< " "<<par[1]<<endl;
683        //  cout << " Chi2    :"<<GetChi2()<<endl;
684        //  cout << " Q axis  :"<<GetQaxis()<<endl;
685        //  cout << " Q strip :"<<GetQstrip()<<endl;
686    
687        float P0;
688        float P1;
689        int dof;
690    
691        // ---------------------------
692        // second : iteration
693        // ---------------------------
694        int niter = 0; // number of iterations
695        int pok[22];
696        // prova
697        // pesi lineari
698        float W22 = 1.;
699        float W1 = 10.;
700        float b = (W22 - W1) / (22. - 1.);// (w22 - w1)/(22-1)
701        float a = W1 - b;
702    
703        for (int i = 0; i < 22; i++) {
704          pok[i] = 0;
705        }
706        do {
707          niter++;
708          float ttoll = toll; //tolerance (cm)
709          if (niter == 1)
710            ttoll = 10. * toll;
711          dof = GetN(); //previous fit
712          P0 = par[0];
713          P1 = par[1];
714          Reset(); //fit reset
715          for (int i = 0; i < 22; i++) {
716            cog[i] = 0;
717            qpl[i] = 0;
718          }
719          for (int ih = 0; ih < cevent->GetN(); ih++) {
720            float x = cevent->xycoord[ih];
721            float z = cevent->zcoord[ih];
722            float q = cevent->q[ih];
723            int ip = cevent->idp[ih];
724            float d = fabs(x - P0 - P1 * z) / sqrt(P1 * P1 + 1);
725            //              cout<<d<<endl;
726            if ((niter == 1 || (niter > 1 && pok[ip] == 1)) && (ip <= ipmax || ip <= 3) &&
727            //                  ((ipmax>3 && ip<=ipmax) || ipmax<=3) &&
728                true) {
729    
730              if (niter == 1 && d < dmin[ip]) {
731                dmin[ip] = d;
732                closest[ip] = ih;
733              }
734    
735              float W = a + b * (ip + 1);
736              //                cout << ip << " "<<W<<endl;
737              if (d < ttoll || (niter == 2 && ih == closest[ip])) {
738                /*                          if(ip==ipmax && ipmax>1)Add(z,x,q,W*q);  */
739                /*                          else Add(z,x,q,W);  */
740                cog[ip] += x * q;
741                qpl[ip] += q;
742                Add(z, x, q, W);
743                pok[ip] = 1;
744                cout << ip << " " << qpl[ip] << endl;
745              }
746            }
747          }
748          // break the track if more than 3 planes are missing
749          int nmissing = 0;
750          for (int ip = 0; ip < 22; ip++) {
751            if (pok[ip] == 0) {
752              nmissing++;
753            }
754            else {
755              nmissing = 0;
756            }
757            if (nmissing == 6) {
758              for (int ipp = ip + 1; ipp < 22; ipp++)
759                pok[ipp] = 0;
760              break;
761            }
762          }
763    
764          for (int ip = 0; ip < 22; ip++)
765            cout << qpl[ip] << " ";
766          cout << endl;
767    
768          for (int ip = 0; ip < 22; ip++)
769            if (qpl[ip] != 0)
770              cog[ip] = cog[ip] / qpl[ip];
771    
772          if (niter == 100)
773            break;
774          if (GetN() < 3)
775            break;
776    
777          Fit();
778    
779        } while (niter == 1 || GetN() != dof);
780        //  cout << " >>> "<<GetN()<<endl;
781        if (GetN() < 3)
782          return 0;
783        Fit();
784        /* //       cout << " n. hits :"<<GetN()<<endl; */
785        /* //       cout << " P0 P1   :"<<par[0]<< " "<<par[1]<<endl; */
786        /* //       cout << " Chi2    :"<<GetChi2()<<endl; */
787        /* //       cout << " Q axis  :"<<GetQaxis()<<endl; */
788        /* //       cout << " Q strip :"<<GetQstrip()<<endl; */
789    
790        /*  // --------------------------------------------- */
791        /*  // third : selecting only closest strip-clusters */
792        /*  // and fit baricenters */
793        /*  // --------------------------------------------- */
794        /*  P0  = par[0]; */
795        /*  P1  = par[1]; */
796        /*  Reset(); */
797    
798        /*  float dmin[22]; */
799        /*  int closest[22]; */
800        /*  for(int ip=0; ip<22; ip++){ */
801        /*      dmin[ip]=999; */
802        /*      closest[ip]=-1; */
803        /*  } */
804        /*  for(int ih=0;ih<cevent->GetN(); ih++){ */
805        /*      float x = cevent->xycoord[ih]; */
806        /*      float z = cevent->zcoord[ih]; */
807        /*      float q = cevent->q[ih]; */
808        /*      int ip = cevent->idp[ih]; */
809        /*      int is = cevent->ids[ih]; */
810        /*      float d = fabs(x-P0-P1*z)/sqrt(P1*P1+1); */
811        /*      if( d < toll && d<dmin[ip] && pok[ip]==1 ){ */
812        /* //               Add(z,x,q,q); */
813        /*          closest[ip]=ih; */
814        /* //               cog[ip] += x*q; */
815        /* //               qpl[ip] += q; */
816        /*      } */
817        /*  } */
818        /*  for(Int_t ip=0; ip<22; ip++){ */
819        /*      if(closest[ip]!=-1){ */
820        /*          float x = cevent->xycoord[closest[ip]]; */
821        /*          float z = cevent->zcoord[closest[ip]]; */
822        /*          float q = cevent->q[closest[ip]]; */
823        /*          int is = cevent->ids[closest[ip]]; */
824        /*          Add(z,x,q,q); */
825        /*          cog[ip] += x*q; */
826        /*          qpl[ip] += q; */
827        /* //               cout << ip << " -o- "<<is<<endl; */
828        /*      } */
829        /*  } */
830        /*  // add +/- one strip */
831        /*  for(int ih=0;ih<cevent->GetN(); ih++){ */
832        /*      float x = cevent->xycoord[ih]; */
833        /*      float z = cevent->zcoord[ih]; */
834        /*      float q = cevent->q[ih]; */
835        /*      int ip = cevent->idp[ih]; */
836        /*      int is = cevent->ids[ih]; */
837        /*      if( closest[ip]!=-1 ){ */
838        /*          int isc = cevent->ids[closest[ip]]; */
839        /*          if( is == isc+1 || is == isc-1 ){ */
840        /*              Add(z,x,q,q); */
841        /*              cog[ip] += x*q; */
842        /*              qpl[ip] += q; */
843        /* //                   cout << ip << " -+- "<<is<<endl; */
844        /*          } */
845        /*      } */
846        /*  } */
847        /*  Fit(); */
848    
849        /*  if(GetN()<3)return 0; */
850        /*  for(int ip=0; ip<22; ip++){ */
851        /*      if(qpl[ip]!=0){ */
852        /*          cog[ip]=cog[ip]/qpl[ip]; */
853        /* //           Add(z,cog[ip],cog[ip],0.7); */
854        /*      } */
855        /*  } */
856    
857        //  cout << " n. hits :"<<GetN()<<endl;
858        //  cout << " P0 P1   :"<<par[0]<< " "<<par[1]<<endl;
859        //  cout << " Chi2    :"<<GetChi2()<<endl;
860        //  cout << " Q axis  :"<<GetQaxis()<<endl;
861        //  cout << " Q strip :"<<GetQstrip()<<endl;
862        //  cout << " ---------------------------"<<endl;
863        return 1;
864      }
865      ;
866    
867      int FitShower(CaloLevel1* calo, Int_t view, Float_t toll) {
868        return FitAxis(calo, view, toll, 0);
869      }
870      ;
871    
872  };  };
873    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.23