/[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.2 by pam-fi, Wed Aug 5 14:07:04 2009 UTC revision 1.3 by pam-fi, Fri Sep 25 15:34:40 2009 UTC
# Line 340  public: Line 340  public:
340    
341      if (nchi < 3)      if (nchi < 3)
342        return -999; // original        return -999; // original
343        //        if(nchi<3 || chi2==0)return -999; // modified by Sergio      //  if(nchi<3 || chi2==0)return -999; // modified by Sergio
344    
345      chi2 = chi2 / (nchi - 2);      chi2 = chi2 / (nchi - 2);
346      return chi2;      return chi2;
# Line 437  public: Line 437  public:
437    //-----------------------------------------------------------------    //-----------------------------------------------------------------
438    // fit the axis (optimized for non-interacting patterns)    // fit the axis (optimized for non-interacting patterns)
439    //-----------------------------------------------------------------    //-----------------------------------------------------------------
440    int FitAxis(CaloLevel1* calo, Int_t view, Float_t toll, Bool_t usemechanicalalignement) {    int FitAxis(CaloLevel1* calo, Int_t view, Float_t toll, Int_t nPlanes = 22, Bool_t usemechanicalalignement = 0) {
441    
442      //  cout << "CaloAxis::FitAxis(...)"<<endl;      //  cout << "CaloAxis::FitAxis(...)"<<endl;
443    
# Line 471  public: Line 471  public:
471      // second : iteration      // second : iteration
472      // ---------------------------      // ---------------------------
473      int niter = 0; // number of iterations      int niter = 0; // number of iterations
474      int pok[22];      int pok[nPlanes];
475    
476      // prova      // prova
477      // pesi lineari      // pesi lineari
478      float W22 = 1.;      float W22 = 1.;
479      float W1 = 10.;      float W1 = 10.;
480      float b = (W22 - W1) / (22. - 1.);// (w22 - w1)/(22-1)      float b = (W22 - W1) / ((float) nPlanes - 1.);// (w22 - w1)/(22-1)
481      float a = W1 - b;      float a = W1 - b;
482    
483      for (int i = 0; i < 22; i++)      for (int i = 0; i < nPlanes; i++)
484        pok[i] = 0;        pok[i] = 0;
485      do {      do {
486        niter++;        niter++;
# Line 492  public: Line 492  public:
492        P1 = par[1];        P1 = par[1];
493        Reset(); //fit reset        Reset(); //fit reset
494        for (int ih = 0; ih < cevent->GetN(); ih++) {//loop over selected hits        for (int ih = 0; ih < cevent->GetN(); ih++) {//loop over selected hits
495          float x = cevent->xycoord[ih];          if (cevent->idp[ih] < nPlanes) {
496          float z = cevent->zcoord[ih];            float x = cevent->xycoord[ih];
497          float q = cevent->q[ih];            float z = cevent->zcoord[ih];
498          int ip = cevent->idp[ih];            float q = cevent->q[ih];
499          float d = fabs(x - P0 - P1 * z) / sqrt(P1 * P1 + 1);//distance to axis            int ip = cevent->idp[ih];
500          //              cout<<d<<endl;            float d = fabs(x - P0 - P1 * z) / sqrt(P1 * P1 + 1);//distance to axis
501          if (d < ttoll && (niter == 1 || (niter > 1 && pok[ip] == 1))) {            //            cout<<d<<endl;
502            //            if( d < ttoll ){            if (d < ttoll && (niter == 1 || (niter > 1 && pok[ip] == 1))) {
503            float W = a + b * (ip + 1);              //          if( d < ttoll ){
504            //                cout << ip << " "<<W<<endl;              float W = a + b * (ip + 1);
505            Add(z, x, q, W);              //              cout << ip << " "<<W<<endl;
506            pok[ip] = 1;              Add(z, x, q, W);
507                pok[ip] = 1;
508              }
509          }          }
510        }        }
511        // break the track if more than 3 planes are missing        // break the track if more than 3 planes are missing
512        int nmissing = 0;        int nmissing = 0;
513        for (int ip = 0; ip < 22; ip++) {        for (int ip = 0; ip < nPlanes; ip++) {
514          if (pok[ip] == 0) {          if (pok[ip] == 0) {
515            nmissing++;            nmissing++;
516          }          }
# Line 516  public: Line 518  public:
518            nmissing = 0;            nmissing = 0;
519          }          }
520          if (nmissing == 3) {          if (nmissing == 3) {
521            for (int ipp = ip + 1; ipp < 22; ipp++)            for (int ipp = ip + 1; ipp < nPlanes; ipp++)
522              pok[ipp] = 0;              pok[ipp] = 0;
523            break;            break;
524          }          }
# Line 547  public: Line 549  public:
549      P1 = par[1];      P1 = par[1];
550      Reset();      Reset();
551    
552      float dmin[22];      float dmin[nPlanes];
553      int closest[22];      int closest[nPlanes];
554      for (int ip = 0; ip < 22; ip++) {      for (int ip = 0; ip < nPlanes; ip++) {
555        dmin[ip] = 999;        dmin[ip] = 999;
556        closest[ip] = -1;        closest[ip] = -1;
557      }      }
558      for (int ih = 0; ih < cevent->GetN(); ih++) {      for (int ih = 0; ih < cevent->GetN(); ih++) {
559        float x = cevent->xycoord[ih];        if (cevent->idp[ih] < nPlanes) {
560        float z = cevent->zcoord[ih];          float x = cevent->xycoord[ih];
561        //float q = cevent->q[ih];          float z = cevent->zcoord[ih];
562        int ip = cevent->idp[ih];          //float q = cevent->q[ih];
563        //int is = cevent->ids[ih];          int ip = cevent->idp[ih];
564        float d = fabs(x - P0 - P1 * z) / sqrt(P1 * P1 + 1);          //int is = cevent->ids[ih];
565        if (d < toll && d < dmin[ip] && pok[ip] == 1) {          float d = fabs(x - P0 - P1 * z) / sqrt(P1 * P1 + 1);
566          //              Add(z,x,q,q);          if (d < toll && d < dmin[ip] && pok[ip] == 1) {
567          closest[ip] = ih;            //            Add(z,x,q,q);
568          //              cog[ip] += x*q;            closest[ip] = ih;
569          //              qpl[ip] += q;            //            cog[ip] += x*q;
570              //            qpl[ip] += q;
571            }
572        }        }
573      }      }
574      for (Int_t ip = 0; ip < 22; ip++) {      for (Int_t ip = 0; ip < nPlanes; ip++) {
575        if (closest[ip] != -1) {        if (closest[ip] != -1) {
576          float x = cevent->xycoord[closest[ip]];          float x = cevent->xycoord[closest[ip]];
577          float z = cevent->zcoord[closest[ip]];          float z = cevent->zcoord[closest[ip]];
# Line 583  public: Line 587  public:
587      // add +/- one strip      // add +/- one strip
588      // -----------------      // -----------------
589      for (int ih = 0; ih < cevent->GetN(); ih++) {      for (int ih = 0; ih < cevent->GetN(); ih++) {
590          if (cevent->idp[ih] < nPlanes) {
591        float x = cevent->xycoord[ih];        float x = cevent->xycoord[ih];
592        float z = cevent->zcoord[ih];        float z = cevent->zcoord[ih];
593        float q = cevent->q[ih];        float q = cevent->q[ih];
# Line 596  public: Line 601  public:
601            qpl[ip] += q;            qpl[ip] += q;
602            //                cout << ip << " -+- "<<is<<endl;            //                cout << ip << " -+- "<<is<<endl;
603          }          }
604        }        }}
605      }      }
606      Fit();      Fit();
607    
608      if (GetN() < 3)      if (GetN() < 3)
609        return 0;        return 0;
610      for (int ip = 0; ip < 22; ip++) {      for (int ip = 0; ip < nPlanes; ip++) {
611        if (qpl[ip] != 0) {        if (qpl[ip] != 0) {
612          cog[ip] = cog[ip] / qpl[ip];          cog[ip] = cog[ip] / qpl[ip];
613          //          Add(z,cog[ip],cog[ip],0.7);          //          Add(z,cog[ip],cog[ip],0.7);
# Line 620  public: Line 625  public:
625    }    }
626    ;    ;
627    
628    int FitAxis(CaloLevel1* calo, Int_t view, Float_t toll) {    /*int FitAxis(CaloLevel1* calo, Int_t view, Float_t toll, Int_t nPlanes = 22) {
629      return FitAxis(calo, view, toll, 0);      return FitAxis(calo, view, toll, 0, nPlanes);
630    }    }
631    ;    ;*/
632    
633    //-----------------------------------------------------------------    //-----------------------------------------------------------------
634    // fit the axis (optimized for interacting patterns?? ...forse)    // fit the axis (optimized for interacting patterns?? ...forse)
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.23