/[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.20 by pam-fi, Mon Apr 20 09:12:05 2015 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  /**  /**
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 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      usepl18x = false;
32      trkAlg = alg;
33    //    //
34  };  };
35    
# Line 41  void CaloNuclei::Clear(){ Line 47  void CaloNuclei::Clear(){
47    dedx3 = 0.;    dedx3 = 0.;
48    qpremean = 0.;    qpremean = 0.;
49    qpremeanN = 0.;    qpremeanN = 0.;
50      maxrel = 0;
51      qNmin1 = 0;
52      qNmin1_w = 0;
53      charge_siegen1 = 0;
54      ZCalo_maxrel_b = 0;
55      ZCalo_dedx_b = 0;
56      ZCalo_dedx_defl= 0;
57      ZCalo_Nmin1_defl= 0;
58    //    //
59    multhit = false;    multhit = false;
60    gap = false;    gap = false;
# Line 53  void CaloNuclei::Print(){ Line 67  void CaloNuclei::Print(){
67    //    //
68    printf("========================================================================\n");    printf("========================================================================\n");
69    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);
70    printf(" interplane [number of available dE/dx before interaction]:.. %i\n",interplane);    printf(" interplane [number of available dE/dx before interaction]:....... %i\n",interplane);
71    printf(" ethr [threshold used to determine interplane]:.............. %f \n",ethr);    printf(" ethr [threshold used to determine interplane]:................... %f \n",ethr);
72    printf(" dedx1 [dE/dx from the first calorimeter plane]:............. %f \n",dedx1);    printf(" dedx1 [dE/dx from the first calorimeter plane]:.................. %f \n",dedx1);
73    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);
74    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);
75    printf(" multhit [true if interplane determined by multiple hits]:... %i \n",multhit);    printf(" multhit [true if interplane determined by multiple hits]:........ %i \n",multhit);
76    printf(" gap [true if interplane determined by a gap]:............... %i \n",gap);    printf(" gap [true if interplane determined by a gap]:.................... %i \n",gap);
77    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);
78    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);
79    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);
80    printf(" N [no of used plane]:....................................... %i \n",N);    printf(" N [no of used plane]:............................................ %i \n",N);
81    printf(" R [no strip used per plane ]:............................... %i \n",R);    printf(" R [no strip used per plane ]:.................................... %i \n",R);
82    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);
83      printf(" qNmin1 [truncated mean using N-1 planes and R strips]: .......... %f \n",qNmin1);
84      printf(" maxrel [dE/dx of strip with maximum release (I plane)]:.......... %f \n",maxrel);
85      printf(" ZCalo_maxrel_b [Z from maximum release in I Calo plane vs beta].. %f \n",ZCalo_maxrel_b);
86      printf(" ZCalo_dedx_b [Z from dedx in I Calo plane vs beta].. ............ %f \n",ZCalo_dedx_b);
87      printf(" ZCalo_dedx_defl [Z from dedx in I Calo plane vs deflection....... %f \n",ZCalo_dedx_defl);
88      printf(" ZCalo_Nmin1_defl [Z from truncated mean (N-1 pl) vs deflection].. %f \n",ZCalo_Nmin1_defl);
89    printf("========================================================================\n");    printf("========================================================================\n");
90    //    //
91  };  };
# Line 112  void CaloNuclei::Process(Int_t ntr){ Line 132  void CaloNuclei::Process(Int_t ntr){
132    //    //
133    if ( debug ) printf(" Always calculate stdedx1 \n");    if ( debug ) printf(" Always calculate stdedx1 \n");
134    //    //
135    // Always calculate stdedx1    // Always calculate stdedx1 and maxrel
136    //    //
137      Int_t cont=0;
138    Int_t view = 0;    Int_t view = 0;
139    Int_t plane = 0;    Int_t plane = 0;
140    Int_t strip = 0;    Int_t strip = 0;
# Line 127  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 142  void CaloNuclei::Process(Int_t ntr){ Line 166  void CaloNuclei::Process(Int_t ntr){
166    // 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
167    //    //
168    if ( indx > 0 ){    if ( indx > 0 ){
169      Int_t mindx = (Int_t)TMath::LocMax(indx,stfpl);      Int_t mindx = (Int_t)TMath::LocMax(indx,vfpl);
170      for (Int_t ii=0; ii<indx; ii++){      for (Int_t ii=0; ii<indx; ii++){
171        if ( stfpl[ii] == stfpl[mindx] ) stdedx1 += vfpl[ii];        if ( stfpl[ii] == stfpl[mindx] ) stdedx1 += vfpl[ii];
172        if ( (mindx-1)>=0 && stfpl[ii] == stfpl[mindx-1] ) stdedx1 += vfpl[ii];        if ( (mindx-1)>=0 && stfpl[ii] == (stfpl[mindx]-1) ) stdedx1 += vfpl[ii];
173        if ( (mindx+1)<96 && stfpl[ii] == stfpl[mindx+1] ) stdedx1 += vfpl[ii];        if ( (mindx+1)<96 && stfpl[ii] == (stfpl[mindx]+1) ) stdedx1 += vfpl[ii];
174          //      if ( (mindx-1)>=0 && stfpl[ii] == stfpl[mindx-1] ) stdedx1 += vfpl[ii];
175          //      if ( (mindx+1)<96 && stfpl[ii] == stfpl[mindx+1] ) stdedx1 += vfpl[ii];
176      };      };
177      maxrel = vfpl[mindx];
178    } else {    } else {
179      stdedx1 = 0.;      stdedx1 = 0.;
180        maxrel = 0.;
181    };    };
182      // cout<<stdedx1<<"      "<<maxrel<<"\n";
183    //    //
184    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);
185    //    //
# Line 162  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 190  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 259  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 370  void CaloNuclei::Process(Int_t ntr){ Line 405  void CaloNuclei::Process(Int_t ntr){
405      //      //
406      // Calculate preq, postq, qpremean      // Calculate preq, postq, qpremean
407      //      //
408        cont++;
409      ii = 0;      ii = 0;
410      Int_t ind = -1;      Int_t ind = -1;
411      Int_t qsplane = -1;      Int_t qsplane = -1;
# Line 386  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 452  void CaloNuclei::Process(Int_t ntr){ Line 491  void CaloNuclei::Process(Int_t ntr){
491        ii++;        ii++;
492        //            //    
493      };      };
494      //  
495    
496     //
497      // 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...
498      //          //    
499      if ( debug ){      if ( debug ){
# Line 471  void CaloNuclei::Process(Int_t ntr){ Line 512  void CaloNuclei::Process(Int_t ntr){
512      Float_t qmt = ethr*0.8; // *0.9      Float_t qmt = ethr*0.8; // *0.9
513      //      //
514      Int_t uplim = TMath::Max(3,N);      Int_t uplim = TMath::Max(3,N);
515        Int_t uplim2 = interplane-1;
516      //      //
517      while ( l < uplim && ind < interplane ){      while ( l < uplim && ind < interplane ){
518        qm = TMath::KOrdStat(interplane,qme,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 486  void CaloNuclei::Process(Int_t ntr){ Line 528  void CaloNuclei::Process(Int_t ntr){
528      };      };
529      //      //
530      qpremean /= (Float_t)RN;      qpremean /= (Float_t)RN;
     //  
531      ind = 0;      ind = 0;
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(interplane,qme2,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 502  void CaloNuclei::Process(Int_t ntr){ Line 543  void CaloNuclei::Process(Int_t ntr){
543        };        };
544        ind++;        ind++;
545      };      };
546      //      ////////////////////////////////////
547        //to calculate qNmin1///////////////
548        ///////////////////////////////////
549        //values under threshold
550        qm2=0;
551        ind = 0;
552        l = 0;
553        RN = 0;
554        S2=0;
555        while ( l < uplim2 && ind<interplane){
556          qm2 = TMath::KOrdStat((Long64_t)interplane,qme2,(Long64_t)ind,work);
557          if ( qm2 < qmt ) S2++;
558          ind++;
559        }
560        qm2=0;
561        ind = 0;
562        l = 0;
563        RN = 0;
564        while ( l < uplim2 && ind < interplane ){
565          qm2 = TMath::KOrdStat((Long64_t)interplane,qme2,(Long64_t)ind,work);
566          if ( qm2 >= qmt ){        
567            if ( l < (interplane - 1 - S2)){  
568              qNmin1_w += qm2;
569              RN++;
570            };
571            l++;
572            if ( debug ) printf(" qm2 value no %i qm %f qmt %f RN %i \n",l,qm2,qmt,RN);
573          };
574          ind++;
575        };
576      qpremeanN /= (Float_t)RN;      qpremeanN /= (Float_t)RN;
577        qNmin1_w /= (Float_t)RN;
578      UN = RN;      UN = RN;
579        ///////set qNmin1 definition///////////
580        if (interplane==1 || interplane==2){  
581          if (dedx1>0) qNmin1=dedx1;
582          else if (stdedx1>0) qNmin1=stdedx1;              
583        }
584        else if (interplane > 2){
585          qNmin1 =  qNmin1_w;
586        }
587        ////////////////////////////////////
588        //////////////////////////////////
589      //      //
590      if ( debug ) printf(" charge is %f \n",sqrt(qpremean));      if ( debug ) printf(" charge is %f \n",sqrt(qpremean));
591      //      //
# Line 522  void CaloNuclei::Process(Int_t ntr){ Line 603  void CaloNuclei::Process(Int_t ntr){
603          postq = 0.;          postq = 0.;
604          qpremean = 0.;          qpremean = 0.;
605          qpremeanN = 0.;          qpremeanN = 0.;
606            qNmin1 = 0;
607          multhit = false;          multhit = false;
608          gap = false;          gap = false;
609          goto retry;          goto retry;
610        };        };
611      };      };
612    };    };
613    
614    
615    
616    //=======================================================================
617    //===========    charge determination stdedx1 vs. beta    ===============
618    //======================     Siegen method    ===========================
619    //=======================================================================
620    
621    // Data from file Calo_Bands_New_7.dat
622    Float_t C0[9] = {0 , 1 , 2 , 3 , 4 , 5 , 6 , 8 , 90 };
623    Float_t B0[9] = {0 , -2.03769 , 7.61781 , 19.7098 , 60.5598 , 57.9226 , 14.8368 , -1358.83 , 8200 };
624    Float_t B1[9] = {0 , 0.0211274 , 9.32057e-010 , 4.47241e-07 , 1.44826e-06 , 2.6189e-05 , 0.00278178 , 55.5445 , 0 };
625    Float_t B2[9] = {0 , -3.91742 , -20.0359 , -16.3043 , -16.9471 , -14.4622 , -10.9594 , -2.38014 , 0 };
626    Float_t B3[9] = {0 , 11.1469 , -6.63105 , -27.8834 , -132.044 , -55.341 , 173.25 , 4115 , 0 };
627    Float_t B4[9] = {0 , -14.3465 , -0.485215 , 18.8122 , 117.533 , -14.0898 , -325.269 , -4388.89 , 0 };
628    Float_t B5[9] = {0 , 6.24281 , 3.96018 , 0 , -26.1881 , 42.9731 , 182.697 , 1661.01 , 0 };
629    
630    Float_t x1[9],y1[9];
631    Int_t n1 = 9;
632    
633    Float_t charge = 1000.;
634    Float_t beta = 100.;
635    
636    //------- First try track dependent beta
637        if( L2->GetNTracks(trkAlg)>=1 ){
638            PamTrack *TRKtrack = L2->GetTrack(0,trkAlg);
639        if (fabs(TRKtrack->GetToFTrack()->beta[12]) < 100.) beta = fabs(TRKtrack->GetToFTrack()->beta[12]);
640                                                    }
641    //------- If no beta found, try standalone beta
642        if (beta == 100.) {
643        ToFTrkVar *ttrack = L2->GetToFStoredTrack(-1);
644        beta = fabs(ttrack->beta[12]);
645                          }
646    
647        if (beta<2.) {  // it makes no sense to allow beta=5 or so...
648    
649        Float_t  mip=0;
650        mip=stdedx1 ;
651    
652        if (mip>0)   {
653    
654        Float_t betahelp =   pow(beta, 1.8);
655        Float_t ym = mip*betahelp;
656        Float_t xb = beta;
657    
658        for ( Int_t jj=0; jj<9; jj++ ){
659        x1[jj] = B0[jj]+B1[jj]*pow(xb,B2[jj])+B3[jj]*xb+B4[jj]*xb*xb+B5[jj]*xb*xb*xb;
660        y1[jj] = C0[jj]*C0[jj] ;
661                                      }
662    
663        TGraph *gr1 = new TGraph(n1,x1,y1);
664        TSpline3 *spl1 = new TSpline3("grs",gr1);    // use a cubic spline
665        Float_t chelp = spl1->Eval(ym);
666        charge = TMath::Sqrt(chelp);
667        gr1->Delete();
668        spl1->Delete();
669    
670                        } // if (mip1>0)
671                        } // beta < 100
672    
673    
674        charge_siegen1 = charge;
675    
676    //=======================  end charge Siegen  ===========================
677    
678    
679    // //=======================================================================
680    // //===========    charge determination Maximum release  vs. beta    ===============
681    // //======================     Rome method    ===========================
682    // //=======================================================================    
683    
684        Float_t D0[9] = {0, 1, 2, 3 , 4 , 5 , 6, 8, 90};
685        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.;
693        mip=0;
694        
695    
696        if (beta<2.) {  // it makes no sense to allow beta=5 or so...
697          
698          
699          mip=maxrel;
700    
701          if (mip>0)   {
702            Float_t ym = mip;
703            Float_t xb = beta;
704            
705            for ( Int_t jj=0; jj<n1; jj++ ){
706              xx1[jj] = E1[jj]*exp(-E2[jj]*xb)+E3[jj];
707              yy1[jj] = D0[jj]*D0[jj] ;
708            }
709            
710            TGraph *gr1 = new TGraph(n1,xx1,yy1);
711            TSpline3 *spl1 = new TSpline3("grs",gr1);    // use a cubic spline
712            Float_t chelp = spl1->Eval(ym);
713            charge = TMath::Sqrt(chelp);
714            gr1->Delete();
715            spl1->Delete();
716    
717            
718          } // if (mip1>0)
719        } // beta < 100
720    
721    
722        ZCalo_maxrel_b = charge;
723    
724        //=======================  end charge Rome: maxril vs beta  ===========================
725    
726    
727    
728    // =======================================================================
729    // ===========    charge determination dedx  vs. beta    ===============
730    // ======================     Rome method    ===========================
731    // =======================================================================  
732    
733        Float_t F0[9] = {0., 1., 2., 3. ,4., 5. , 6., 8, 90};
734        Float_t G1[9] = {0, 500, 500, 642.935 , 848.684, 1346.05, 3238.82, 3468.6, 0};
735        Float_t G2[9] = {0, 11, 7.5, 6.2038 , 5.51723, 5.65265, 6.54089, 5.72723, 0};
736        Float_t G3[9] = {0, 1.2, 4, 9.2421 , 13.9858, 25.3912, 39.6332, 64.5674, 80000};
737    
738      
739        charge = 1000.;
740        mip=0;
741    
742    
743        if (beta<2.) { // it makes no sense to allow beta=5 or so...
744          
745        
746          if( L2->GetNTracks(trkAlg)>=1 ){
747            mip=dedx1;
748          }
749          if (mip==0) mip=stdedx1;
750          
751    
752          if (mip>0)   {
753            
754            Float_t ym = mip;
755            Float_t xb = beta;
756            
757            for ( Int_t jj=0; jj<n1; jj++ ){
758             xx1[jj] = G1[jj]*exp(-G2[jj]*xb)+G3[jj];
759             yy1[jj] = F0[jj]*F0[jj] ;
760            }
761            
762            TGraph *gr1 = new TGraph(n1,xx1,yy1);
763            TSpline3 *spl1 = new TSpline3("grs",gr1);    // use a cubic spline
764            Float_t chelp = spl1->Eval(ym);
765            charge = TMath::Sqrt(chelp);
766            gr1->Delete();
767            spl1->Delete();
768            
769          } //if (mip1>0)
770        } //beta < 100
771        
772        ZCalo_dedx_b = charge;
773    
774        //=======================  end charge Rome: dedx vs beta  ===========================
775        
776        
777    //=======================================================================
778    //===========    charge determination dedx  vs. defl    ===============
779    //======================     Rome method    ===========================
780    //=======================================================================    
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.;
789        mip=0;
790        Float_t defl=0;
791        
792    
793        if (beta<2.) { // it makes no sense to allow beta=5 or so...
794    
795          if( L2->GetNTracks(trkAlg)>=1 ){
796              PamTrack *TRKtrack = L2->GetTrack(0,trkAlg);
797            mip=dedx1;
798            if (mip==0) mip=stdedx1;
799            defl=TRKtrack->GetTrkTrack()->al[4];
800            
801            
802            if (mip>0 && defl<0.7 && defl>0)   {
803              
804              Float_t ym = mip;
805              Float_t xb = defl;
806              
807              for ( Int_t jj=0; jj<n1; jj++ ){
808                xx1[jj] = I1[jj]*xb*xb+I2[jj]*xb+I3[jj];
809                yy1[jj] = H0[jj]*H0[jj] ;
810              }
811              
812              TGraph *gr1 = new TGraph(n1,xx1,yy1);
813              TSpline3 *spl1 = new TSpline3("grs",gr1);    // use a cubic spline
814              Float_t chelp = spl1->Eval(ym);
815              charge = TMath::Sqrt(chelp);
816              gr1->Delete();
817              spl1->Delete();
818              
819            } // if (mip1>0 && defl<0.5 && defl>0)
820          }//Ntrack>=1
821        } //beta < 100
822    
823        ZCalo_dedx_defl = charge;
824    
825    //=======================  end charge Rome: dedx vs defl  ===========================
826    
827    
828    //============================================================================================
829    //===========    charge determination Truncated mean (N-1 planes) vs. defl ===================
830    //================================     Rome method    ========================================
831    //============================================================================================  
832    
833        Float_t L0[9] = {0, 1, 2, 3 , 4 , 5 , 6, 8, 90};
834        Float_t M1[9] = {0, 3.5, 27, 63.0145, 120.504, 173.663, 245.33, 236.517, 0};
835        Float_t M2[9] = {0, -1, -10.6, -15.005, -31.0635, -39.4988, -60.5011, -46.3992, 0};
836        Float_t M3[9] = {0, 1, 7, 12.5037, 22.8652, 35.2907, 51.4678, 86.4155, 80000};
837    
838        charge = 1000.;
839        mip=0;
840      
841    
842        if (beta<2.) { // it makes no sense to allow beta=5 or so...
843    
844          if( L2->GetNTracks(trkAlg)>=1 ){
845            mip=qNmin1;
846            
847            if (mip>0 && defl<0.7 && defl>0)   {
848              
849              Float_t ym = mip;
850              Float_t xb = defl;
851              
852              for ( Int_t jj=0; jj<n1; jj++ ){
853                xx1[jj] = M1[jj]*xb*xb+M2[jj]*xb+M3[jj];
854                yy1[jj] = L0[jj]*L0[jj] ;
855              }
856              
857              TGraph *gr1 = new TGraph(n1,xx1,yy1);
858              TSpline3 *spl1 = new TSpline3("grs",gr1);    // use a cubic spline
859              Float_t chelp = spl1->Eval(ym);
860              charge = TMath::Sqrt(chelp);
861              gr1->Delete();
862              spl1->Delete();
863              
864            } // if (mip1>0 && defl<0.5 && defl>0)
865          }//Ntrack>=1
866        } //beta < 100
867        
868        ZCalo_Nmin1_defl = charge;
869    
870    //=======================  end charge Rome: Nmin1 vs defl  ===========================
871    
872    
873    
874    
875    
876    //    //
877    if ( debug ) this->Print();    if ( debug ) this->Print();
878    if ( debug ) printf(" esci \n");    if ( debug ) printf(" esci \n");
Legend:
Removed from v.1.11  
changed lines
  Added in v.1.20
  ViewVC Help
Powered by ViewVC 1.1.23