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

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

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

revision 1.15 by mocchiut, Fri Nov 28 16:00:19 2008 UTC revision 1.20 by pam-fi, Mon Apr 20 09:12:05 2015 UTC
# Line 7  Line 7 
7  /**  /**
8   * Default constructor   * Default constructor
9   */   */
10  CaloNuclei::CaloNuclei(){  // CaloNuclei::CaloNuclei(){
11    Clear();  //   Clear();
12  };  // };
13    
14  CaloNuclei::CaloNuclei(PamLevel2 *l2p){    CaloNuclei::CaloNuclei(PamLevel2 *l2p,const char* alg){  
15    //    //
16    Clear();    Clear();
17    //    //
# Line 28  CaloNuclei::CaloNuclei(PamLevel2 *l2p){ Line 28  CaloNuclei::CaloNuclei(PamLevel2 *l2p){
28    debug = false;    debug = false;
29    // debug = true;    // debug = true;
30    usetrack = true;    usetrack = true;
31      usepl18x = false;
32      trkAlg = alg;
33    //    //
34  };  };
35    
# Line 146  void CaloNuclei::Process(Int_t ntr){ Line 148  void CaloNuclei::Process(Int_t ntr){
148      //      //
149      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
150      //      //
151        if ( !usepl18x && view==0 && plane==18 ) mip = 0.;
152        //
153        //
154      // put in vfpl vector the energy release on the first plane      // put in vfpl vector the energy release on the first plane
155      //      //
156      if ( strip != -1 && view == 1 && plane == 0 ) {      if ( strip != -1 && view == 1 && plane == 0 ) {
# Line 186  void CaloNuclei::Process(Int_t ntr){ Line 191  void CaloNuclei::Process(Int_t ntr){
191    //    //
192    if ( usetrack ){    if ( usetrack ){
193      if ( ntr >= 0 ){      if ( ntr >= 0 ){
194        ptrack = L2->GetTrack(ntr);          ptrack = L2->GetTrack(ntr,trkAlg);
195        if ( ptrack ) track = ptrack->GetCaloTrack();        if ( ptrack ) track = ptrack->GetCaloTrack();
196      } else {      } else {
197        track = L2->GetCaloStoredTrack(ntr);        track = L2->GetCaloStoredTrack(ntr);
198      };      };
199      //      //
200      if ( !track && ntr >= 0 ){      if ( !track && ntr >= 0 ){
201        printf(" ERROR: cannot find any track!\n");        printf(" ERROR: cannot find any track! \n");      
202          cout << "ERROR: trk.algorythm --> "<<trkAlg<<endl;
203        printf(" ERROR: CaloNuclei variables not completely filled \n");        printf(" ERROR: CaloNuclei variables not completely filled \n");
204        return;          return;  
205      };      };
# Line 214  void CaloNuclei::Process(Int_t ntr){ Line 220  void CaloNuclei::Process(Int_t ntr){
220      //      //
221      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
222      //      //
223        if ( !usepl18x && view==0 && plane==18 ) mip = 0.;
224        //
225      if ( ntr >= 0 ){      if ( ntr >= 0 ){
226        //        //
227        if ( strip != -1 &&        if ( strip != -1 &&
# Line 283  void CaloNuclei::Process(Int_t ntr){ Line 291  void CaloNuclei::Process(Int_t ntr){
291      //      //
292      mip = L2->GetCaloLevel1()->DecodeEstrip(ii,view,plane,strip);          mip = L2->GetCaloLevel1()->DecodeEstrip(ii,view,plane,strip);    
293      //      //
294        if ( !usepl18x && view==0 && plane==18 ) mip = 0.;
295        //
296        //
297      if ( ntr >= 0 ){      if ( ntr >= 0 ){
298        if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] &&        if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] &&
299             ( strip == (track->tibar[plane][view]-1) || strip == (track->tibar[plane][view]-2) || strip == (track->tibar[plane][view]) )             ( strip == (track->tibar[plane][view]-1) || strip == (track->tibar[plane][view]-2) || strip == (track->tibar[plane][view]) )
# Line 411  void CaloNuclei::Process(Int_t ntr){ Line 422  void CaloNuclei::Process(Int_t ntr){
422        //        //
423        mip = L2->GetCaloLevel1()->DecodeEstrip(ii,view,plane,strip);            mip = L2->GetCaloLevel1()->DecodeEstrip(ii,view,plane,strip);    
424        //        //
425          if ( !usepl18x && view==0 && plane==18 ) mip = 0.;
426          //
427          //
428        if ( strip != -1 ){        if ( strip != -1 ){
429          if ( view == 0 ){          if ( view == 0 ){
430            ipl = (1 + plane) * 2;            ipl = (1 + plane) * 2;
# Line 501  void CaloNuclei::Process(Int_t ntr){ Line 515  void CaloNuclei::Process(Int_t ntr){
515      Int_t uplim2 = interplane-1;      Int_t uplim2 = interplane-1;
516      //      //
517      while ( l < uplim && ind < interplane ){      while ( l < uplim && ind < interplane ){
518        qm = TMath::KOrdStat((Long_64_t)interplane,qme,(Long64_t)ind,work);        qm = TMath::KOrdStat((Long64_t)interplane,qme,(Long64_t)ind,work);
519        if ( qm >= qmt ){        if ( qm >= qmt ){
520          if ( l < 3 ){          if ( l < 3 ){
521            qpremean += qm;            qpremean += qm;
# Line 518  void CaloNuclei::Process(Int_t ntr){ Line 532  void CaloNuclei::Process(Int_t ntr){
532      l = 0;      l = 0;
533      RN = 0;      RN = 0;
534      while ( l < uplim && ind < interplane ){      while ( l < uplim && ind < interplane ){
535        qm2 = TMath::KOrdStat((Long_64_t)interplane,qme2,(Long_64_t)ind,work);        qm2 = TMath::KOrdStat((Long64_t)interplane,qme2,(Long64_t)ind,work);
536        if ( qm2 >= qmt ){                if ( qm2 >= qmt ){        
537          if ( l < N ){          if ( l < N ){
538            qpremeanN += qm2;            qpremeanN += qm2;
# Line 539  void CaloNuclei::Process(Int_t ntr){ Line 553  void CaloNuclei::Process(Int_t ntr){
553      RN = 0;      RN = 0;
554      S2=0;      S2=0;
555      while ( l < uplim2 && ind<interplane){      while ( l < uplim2 && ind<interplane){
556        qm2 = TMath::KOrdStat((Long_64_t)interplane,qme2,(Long_64_t)ind,work);        qm2 = TMath::KOrdStat((Long64_t)interplane,qme2,(Long64_t)ind,work);
557        if ( qm2 < qmt ) S2++;        if ( qm2 < qmt ) S2++;
558        ind++;        ind++;
559      }      }
# Line 548  void CaloNuclei::Process(Int_t ntr){ Line 562  void CaloNuclei::Process(Int_t ntr){
562      l = 0;      l = 0;
563      RN = 0;      RN = 0;
564      while ( l < uplim2 && ind < interplane ){      while ( l < uplim2 && ind < interplane ){
565        qm2 = TMath::KOrdStat((Long_64_t)interplane,qme2,(Long_64_t)ind,work);        qm2 = TMath::KOrdStat((Long64_t)interplane,qme2,(Long64_t)ind,work);
566        if ( qm2 >= qmt ){                if ( qm2 >= qmt ){        
567          if ( l < (interplane - 1 - S2)){            if ( l < (interplane - 1 - S2)){  
568            qNmin1_w += qm2;            qNmin1_w += qm2;
# Line 620  Float_t charge = 1000.; Line 634  Float_t charge = 1000.;
634  Float_t beta = 100.;  Float_t beta = 100.;
635    
636  //------- First try track dependent beta  //------- First try track dependent beta
637      if( L2->GetTrkLevel2()->GetNTracks()>=1 ){      if( L2->GetNTracks(trkAlg)>=1 ){
638      PamTrack *TRKtrack = L2->GetTrack(0);          PamTrack *TRKtrack = L2->GetTrack(0,trkAlg);
639      if (fabs(TRKtrack->GetToFTrack()->beta[12]) < 100.) beta = fabs(TRKtrack->GetToFTrack()->beta[12]);      if (fabs(TRKtrack->GetToFTrack()->beta[12]) < 100.) beta = fabs(TRKtrack->GetToFTrack()->beta[12]);
640                                                  }                                                  }
641  //------- If no beta found, try standalone beta  //------- If no beta found, try standalone beta
# Line 665  Float_t beta = 100.; Line 679  Float_t beta = 100.;
679  // //=======================================================================  // //=======================================================================
680  // //===========    charge determination Maximum release  vs. beta    ===============  // //===========    charge determination Maximum release  vs. beta    ===============
681  // //======================     Rome method    ===========================  // //======================     Rome method    ===========================
682  // //=======================================================================  // //=======================================================================    
       
     Float_t D0[7] = {0, 3 , 4 , 5 , 6, 8, 90};  
     Float_t E1[7] = {0 ,923.553 , 659.842, 1113.97, 3037.25, 3034.84, 0};  
     Float_t E2[7] = {0 ,6.92574 , 5.08865, 5.29349, 6.41442, 5.52969, 0};  
     Float_t E3[7] = {0 ,9.7227  , 13.18, 23.5444, 38.2057, 63.6784, 80000};  
683    
684      Float_t xx1[7],yy1[7];      Float_t D0[9] = {0, 1, 2, 3 , 4 , 5 , 6, 8, 90};
685      n1 = 7;      Float_t E1[9] = {0, 500, 500, 923.553 , 659.842, 1113.97, 3037.25, 3034.84, 0};
686        Float_t E2[9] = {0, 11.0, 7.5, 6.92574 , 5.08865, 5.29349, 6.41442, 5.52969, 0};
687        Float_t E3[9] = {0, 1.2, 4, 9.7227  , 13.18, 23.5444, 38.2057, 63.6784, 80000};
688    
689        Float_t xx1[9],yy1[9];
690        n1 = 9;
691            
692      charge = 1000.;      charge = 1000.;
693      mip=0;      mip=0;
# Line 714  Float_t beta = 100.; Line 728  Float_t beta = 100.;
728  // =======================================================================  // =======================================================================
729  // ===========    charge determination dedx  vs. beta    ===============  // ===========    charge determination dedx  vs. beta    ===============
730  // ======================     Rome method    ===========================  // ======================     Rome method    ===========================
731  // =======================================================================  // =======================================================================  
732        
733      Float_t F0[7] = {0.,3. ,4., 5. , 6., 8, 90};      Float_t F0[9] = {0., 1., 2., 3. ,4., 5. , 6., 8, 90};
734      Float_t G1[7] = {0 ,642.935 , 848.684, 1346.05, 3238.82, 3468.6, 0};      Float_t G1[9] = {0, 500, 500, 642.935 , 848.684, 1346.05, 3238.82, 3468.6, 0};
735      Float_t G2[7] = {0 ,6.2038 , 5.51723, 5.65265, 6.54089, 5.72723, 0};      Float_t G2[9] = {0, 11, 7.5, 6.2038 , 5.51723, 5.65265, 6.54089, 5.72723, 0};
736      Float_t G3[7] = {0 ,9.2421 , 13.9858, 25.3912, 39.6332, 64.5674, 80000};      Float_t G3[9] = {0, 1.2, 4, 9.2421 , 13.9858, 25.3912, 39.6332, 64.5674, 80000};
737    
738        
739      charge = 1000.;      charge = 1000.;
# Line 729  Float_t beta = 100.; Line 743  Float_t beta = 100.;
743      if (beta<2.) { // it makes no sense to allow beta=5 or so...      if (beta<2.) { // it makes no sense to allow beta=5 or so...
744                
745            
746        if( L2->GetTrkLevel2()->GetNTracks()>=1 ){        if( L2->GetNTracks(trkAlg)>=1 ){
747          mip=dedx1;          mip=dedx1;
748        }        }
749        if (mip==0) mip=stdedx1;        if (mip==0) mip=stdedx1;
# Line 763  Float_t beta = 100.; Line 777  Float_t beta = 100.;
777  //=======================================================================  //=======================================================================
778  //===========    charge determination dedx  vs. defl    ===============  //===========    charge determination dedx  vs. defl    ===============
779  //======================     Rome method    ===========================  //======================     Rome method    ===========================
780  //=======================================================================  //=======================================================================    
       
     //new  
      Float_t H0[7] = {0, 3 , 4 , 5 , 6, 8, 90 };  
      Float_t I1[7] = {0 , 56.1019, 101.673, 109.225, 150.599, 388.531, 0};  
      Float_t I2[7] = {0 , -12.5581, -22.5543, -15.9823, -28.2207, -93.6871, 0};  
      Float_t I3[7] = {0 , 11.6218, 19.664, 32.1817, 45.7527, 84.5992, 80000};  
       
   
 //     Float_t H0[7] = {0, 3 , 4 , 5 , 6, 8, 90 };  
 //     Float_t I1[7] = {0 , 56.1019, 101.673, 155, 150.599, 388.531, 0};  
 //     Float_t I2[7] = {0 , -12.5581, -22.5543, -35.6217, -28.2207, -93.6871, 0};  
 //     Float_t I3[7] = {0 , 11.6218, 19.664, 34.3311, 45.7527, 84.5992, 8000};  
781    
782         Float_t H0[9] = {0, 1, 2, 3 , 4 , 5 , 6, 8, 90 };
783         Float_t I1[9] = {0, 3.5, 40, 56.1019, 101.673, 109.225, 150.599, 388.531, 0};
784         Float_t I2[9] = {0, -1, -13.6, -12.5581, -22.5543, -15.9823, -28.2207, -93.6871, 0};
785         Float_t I3[9] = {0, 1, 5.3, 11.6218, 19.664, 32.1817, 45.7527, 84.5992, 80000};
786    
787            
788      charge = 1000.;      charge = 1000.;
# Line 786  Float_t beta = 100.; Line 792  Float_t beta = 100.;
792    
793      if (beta<2.) { // it makes no sense to allow beta=5 or so...      if (beta<2.) { // it makes no sense to allow beta=5 or so...
794    
795        if( L2->GetTrkLevel2()->GetNTracks()>=1 ){        if( L2->GetNTracks(trkAlg)>=1 ){
796          PamTrack *TRKtrack = L2->GetTrack(0);            PamTrack *TRKtrack = L2->GetTrack(0,trkAlg);
797          mip=dedx1;          mip=dedx1;
798          if (mip==0) mip=stdedx1;          if (mip==0) mip=stdedx1;
799          defl=TRKtrack->GetTrkTrack()->al[4];          defl=TRKtrack->GetTrkTrack()->al[4];
# Line 822  Float_t beta = 100.; Line 828  Float_t beta = 100.;
828  //============================================================================================  //============================================================================================
829  //===========    charge determination Truncated mean (N-1 planes) vs. defl ===================  //===========    charge determination Truncated mean (N-1 planes) vs. defl ===================
830  //================================     Rome method    ========================================  //================================     Rome method    ========================================
831  //============================================================================================  //============================================================================================  
832        
833      Float_t L0[7] = {0, 3 , 4 , 5 , 6, 8, 90};      Float_t L0[9] = {0, 1, 2, 3 , 4 , 5 , 6, 8, 90};
834      Float_t M1[7] = {0 , 63.0145, 120.504, 173.663, 245.33, 236.517, 0};      Float_t M1[9] = {0, 3.5, 27, 63.0145, 120.504, 173.663, 245.33, 236.517, 0};
835      Float_t M2[7] = {0 , -15.005, -31.0635, -39.4988, -60.5011, -46.3992, 0};      Float_t M2[9] = {0, -1, -10.6, -15.005, -31.0635, -39.4988, -60.5011, -46.3992, 0};
836      Float_t M3[7] = {0 , 12.5037, 22.8652, 35.2907, 51.4678, 86.4155, 8000};      Float_t M3[9] = {0, 1, 7, 12.5037, 22.8652, 35.2907, 51.4678, 86.4155, 80000};
837    
838      charge = 1000.;      charge = 1000.;
839      mip=0;      mip=0;
# Line 835  Float_t beta = 100.; Line 841  Float_t beta = 100.;
841    
842      if (beta<2.) { // it makes no sense to allow beta=5 or so...      if (beta<2.) { // it makes no sense to allow beta=5 or so...
843    
844        if( L2->GetTrkLevel2()->GetNTracks()>=1 ){        if( L2->GetNTracks(trkAlg)>=1 ){
845          mip=qNmin1;          mip=qNmin1;
846                    
847          if (mip>0 && defl<0.7 && defl>0)   {          if (mip>0 && defl<0.7 && defl>0)   {
Legend:
Removed from v.1.15  
changed lines
  Added in v.1.20
  ViewVC Help
Powered by ViewVC 1.1.23