/[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.11 by mocchiut, Mon Nov 26 08:48:24 2007 UTC revision 1.14 by malvezzi, Fri Nov 7 09:06:15 2008 UTC
# Line 1  Line 1 
1  #include <CaloNuclei.h>  #include <CaloNuclei.h>
2    #include <TGraph.h>
3    #include <TSpline.h>
4    #include <TMVA/TSpline1.h>
5    
6  //--------------------------------------  //--------------------------------------
7  /**  /**
# Line 23  CaloNuclei::CaloNuclei(PamLevel2 *l2p){ Line 26  CaloNuclei::CaloNuclei(PamLevel2 *l2p){
26    R = 3;    R = 3;
27    //    //
28    debug = false;    debug = false;
29      // debug = true;
30    usetrack = true;    usetrack = true;
31    //    //
32  };  };
# Line 41  void CaloNuclei::Clear(){ Line 45  void CaloNuclei::Clear(){
45    dedx3 = 0.;    dedx3 = 0.;
46    qpremean = 0.;    qpremean = 0.;
47    qpremeanN = 0.;    qpremeanN = 0.;
48      maxrel = 0;
49      qNmin1 = 0;
50      qNmin1_w = 0;
51      charge_siegen1 = 0;
52      ZCalo_maxrel_b = 0;
53      ZCalo_dedx_b = 0;
54      ZCalo_dedx_defl= 0;
55      ZCalo_Nmin1_defl= 0;
56    //    //
57    multhit = false;    multhit = false;
58    gap = false;    gap = false;
# Line 53  void CaloNuclei::Print(){ Line 65  void CaloNuclei::Print(){
65    //    //
66    printf("========================================================================\n");    printf("========================================================================\n");
67    printf(" OBT: %u PKT: %u ATIME: %u Track %i Use track %i \n",OBT,PKT,atime,tr,usetrack);    printf(" OBT: %u PKT: %u ATIME: %u Track %i Use track %i \n",OBT,PKT,atime,tr,usetrack);
68    printf(" interplane [number of available dE/dx before interaction]:.. %i\n",interplane);    printf(" interplane [number of available dE/dx before interaction]:....... %i\n",interplane);
69    printf(" ethr [threshold used to determine interplane]:.............. %f \n",ethr);    printf(" ethr [threshold used to determine interplane]:................... %f \n",ethr);
70    printf(" dedx1 [dE/dx from the first calorimeter plane]:............. %f \n",dedx1);    printf(" dedx1 [dE/dx from the first calorimeter plane]:.................. %f \n",dedx1);
71    printf(" stdedx1 [dE/dx from the first calorimeter plane standalone]: %f \n",stdedx1);    printf(" stdedx1 [dE/dx from the first calorimeter plane standalone]:..... %f \n",stdedx1);
72    printf(" dedx3 [dE/dx (average) if the first 3 Si planes]:........... %f \n",dedx3);    printf(" dedx3 [dE/dx (average) if the first 3 Si planes]:................ %f \n",dedx3);
73    printf(" multhit [true if interplane determined by multiple hits]:... %i \n",multhit);    printf(" multhit [true if interplane determined by multiple hits]:........ %i \n",multhit);
74    printf(" gap [true if interplane determined by a gap]:............... %i \n",gap);    printf(" gap [true if interplane determined by a gap]:.................... %i \n",gap);
75    printf(" preq [total energy in MIP before the interaction plane]:.... %f \n",preq);    printf(" preq [total energy in MIP before the interaction plane]:......... %f \n",preq);
76    printf(" postq [total energy in MIP after the interaction plane]:.... %f \n",postq);    printf(" postq [total energy in MIP after the interaction plane]:......... %f \n",postq);
77    printf(" qpremean [truncated mean using 3 planes and 3 strips]:...... %f \n",qpremean);    printf(" qpremean [truncated mean using 3 planes and 3 strips]:........... %f \n",qpremean);
78    printf(" N [no of used plane]:....................................... %i \n",N);    printf(" N [no of used plane]:............................................ %i \n",N);
79    printf(" R [no strip used per plane ]:............................... %i \n",R);    printf(" R [no strip used per plane ]:.................................... %i \n",R);
80    printf(" qpremeanN [truncated mean using N planes and R strips]:..... %f \n",qpremeanN);    printf(" qpremeanN [truncated mean using N planes and R strips]:.......... %f \n",qpremeanN);
81      printf(" qNmin1 [truncated mean using N-1 planes and R strips]: .......... %f \n",qNmin1);
82      printf(" maxrel [dE/dx of strip with maximum release (I plane)]:.......... %f \n",maxrel);
83      printf(" ZCalo_maxrel_b [Z from maximum release in I Calo plane vs beta].. %f \n",ZCalo_maxrel_b);
84      printf(" ZCalo_dedx_b [Z from dedx in I Calo plane vs beta].. ............ %f \n",ZCalo_dedx_b);
85      printf(" ZCalo_dedx_defl [Z from dedx in I Calo plane vs deflection....... %f \n",ZCalo_dedx_defl);
86      printf(" ZCalo_Nmin1_defl [Z from truncated mean (N-1 pl) vs deflection].. %f \n",ZCalo_Nmin1_defl);
87    printf("========================================================================\n");    printf("========================================================================\n");
88    //    //
89  };  };
# Line 112  void CaloNuclei::Process(Int_t ntr){ Line 130  void CaloNuclei::Process(Int_t ntr){
130    //    //
131    if ( debug ) printf(" Always calculate stdedx1 \n");    if ( debug ) printf(" Always calculate stdedx1 \n");
132    //    //
133    // Always calculate stdedx1    // Always calculate stdedx1 and maxrel
134    //    //
135      Int_t cont=0;
136    Int_t view = 0;    Int_t view = 0;
137    Int_t plane = 0;    Int_t plane = 0;
138    Int_t strip = 0;    Int_t strip = 0;
# Line 142  void CaloNuclei::Process(Int_t ntr){ Line 161  void CaloNuclei::Process(Int_t ntr){
161    // find energy released along the strip of maximum on the first plane and on the two neighbour strips    // find energy released along the strip of maximum on the first plane and on the two neighbour strips
162    //    //
163    if ( indx > 0 ){    if ( indx > 0 ){
164      Int_t mindx = (Int_t)TMath::LocMax(indx,stfpl);      Int_t mindx = (Int_t)TMath::LocMax(indx,vfpl);
165      for (Int_t ii=0; ii<indx; ii++){      for (Int_t ii=0; ii<indx; ii++){
166        if ( stfpl[ii] == stfpl[mindx] ) stdedx1 += vfpl[ii];        if ( stfpl[ii] == stfpl[mindx] ) stdedx1 += vfpl[ii];
167        if ( (mindx-1)>=0 && stfpl[ii] == stfpl[mindx-1] ) stdedx1 += vfpl[ii];        if ( (mindx-1)>=0 && stfpl[ii] == (stfpl[mindx]-1) ) stdedx1 += vfpl[ii];
168        if ( (mindx+1)<96 && stfpl[ii] == stfpl[mindx+1] ) stdedx1 += vfpl[ii];        if ( (mindx+1)<96 && stfpl[ii] == (stfpl[mindx]+1) ) stdedx1 += vfpl[ii];
169          //      if ( (mindx-1)>=0 && stfpl[ii] == stfpl[mindx-1] ) stdedx1 += vfpl[ii];
170          //      if ( (mindx+1)<96 && stfpl[ii] == stfpl[mindx+1] ) stdedx1 += vfpl[ii];
171      };      };
172      maxrel = vfpl[mindx];
173    } else {    } else {
174      stdedx1 = 0.;      stdedx1 = 0.;
175        maxrel = 0.;
176    };    };
177      // cout<<stdedx1<<"      "<<maxrel<<"\n";
178    //    //
179    if ( debug ) printf(" if ( !usetrack ) return: usetrack %i ntr %i \n",usetrack,ntr);    if ( debug ) printf(" if ( !usetrack ) return: usetrack %i ntr %i \n",usetrack,ntr);
180    //    //
# Line 370  void CaloNuclei::Process(Int_t ntr){ Line 394  void CaloNuclei::Process(Int_t ntr){
394      //      //
395      // Calculate preq, postq, qpremean      // Calculate preq, postq, qpremean
396      //      //
397        cont++;
398      ii = 0;      ii = 0;
399      Int_t ind = -1;      Int_t ind = -1;
400      Int_t qsplane = -1;      Int_t qsplane = -1;
# Line 452  void CaloNuclei::Process(Int_t ntr){ Line 477  void CaloNuclei::Process(Int_t ntr){
477        ii++;        ii++;
478        //            //    
479      };      };
480      //  
481    
482     //
483      // here we must calculate qpremean, order vector qme, select 3 lowest measurements and caculate the mean...      // here we must calculate qpremean, order vector qme, select 3 lowest measurements and caculate the mean...
484      //          //    
485      if ( debug ){      if ( debug ){
# Line 471  void CaloNuclei::Process(Int_t ntr){ Line 498  void CaloNuclei::Process(Int_t ntr){
498      Float_t qmt = ethr*0.8; // *0.9      Float_t qmt = ethr*0.8; // *0.9
499      //      //
500      Int_t uplim = TMath::Max(3,N);      Int_t uplim = TMath::Max(3,N);
501        Int_t uplim2 = interplane-1;
502      //      //
503      while ( l < uplim && ind < interplane ){      while ( l < uplim && ind < interplane ){
504        qm = TMath::KOrdStat(interplane,qme,ind,work);        qm = TMath::KOrdStat(interplane,qme,ind,work);
# Line 486  void CaloNuclei::Process(Int_t ntr){ Line 514  void CaloNuclei::Process(Int_t ntr){
514      };      };
515      //      //
516      qpremean /= (Float_t)RN;      qpremean /= (Float_t)RN;
     //  
517      ind = 0;      ind = 0;
518      l = 0;      l = 0;
519      RN = 0;      RN = 0;
# Line 502  void CaloNuclei::Process(Int_t ntr){ Line 529  void CaloNuclei::Process(Int_t ntr){
529        };        };
530        ind++;        ind++;
531      };      };
532      //      ////////////////////////////////////
533        //to calculate qNmin1///////////////
534        ///////////////////////////////////
535        //values under threshold
536        qm2=0;
537        ind = 0;
538        l = 0;
539        RN = 0;
540        S2=0;
541        while ( l < uplim2 && ind<interplane){
542          qm2 = TMath::KOrdStat(interplane,qme2,ind,work);
543          if ( qm2 < qmt ) S2++;
544          ind++;
545        }
546        qm2=0;
547        ind = 0;
548        l = 0;
549        RN = 0;
550        while ( l < uplim2 && ind < interplane ){
551          qm2 = TMath::KOrdStat(interplane,qme2,ind,work);
552          if ( qm2 >= qmt ){        
553            if ( l < (interplane - 1 - S2)){  
554              qNmin1_w += qm2;
555              RN++;
556            };
557            l++;
558            if ( debug ) printf(" qm2 value no %i qm %f qmt %f RN %i \n",l,qm2,qmt,RN);
559          };
560          ind++;
561        };
562      qpremeanN /= (Float_t)RN;      qpremeanN /= (Float_t)RN;
563        qNmin1_w /= (Float_t)RN;
564      UN = RN;      UN = RN;
565        ///////set qNmin1 definition///////////
566        if (interplane==1 || interplane==2){  
567          if (dedx1>0) qNmin1=dedx1;
568          else if (stdedx1>0) qNmin1=stdedx1;              
569        }
570        else if (interplane > 2){
571          qNmin1 =  qNmin1_w;
572        }
573        ////////////////////////////////////
574        //////////////////////////////////
575      //      //
576      if ( debug ) printf(" charge is %f \n",sqrt(qpremean));      if ( debug ) printf(" charge is %f \n",sqrt(qpremean));
577      //      //
# Line 522  void CaloNuclei::Process(Int_t ntr){ Line 589  void CaloNuclei::Process(Int_t ntr){
589          postq = 0.;          postq = 0.;
590          qpremean = 0.;          qpremean = 0.;
591          qpremeanN = 0.;          qpremeanN = 0.;
592            qNmin1 = 0;
593          multhit = false;          multhit = false;
594          gap = false;          gap = false;
595          goto retry;          goto retry;
596        };        };
597      };      };
598    };    };
599    
600    
601    
602    //=======================================================================
603    //===========    charge determination stdedx1 vs. beta    ===============
604    //======================     Siegen method    ===========================
605    //=======================================================================
606    
607    // Data from file Calo_Bands_New_7.dat
608    Float_t C0[9] = {0 , 1 , 2 , 3 , 4 , 5 , 6 , 8 , 90 };
609    Float_t B0[9] = {0 , -2.03769 , 7.61781 , 19.7098 , 60.5598 , 57.9226 , 14.8368 , -1358.83 , 8200 };
610    Float_t B1[9] = {0 , 0.0211274 , 9.32057e-010 , 4.47241e-07 , 1.44826e-06 , 2.6189e-05 , 0.00278178 , 55.5445 , 0 };
611    Float_t B2[9] = {0 , -3.91742 , -20.0359 , -16.3043 , -16.9471 , -14.4622 , -10.9594 , -2.38014 , 0 };
612    Float_t B3[9] = {0 , 11.1469 , -6.63105 , -27.8834 , -132.044 , -55.341 , 173.25 , 4115 , 0 };
613    Float_t B4[9] = {0 , -14.3465 , -0.485215 , 18.8122 , 117.533 , -14.0898 , -325.269 , -4388.89 , 0 };
614    Float_t B5[9] = {0 , 6.24281 , 3.96018 , 0 , -26.1881 , 42.9731 , 182.697 , 1661.01 , 0 };
615    
616    Float_t x1[9],y1[9];
617    Int_t n1 = 9;
618    
619    Float_t charge = 1000.;
620    Float_t beta = 100.;
621    
622    //------- First try track dependent beta
623        if( L2->GetTrkLevel2()->GetNTracks()>=1 ){
624        PamTrack *TRKtrack = L2->GetTrack(0);
625        if (fabs(TRKtrack->GetToFTrack()->beta[12]) < 100.) beta = fabs(TRKtrack->GetToFTrack()->beta[12]);
626                                                    }
627    //------- If no beta found, try standalone beta
628        if (beta == 100.) {
629        ToFTrkVar *ttrack = L2->GetToFStoredTrack(-1);
630        beta = fabs(ttrack->beta[12]);
631                          }
632    
633        if (beta<2.) {  // it makes no sense to allow beta=5 or so...
634    
635        Float_t  mip=0;
636        mip=stdedx1 ;
637    
638        if (mip>0)   {
639    
640        Float_t betahelp =   pow(beta, 1.8);
641        Float_t ym = mip*betahelp;
642        Float_t xb = beta;
643    
644        for ( Int_t jj=0; jj<9; jj++ ){
645        x1[jj] = B0[jj]+B1[jj]*pow(xb,B2[jj])+B3[jj]*xb+B4[jj]*xb*xb+B5[jj]*xb*xb*xb;
646        y1[jj] = C0[jj]*C0[jj] ;
647                                      }
648    
649        TGraph *gr1 = new TGraph(n1,x1,y1);
650        TSpline3 *spl1 = new TSpline3("grs",gr1);    // use a cubic spline
651        Float_t chelp = spl1->Eval(ym);
652        charge = TMath::Sqrt(chelp);
653        gr1->Delete();
654        spl1->Delete();
655    
656                        } // if (mip1>0)
657                        } // beta < 100
658    
659    
660        charge_siegen1 = charge;
661    
662    //=======================  end charge Siegen  ===========================
663    
664    
665    // //=======================================================================
666    // //===========    charge determination Maximum release  vs. beta    ===============
667    // //======================     Rome method    ===========================
668    // //=======================================================================
669        
670        Float_t D0[7] = {0, 3 , 4 , 5 , 6, 8, 90};
671        Float_t E1[7] = {0 ,923.553 , 659.842, 1113.97, 3037.25, 3034.84, 0};
672        Float_t E2[7] = {0 ,6.92574 , 5.08865, 5.29349, 6.41442, 5.52969, 0};
673        Float_t E3[7] = {0 ,9.7227  , 13.18, 23.5444, 38.2057, 63.6784, 80000};
674    
675        Float_t xx1[7],yy1[7];
676        n1 = 7;
677        
678        charge = 1000.;
679        mip=0;
680        
681    
682        if (beta<2.) {  // it makes no sense to allow beta=5 or so...
683          
684          
685          mip=maxrel;
686    
687          if (mip>0)   {
688            Float_t ym = mip;
689            Float_t xb = beta;
690            
691            for ( Int_t jj=0; jj<n1; jj++ ){
692              xx1[jj] = E1[jj]*exp(-E2[jj]*xb)+E3[jj];
693              yy1[jj] = D0[jj]*D0[jj] ;
694            }
695            
696            TGraph *gr1 = new TGraph(n1,xx1,yy1);
697            TSpline3 *spl1 = new TSpline3("grs",gr1);    // use a cubic spline
698            Float_t chelp = spl1->Eval(ym);
699            charge = TMath::Sqrt(chelp);
700            gr1->Delete();
701            spl1->Delete();
702    
703            
704          } // if (mip1>0)
705        } // beta < 100
706    
707    
708        ZCalo_maxrel_b = charge;
709    
710        //=======================  end charge Rome: maxril vs beta  ===========================
711    
712    
713    
714    // =======================================================================
715    // ===========    charge determination dedx  vs. beta    ===============
716    // ======================     Rome method    ===========================
717    // =======================================================================
718        
719        Float_t F0[7] = {0.,3. ,4., 5. , 6., 8, 90};
720        Float_t G1[7] = {0 ,642.935 , 848.684, 1346.05, 3238.82, 3468.6, 0};
721        Float_t G2[7] = {0 ,6.2038 , 5.51723, 5.65265, 6.54089, 5.72723, 0};
722        Float_t G3[7] = {0 ,9.2421 , 13.9858, 25.3912, 39.6332, 64.5674, 80000};
723    
724      
725        charge = 1000.;
726        mip=0;
727    
728    
729        if (beta<2.) { // it makes no sense to allow beta=5 or so...
730          
731        
732          if( L2->GetTrkLevel2()->GetNTracks()>=1 ){
733            mip=dedx1;
734          }
735          if (mip==0) mip=stdedx1;
736          
737    
738          if (mip>0)   {
739            
740            Float_t ym = mip;
741            Float_t xb = beta;
742            
743            for ( Int_t jj=0; jj<n1; jj++ ){
744             xx1[jj] = G1[jj]*exp(-G2[jj]*xb)+G3[jj];
745             yy1[jj] = F0[jj]*F0[jj] ;
746            }
747            
748            TGraph *gr1 = new TGraph(n1,xx1,yy1);
749            TSpline3 *spl1 = new TSpline3("grs",gr1);    // use a cubic spline
750            Float_t chelp = spl1->Eval(ym);
751            charge = TMath::Sqrt(chelp);
752            gr1->Delete();
753            spl1->Delete();
754            
755          } //if (mip1>0)
756        } //beta < 100
757        
758        ZCalo_dedx_b = charge;
759    
760        //=======================  end charge Rome: dedx vs beta  ===========================
761        
762        
763    //=======================================================================
764    //===========    charge determination dedx  vs. defl    ===============
765    //======================     Rome method    ===========================
766    //=======================================================================
767        
768        //new
769         Float_t H0[7] = {0, 3 , 4 , 5 , 6, 8, 90 };
770         Float_t I1[7] = {0 , 56.1019, 101.673, 109.225, 150.599, 388.531, 0};
771         Float_t I2[7] = {0 , -12.5581, -22.5543, -15.9823, -28.2207, -93.6871, 0};
772         Float_t I3[7] = {0 , 11.6218, 19.664, 32.1817, 45.7527, 84.5992, 80000};
773        
774    
775    //     Float_t H0[7] = {0, 3 , 4 , 5 , 6, 8, 90 };
776    //     Float_t I1[7] = {0 , 56.1019, 101.673, 155, 150.599, 388.531, 0};
777    //     Float_t I2[7] = {0 , -12.5581, -22.5543, -35.6217, -28.2207, -93.6871, 0};
778    //     Float_t I3[7] = {0 , 11.6218, 19.664, 34.3311, 45.7527, 84.5992, 8000};
779    
780    
781        
782        charge = 1000.;
783        mip=0;
784        Float_t defl=0;
785        
786    
787        if (beta<2.) { // it makes no sense to allow beta=5 or so...
788    
789          if( L2->GetTrkLevel2()->GetNTracks()>=1 ){
790            PamTrack *TRKtrack = L2->GetTrack(0);
791            mip=dedx1;
792            if (mip==0) mip=stdedx1;
793            defl=TRKtrack->GetTrkTrack()->al[4];
794            
795            
796            if (mip>0 && defl<0.7 && defl>0)   {
797              
798              Float_t ym = mip;
799              Float_t xb = defl;
800              
801              for ( Int_t jj=0; jj<n1; jj++ ){
802                xx1[jj] = I1[jj]*xb*xb+I2[jj]*xb+I3[jj];
803                yy1[jj] = H0[jj]*H0[jj] ;
804              }
805              
806              TGraph *gr1 = new TGraph(n1,xx1,yy1);
807              TSpline3 *spl1 = new TSpline3("grs",gr1);    // use a cubic spline
808              Float_t chelp = spl1->Eval(ym);
809              charge = TMath::Sqrt(chelp);
810              gr1->Delete();
811              spl1->Delete();
812              
813            } // if (mip1>0 && defl<0.5 && defl>0)
814          }//Ntrack>=1
815        } //beta < 100
816    
817        ZCalo_dedx_defl = charge;
818    
819    //=======================  end charge Rome: dedx vs defl  ===========================
820    
821    
822    //============================================================================================
823    //===========    charge determination Truncated mean (N-1 planes) vs. defl ===================
824    //================================     Rome method    ========================================
825    //============================================================================================
826        
827        Float_t L0[7] = {0, 3 , 4 , 5 , 6, 8, 90};
828        Float_t M1[7] = {0 , 63.0145, 120.504, 173.663, 245.33, 236.517, 0};
829        Float_t M2[7] = {0 , -15.005, -31.0635, -39.4988, -60.5011, -46.3992, 0};
830        Float_t M3[7] = {0 , 12.5037, 22.8652, 35.2907, 51.4678, 86.4155, 8000};
831    
832        charge = 1000.;
833        mip=0;
834      
835    
836        if (beta<2.) { // it makes no sense to allow beta=5 or so...
837    
838          if( L2->GetTrkLevel2()->GetNTracks()>=1 ){
839            mip=qNmin1;
840            
841            if (mip>0 && defl<0.7 && defl>0)   {
842              
843              Float_t ym = mip;
844              Float_t xb = defl;
845              
846              for ( Int_t jj=0; jj<n1; jj++ ){
847                xx1[jj] = M1[jj]*xb*xb+M2[jj]*xb+M3[jj];
848                yy1[jj] = L0[jj]*L0[jj] ;
849              }
850              
851              TGraph *gr1 = new TGraph(n1,xx1,yy1);
852              TSpline3 *spl1 = new TSpline3("grs",gr1);    // use a cubic spline
853              Float_t chelp = spl1->Eval(ym);
854              charge = TMath::Sqrt(chelp);
855              gr1->Delete();
856              spl1->Delete();
857              
858            } // if (mip1>0 && defl<0.5 && defl>0)
859          }//Ntrack>=1
860        } //beta < 100
861        
862        ZCalo_Nmin1_defl = charge;
863    
864    //=======================  end charge Rome: Nmin1 vs defl  ===========================
865    
866    
867    
868    
869    
870    //    //
871    if ( debug ) this->Print();    if ( debug ) this->Print();
872    if ( debug ) printf(" esci \n");    if ( debug ) printf(" esci \n");
Legend:
Removed from v.1.11  
changed lines
  Added in v.1.14
  ViewVC Help
Powered by ViewVC 1.1.23