/[PAMELA software]/calo/flight/CaloEnergy/src/CaloEnergy.cpp
ViewVC logotype

Diff of /calo/flight/CaloEnergy/src/CaloEnergy.cpp

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

revision 1.4 by mocchiut, Wed Jul 29 12:58:27 2009 UTC revision 1.6 by mocchiut, Tue Aug 11 13:32:04 2009 UTC
# Line 39  void CaloEnergy::UseCaloPreSampler(){ Line 39  void CaloEnergy::UseCaloPreSampler(){
39    //    //
40    // use the presampler setting forcefitmode to 1000 means to force the DV routine to find the track inside the calorimeter using the "shower" approach developed for electrons    // use the presampler setting forcefitmode to 1000 means to force the DV routine to find the track inside the calorimeter using the "shower" approach developed for electrons
41    //    //
42    cp = new CaloPreSampler(L2);    if ( !cp ) cp = new CaloPreSampler(L2);
43    cp->SplitInto(0,22);    cp->SplitInto(0,22);
44    cp->SetForceFitMode(1000);    cp->SetForceFitMode(1000);
45  //   cp->UseTracker(false);  //   cp->UseTracker(false);
# Line 48  void CaloEnergy::UseCaloPreSampler(){ Line 48  void CaloEnergy::UseCaloPreSampler(){
48  //   cp->Process();  //   cp->Process();
49  }  }
50    
51    
52    void CaloEnergy::UseLongFit(){
53      fPl = 0;
54      fLong = true;
55      if ( !clong ){
56        clong = new CaloLong(L2);
57        if ( cp ) clong->SetCaloLevel2Pointer(cp->GetLevel2Pointer());
58        clong->SplitInto(0,22);
59      };
60      //
61    }
62    
63  void CaloEnergy::Set(){  void CaloEnergy::Set(){
64    //    //
65    // set default values, NB default conversion factor for energy is just very approximated!    // set default values, NB default conversion factor for energy is just very approximated!
# Line 70  void CaloEnergy::Set(){ Line 82  void CaloEnergy::Set(){
82    fPl = 1;    fPl = 1;
83    fRad = -1;    fRad = -1;
84    cp = NULL;    cp = NULL;
85      clong = NULL;
86      x0max = -1.;
87    //    //
88    this->DefineGeometry();    this->DefineGeometry();
89    fXosel =true;    fXosel =true;
# Line 139  Int_t CaloEnergy::GetMaxplane(TString se Line 153  Int_t CaloEnergy::GetMaxplane(TString se
153    return(-1);    return(-1);
154  }  }
155    
156    Float_t CaloEnergy::GetEnergyAtMaxplane(TString section){
157      section.ToUpper();
158      if ( section.Contains("XO") ) return xomax_en;
159      if ( section.Contains("XE") ) return xemax_en;
160      if ( section.Contains("YO") ) return yomax_en;
161      if ( section.Contains("YE") ) return yemax_en;
162      return(-1);
163    }
164    
165  Float_t CaloEnergy::GetMaxEnergy(TString section){  Float_t CaloEnergy::GetMaxEnergy(TString section){
166    section.ToUpper();    section.ToUpper();
167    if ( section.Contains("XO") ) return fXOen_maxplane;    if ( fLong ){
168    if ( section.Contains("XE") ) return fXEen_maxplane;      this->Process(section);    
169    if ( section.Contains("YO") ) return fYOen_maxplane;      return fXOen_maxplane;
170    if ( section.Contains("YE") ) return fYEen_maxplane;    } else {
171        if ( section.Contains("XO") ) return fXOen_maxplane;
172        if ( section.Contains("XE") ) return fXEen_maxplane;
173        if ( section.Contains("YO") ) return fYOen_maxplane;
174        if ( section.Contains("YE") ) return fYEen_maxplane;
175      };
176    return(-1);    return(-1);
177  }  }
178    
179    Float_t CaloEnergy::GetMaxEnergy(){  
180      if ( fLong ){
181        if ( debug ) printf(" oh! call process! with asntr %s and sntr %s \n",asntr.Data(),sntr.Data());
182        this->Process(asntr);
183      };
184      return((fXEen_maxplane+fYOen_maxplane+fYEen_maxplane+fXOen_maxplane));
185    }
186    
187  void CaloEnergy::DefineGeometry(){  void CaloEnergy::DefineGeometry(){
188    //    //
189    // Use CaloStrip to determine once the position of border strips for each section    // Use CaloStrip to determine once the position of border strips for each section
# Line 272  void CaloEnergy::Clear(){ Line 308  void CaloEnergy::Clear(){
308    fMax_planexo = 0;    fMax_planexo = 0;
309    fMax_planexe = 0;    fMax_planexe = 0;
310    fMax_planeyo = 0;    fMax_planeyo = 0;
311    fMax_planeye = 0;    fMax_planeye = 0;
312    memset(enstrip,0,2*22*96*(sizeof(Float_t)));      xomax_en= 0.;
313      xemax_en= 0.;
314      yomax_en= 0.;
315      yemax_en= 0.;
316      //
317     memset(enstrip,0,2*22*96*(sizeof(Float_t)));  
318    en = 0.;    en = 0.;
319    view = 0;    view = 0;
320    plane = 0;    plane = 0;
# Line 336  void CaloEnergy::Print(){ Line 377  void CaloEnergy::Print(){
377    printf(" energyxo:.............. %f \n",energyxo);    printf(" energyxo:.............. %f \n",energyxo);
378    printf(" energyye:.............. %f \n",energyye);    printf(" energyye:.............. %f \n",energyye);
379    printf(" energyyo:.............. %f \n",energyyo);    printf(" energyyo:.............. %f \n",energyyo);
380      printf(" x0max   :.............. %f \n",x0max);
381    printf(" debug   :.............. %i \n",debug);    printf(" debug   :.............. %i \n",debug);
382    
383    printf("========================================================================\n");    printf("========================================================================\n");
# Line 421  Bool_t CaloEnergy::IsInsideAcceptance(TS Line 463  Bool_t CaloEnergy::IsInsideAcceptance(TS
463          //          //
464          en_xep[i] += enstrip[1][2*i][strip];          en_xep[i] += enstrip[1][2*i][strip];
465          en_yop[i] += enstrip[0][2*i][strip];          en_yop[i] += enstrip[0][2*i][strip];
466          en_xop[i] += enstrip[1][2*i+1][strip];          en_xop[i] += enstrip[1][(2*i)+1][strip];
467          en_yep[i] += enstrip[0][2*i+1][strip];          en_yep[i] += enstrip[0][(2*i)+1][strip];
468        } else {        } else {
469          //          //
470          // use only the strips inside a cylinder of given radius fRad          // use only the strips inside a cylinder of given radius fRad
471          //          //
472          if ( strip >= cl2->cibar[2*i][1]-1-fRad &&  strip <= cl2->cibar[2*i][1]-1+fRad ) en_xep[i] += enstrip[1][2*i][strip];          if ( strip >= cl2->cibar[2*i][1]-1-fRad &&  strip <= cl2->cibar[2*i][1]-1+fRad ) en_xep[i] += enstrip[1][2*i][strip];
473          if ( strip >= cl2->cibar[2*i][0]-1-fRad &&  strip <= cl2->cibar[2*i][0]-1+fRad ) en_yop[i] += enstrip[0][2*i][strip];          if ( strip >= cl2->cibar[2*i][0]-1-fRad &&  strip <= cl2->cibar[2*i][0]-1+fRad ) en_yop[i] += enstrip[0][2*i][strip];
474          if ( strip >= cl2->cibar[2*i+1][1]-1-fRad &&  strip <= cl2->cibar[2*i+1][1]-1+fRad ) en_xop[i] += enstrip[1][2*i+1][strip];          if ( strip >= cl2->cibar[(2*i)+1][1]-1-fRad &&  strip <= cl2->cibar[(2*i)+1][1]-1+fRad ) en_xop[i] += enstrip[1][(2*i)+1][strip];
475          if ( strip >= cl2->cibar[2*i+1][0]-1-fRad &&  strip <= cl2->cibar[2*i+1][0]-1+fRad ) en_yep[i] += enstrip[0][2*i+1][strip];          if ( strip >= cl2->cibar[(2*i)+1][0]-1-fRad &&  strip <= cl2->cibar[(2*i)+1][0]-1+fRad ) en_yep[i] += enstrip[0][(2*i)+1][strip];
476        };        };
477      };      };
478      energyxe += en_xep[i];      energyxe += en_xep[i];
# Line 441  Bool_t CaloEnergy::IsInsideAcceptance(TS Line 483  Bool_t CaloEnergy::IsInsideAcceptance(TS
483    //    //
484    // Find the plane of maximum for each section    // Find the plane of maximum for each section
485    //    //
   Float_t xomax_en= 0.;  
   Float_t xemax_en= 0.;  
   Float_t yomax_en= 0.;  
   Float_t yemax_en= 0.;  
   //  
486    //    //
487    Int_t xen = 0;    Int_t xen = 0;
488    Int_t yon = 0;    Int_t yon = 0;
# Line 494  Bool_t CaloEnergy::IsInsideAcceptance(TS Line 531  Bool_t CaloEnergy::IsInsideAcceptance(TS
531      };      };
532    };    };
533    //    //
534      // the maximum is given externally as X0, convert it to plane and section
535      //
536      if ( x0max > 0. ){
537        if ( debug ) printf(" CALCULATE MAX PLANE GIVEN X0 ASSUMING PERPENDICULAR TRACK \n");
538        Int_t wpl = (Int_t)roundf(x0max/0.76);
539        Bool_t isY = false;
540        if ( ((x0max/0.76)-(Float_t)wpl) > 0. ) isY = true;
541        if ( !(wpl%2) ){
542          // 0, 2, 4, ...
543          if ( isY ){
544            yomax_en = 1000.;
545            xemax_en = 500.;
546            fMax_planeyo=wpl/2;
547            fMax_planexe=wpl/2;
548          } else {
549            yomax_en = 500.;
550            xemax_en = 1000.;
551            fMax_planeyo=wpl/2;
552            fMax_planexe=wpl/2;
553          };
554          yomax_en = 0.;
555          xemax_en = 0.;
556        } else {
557          // 1, 3, 5, ...
558          if ( isY ){
559            yemax_en = 1000.;
560            xomax_en = 500.;
561            fMax_planeye=(wpl-1)/2;
562            fMax_planexo=(wpl-1)/2;
563          } else {
564            yemax_en = 500.;
565            xomax_en = 1000.;
566            fMax_planeye=(wpl-1)/2;
567            fMax_planexo=(wpl-1)/2;
568          };
569          yemax_en = 0.;
570          xomax_en = 0.;
571        };
572        if ( debug ) printf(" x0max %f wpl %i isY %i yomax_en %f xemax_en %f yemax_en %f xomax_en %f fMaxplane %i %i %i %i\n",x0max,wpl,isY,yomax_en,xemax_en,yemax_en,xomax_en,fMax_planeyo,fMax_planexe,fMax_planeye,fMax_planexo);
573      };
574      //
575    Int_t nPl = fPl;    Int_t nPl = fPl;
576    //    //
577    // Set the maximum in case of coherent mode was selected    // Set the maximum in case of coherent mode was selected
# Line 588  Bool_t CaloEnergy::IsInsideAcceptance(TS Line 666  Bool_t CaloEnergy::IsInsideAcceptance(TS
666      for (Int_t i=0; i<11; i++) {      for (Int_t i=0; i<11; i++) {
667        if        if
668          (          (
669           (((track_coordx[2*i+1][1]>=(-12.054+fM))&&           (((track_coordx[(2*i)+1][1]>=(-12.054+fM))&&
670             (track_coordx[2*i+1][1]<=(-4.246-fM)))&&             (track_coordx[(2*i)+1][1]<=(-4.246-fM)))&&
671            (((cl2->cbar[2*i+1][1]>=xo1 + fM1)&&            (((cl2->cbar[(2*i)+1][1]>=xo1 + fM1)&&
672              (cl2->cbar[2*i+1][1]<=xo2 - fM1 ))||              (cl2->cbar[(2*i)+1][1]<=xo2 - fM1 ))||
673             ((cl2->cbar[2*i+1][1]>=xo3 + fM1)&&             ((cl2->cbar[(2*i)+1][1]>=xo3 + fM1)&&
674              (cl2->cbar[2*i+1][1]<=xo4 - fM1 ))||              (cl2->cbar[(2*i)+1][1]<=xo4 - fM1 ))||
675             ((cl2->cbar[2*i+1][1]>=xo5 + fM1)&&             ((cl2->cbar[(2*i)+1][1]>=xo5 + fM1)&&
676              (cl2->cbar[2*i+1][1]<=xo6 - fM1 ))))||              (cl2->cbar[(2*i)+1][1]<=xo6 - fM1 ))))||
677                    
678           (((track_coordx[2*i+1][1]>=(-4.004+fM))&&           (((track_coordx[(2*i)+1][1]>=(-4.004+fM))&&
679             (track_coordx[2*i+1][1]<=(3.804-fM)))&&             (track_coordx[(2*i)+1][1]<=(3.804-fM)))&&
680            (((cl2->cbar[2*i+1][1]>=xo1 + fM1)&&            (((cl2->cbar[(2*i)+1][1]>=xo1 + fM1)&&
681              (cl2->cbar[2*i+1][1]<=xo2 - fM1 ))||              (cl2->cbar[(2*i)+1][1]<=xo2 - fM1 ))||
682             ((cl2->cbar[2*i+1][1]>=xo3 + fM1)&&             ((cl2->cbar[(2*i)+1][1]>=xo3 + fM1)&&
683              (cl2->cbar[2*i+1][1]<=xo4 - fM1))||              (cl2->cbar[(2*i)+1][1]<=xo4 - fM1))||
684             ((cl2->cbar[2*i+1][1]>=xo5 + fM1)&&             ((cl2->cbar[(2*i)+1][1]>=xo5 + fM1)&&
685              (cl2->cbar[2*i+1][1]<=xo6 - fM1 ))))||              (cl2->cbar[(2*i)+1][1]<=xo6 - fM1 ))))||
686                    
687           (((track_coordx[2*i+1][1]>=(4.046+fM))&&           (((track_coordx[(2*i)+1][1]>=(4.046+fM))&&
688             (track_coordx[2*i+1][1]<=(11.854-fM)))&&             (track_coordx[(2*i)+1][1]<=(11.854-fM)))&&
689            (((cl2->cbar[2*i+1][1]>=xo1 + fM1)&&            (((cl2->cbar[(2*i)+1][1]>=xo1 + fM1)&&
690              (cl2->cbar[2*i+1][1]<=xo2 - fM1))||              (cl2->cbar[(2*i)+1][1]<=xo2 - fM1))||
691             ((cl2->cbar[2*i+1][1]>=xo3 + fM1)&&             ((cl2->cbar[(2*i)+1][1]>=xo3 + fM1)&&
692              (cl2->cbar[2*i+1][1]<=xo4 - fM1 ))||              (cl2->cbar[(2*i)+1][1]<=xo4 - fM1 ))||
693             ((cl2->cbar[2*i+1][1]>=xo5 + fM1)&&             ((cl2->cbar[(2*i)+1][1]>=xo5 + fM1)&&
694              (cl2->cbar[2*i+1][1]<=xo6 - fM1 ))))              (cl2->cbar[(2*i)+1][1]<=xo6 - fM1 ))))
695           ){           ){
696          fXosel = true;          fXosel = true;
697          fXoout = i;          fXoout = i;
# Line 692  Bool_t CaloEnergy::IsInsideAcceptance(TS Line 770  Bool_t CaloEnergy::IsInsideAcceptance(TS
770      for (Int_t i=0; i<11; i++) {      for (Int_t i=0; i<11; i++) {
771        if        if
772          (          (
773           (((track_coordy[2*i+1][0]>=(-12.154+fM))&&           (((track_coordy[(2*i)+1][0]>=(-12.154+fM))&&
774             (track_coordy[2*i+1][0]<=(-4.346-fM)))&&             (track_coordy[(2*i)+1][0]<=(-4.346-fM)))&&
775            (((cl2->cbar[2*i+1][0]>=ye1 + fM1)&&            (((cl2->cbar[(2*i)+1][0]>=ye1 + fM1)&&
776              (cl2->cbar[2*i+1][0]<=ye2 - fM1 ))||              (cl2->cbar[(2*i)+1][0]<=ye2 - fM1 ))||
777             ((cl2->cbar[2*i+1][0]>=ye3 + fM1)&&             ((cl2->cbar[(2*i)+1][0]>=ye3 + fM1)&&
778              (cl2->cbar[2*i+1][0]<=ye4 - fM1 ))||              (cl2->cbar[(2*i)+1][0]<=ye4 - fM1 ))||
779             ((cl2->cbar[2*i+1][0]>=ye5 + fM1)&&             ((cl2->cbar[(2*i)+1][0]>=ye5 + fM1)&&
780              (cl2->cbar[2*i+1][0]<=ye6 - fM1 ))))||              (cl2->cbar[(2*i)+1][0]<=ye6 - fM1 ))))||
781                    
782           (((track_coordy[2*i+1][0]>=(-4.104+fM))&&           (((track_coordy[(2*i)+1][0]>=(-4.104+fM))&&
783             (track_coordy[2*i+1][0]<=(3.704-fM)))&&             (track_coordy[(2*i)+1][0]<=(3.704-fM)))&&
784            (((cl2->cbar[2*i+1][0]>=ye1 + fM1)&&            (((cl2->cbar[(2*i)+1][0]>=ye1 + fM1)&&
785              (cl2->cbar[2*i+1][0]<=ye2 - fM1 ))||              (cl2->cbar[(2*i)+1][0]<=ye2 - fM1 ))||
786             ((cl2->cbar[2*i+1][0]>=ye3 + fM1)&&             ((cl2->cbar[(2*i)+1][0]>=ye3 + fM1)&&
787              (cl2->cbar[2*i+1][0]<=ye4 - fM1))||              (cl2->cbar[(2*i)+1][0]<=ye4 - fM1))||
788             ((cl2->cbar[2*i+1][0]>=ye5 + fM1)&&             ((cl2->cbar[(2*i)+1][0]>=ye5 + fM1)&&
789              (cl2->cbar[2*i+1][0]<=ye6 - fM1 ))))||              (cl2->cbar[(2*i)+1][0]<=ye6 - fM1 ))))||
790                    
791           (((track_coordy[2*i+1][0]>=(3.946+fM))&&           (((track_coordy[(2*i)+1][0]>=(3.946+fM))&&
792             (track_coordy[2*i+1][0]<=(11.754-fM)))&&             (track_coordy[(2*i)+1][0]<=(11.754-fM)))&&
793            (((cl2->cbar[2*i+1][0]>=ye1 + fM1)&&            (((cl2->cbar[(2*i)+1][0]>=ye1 + fM1)&&
794              (cl2->cbar[2*i+1][0]<=ye2 - fM1))||              (cl2->cbar[(2*i)+1][0]<=ye2 - fM1))||
795             ((cl2->cbar[2*i+1][0]>=ye3 + fM1)&&             ((cl2->cbar[(2*i)+1][0]>=ye3 + fM1)&&
796              (cl2->cbar[2*i+1][0]<=ye4 - fM1 ))||              (cl2->cbar[(2*i)+1][0]<=ye4 - fM1 ))||
797             ((cl2->cbar[2*i+1][0]>=ye5 + fM1)&&             ((cl2->cbar[(2*i)+1][0]>=ye5 + fM1)&&
798              (cl2->cbar[2*i+1][0]<=ye6 - fM1 ))))              (cl2->cbar[(2*i)+1][0]<=ye6 - fM1 ))))
799           ){           ){
800          fYesel = true;          fYesel = true;
801          fYeout = i;          fYeout = i;
# Line 791  Bool_t CaloEnergy::IsInsideAcceptance(TS Line 869  Bool_t CaloEnergy::IsInsideAcceptance(TS
869    // Count the number of sections in which the event is contained    // Count the number of sections in which the event is contained
870    //    //
871    fCount = (Float_t)((Int_t)fXesel+(Int_t)fXosel+(Int_t)fYesel+(Int_t)fYosel);    fCount = (Float_t)((Int_t)fXesel+(Int_t)fXosel+(Int_t)fYesel+(Int_t)fYosel);
   if ( debug ) printf(" IsInside XE %i XO %i YE %i YO %i => SEL %i \n",fXesel,fXosel,fYesel,fYosel,fSel);  
872    //    //
873    if ( indep ){    if ( indep ){
874      //      //
# Line 814  Bool_t CaloEnergy::IsInsideAcceptance(TS Line 891  Bool_t CaloEnergy::IsInsideAcceptance(TS
891    };    };
892    //    //
893    if ( debug ) printf("sel %i indep %i fMax_plane %f conv_r %f en_maxplane %f encalo %f  \n",fSel,indep,fMax_plane,fConv_rxo,fXOen_maxplane,fEnergy);    if ( debug ) printf("sel %i indep %i fMax_plane %f conv_r %f en_maxplane %f encalo %f  \n",fSel,indep,fMax_plane,fConv_rxo,fXOen_maxplane,fEnergy);
894      if ( debug ) printf(" IsInside XE %i XO %i YE %i YO %i => SEL %i \n",fXesel,fXosel,fYesel,fYosel,fSel);
895    //    //
896    // finish exit    // finish exit
897    //    //
# Line 865  void CaloEnergy::Process(TString section Line 943  void CaloEnergy::Process(TString section
943    if ( (fM1+0.122-0.244*(Float_t)fRad) < 0. ) printf("Error: (fM1+0.122-0.244*(Float_t)fRad) < 0. fM1 %f fRad %i %f \n",fM1,fRad,(fM1+0.122-0.244*(Float_t)fRad));    if ( (fM1+0.122-0.244*(Float_t)fRad) < 0. ) printf("Error: (fM1+0.122-0.244*(Float_t)fRad) < 0. fM1 %f fRad %i %f \n",fM1,fRad,(fM1+0.122-0.244*(Float_t)fRad));
944    //    //
945    if ( fLong ){    if ( fLong ){
946        if ( debug ) printf(" ==================================================================> LONGITUDINAL FIT! \n");
947      //      //
948      // use long fit to measure energy NOT IMPLEMENTED YET      // use long fit to measure energy
949        //
950        if ( this->IsInsideAcceptance(section) ){
951          //
952          if ( debug ) printf(" ==================================================================> LONG INSIDE! \n");
953          //
954          Float_t myene[2][22];
955          memset(myene,0,(sizeof(Float_t))*2*22);
956          for (Int_t j=0; j<11; j++){
957            if ( section.Contains("XE") ) myene[1][2*j] = en_xep[j];
958            if ( section.Contains("YO") ) myene[0][2*j] = en_yop[j];
959            if ( section.Contains("XO") ) myene[1][(2*j)+1] = en_xop[j];
960            if ( section.Contains("YE") ) myene[0][(2*j)+1] = en_yep[j];    
961          };
962          clong->UnMaskSections();
963          if ( !(section.Contains("YE")) ) clong->MaskSection("YE");
964          if ( !(section.Contains("YO")) ) clong->MaskSection("YO");
965          if ( !(section.Contains("XO")) ) clong->MaskSection("XO");
966          if ( !(section.Contains("XE")) ) clong->MaskSection("XE");
967          clong->ForceNextFit();
968          clong->SetEnergies(myene);
969          if ( debug ){
970            clong->Fit(true);
971          } else {
972            clong->Fit();
973          };      
974          fXOen_maxplane = clong->Get_E0();
975          fYOen_maxplane  = 0.;
976          fYEen_maxplane = 0.;
977          fXEen_maxplane = 0.;
978          fEnergy=fXOen_maxplane/fConv_rxo;
979          if ( fEnergy != fEnergy || clong->Get_fitresult() != 0 ) fEnergy = -1.;
980          //      if ( fEnergy != fEnergy ) fEnergy = 1.;
981          //
982        } else {
983          //
984          // if the event is not in the acceptance, return a negative energy.
985          //
986          if ( debug ) printf(" Outside acceptance \n");
987          fEnergy *= -1.;
988          //
989        };
990      //      //
991    } else {    } else {
992      //      //
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.23