CaloNuclei/src/CaloNuclei.cpp

#include <CaloNuclei.h>
#include <TGraph.h>
#include <TSpline.h>
#include <TMVA/TSpline1.h>

//--------------------------------------
/**
 * Default constructor 
 */
// CaloNuclei::CaloNuclei(){
//   Clear();
// };

CaloNuclei::CaloNuclei(PamLevel2 *l2p,const char* alg){  
  //
  Clear();
  //
  L2 = l2p;
  //
  if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
  //
  OBT = 0;
  PKT = 0;
  atime = 0;
  N = 5;
  R = 3;
  //
  debug = false;
  // debug = true;
  usetrack = true;
  usepl18x = false;
  trkAlg = alg;
  //
};

void CaloNuclei::Clear(){
  //
  UN = 0;
  tr = 0;
  sntr = 0;
  interplane = 0;
  preq = 0.;
  postq = 0.;
  stdedx1 = 0.;
  ethr = 0.;
  dedx1 = 0.;
  dedx3 = 0.;
  qpremean = 0.; 
  qpremeanN = 0.; 
  maxrel = 0;
  qNmin1 = 0;
  qNmin1_w = 0;
  charge_siegen1 = 0;
  ZCalo_maxrel_b = 0;
  ZCalo_dedx_b = 0;
  ZCalo_dedx_defl= 0;
  ZCalo_Nmin1_defl= 0;
  //
  multhit = false;
  gap = false;
  //
};

void CaloNuclei::Print(){
  //
  Process();
  //
  printf("========================================================================\n");
  printf(" OBT: %u PKT: %u ATIME: %u Track %i Use track %i \n",OBT,PKT,atime,tr,usetrack);
  printf(" interplane [number of available dE/dx before interaction]:....... %i\n",interplane);
  printf(" ethr [threshold used to determine interplane]:................... %f \n",ethr); 
  printf(" dedx1 [dE/dx from the first calorimeter plane]:.................. %f \n",dedx1);
  printf(" stdedx1 [dE/dx from the first calorimeter plane standalone]:..... %f \n",stdedx1);
  printf(" dedx3 [dE/dx (average) if the first 3 Si planes]:................ %f \n",dedx3);
  printf(" multhit [true if interplane determined by multiple hits]:........ %i \n",multhit);
  printf(" gap [true if interplane determined by a gap]:.................... %i \n",gap);
  printf(" preq [total energy in MIP before the interaction plane]:......... %f \n",preq);
  printf(" postq [total energy in MIP after the interaction plane]:......... %f \n",postq);
  printf(" qpremean [truncated mean using 3 planes and 3 strips]:........... %f \n",qpremean); 
  printf(" N [no of used plane]:............................................ %i \n",N); 
  printf(" R [no strip used per plane ]:.................................... %i \n",R); 
  printf(" qpremeanN [truncated mean using N planes and R strips]:.......... %f \n",qpremeanN); 
  printf(" qNmin1 [truncated mean using N-1 planes and R strips]: .......... %f \n",qNmin1); 
  printf(" maxrel [dE/dx of strip with maximum release (I plane)]:.......... %f \n",maxrel); 
  printf(" ZCalo_maxrel_b [Z from maximum release in I Calo plane vs beta].. %f \n",ZCalo_maxrel_b);
  printf(" ZCalo_dedx_b [Z from dedx in I Calo plane vs beta].. ............ %f \n",ZCalo_dedx_b);
  printf(" ZCalo_dedx_defl [Z from dedx in I Calo plane vs deflection....... %f \n",ZCalo_dedx_defl);
  printf(" ZCalo_Nmin1_defl [Z from truncated mean (N-1 pl) vs deflection].. %f \n",ZCalo_Nmin1_defl);
  printf("========================================================================\n");
  //
};

void CaloNuclei::Delete(){
  Clear();
  //delete this;
};


void CaloNuclei::Process(){
  Process(0);
};

void CaloNuclei::Process(Int_t ntr){
  //  
  if ( !L2 ){
    printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
    printf(" ERROR: CaloNuclei variables not filled \n");
    return;
  };
  //
  Bool_t newentry = false;
  //
  if ( L2->IsORB() ){
    if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime || ntr != sntr ){
      newentry = true;
      OBT = L2->GetOrbitalInfo()->OBT;
      PKT = L2->GetOrbitalInfo()->pkt_num;
      atime = L2->GetOrbitalInfo()->absTime;
      sntr = ntr;
    };
  } else {
    newentry = true;
  };
  //
  if ( !newentry ) return;
  //
  tr = ntr;
  //
  if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
  //
  Clear();
  //
  if ( debug ) printf(" Always calculate stdedx1 \n");
  //
  // Always calculate stdedx1 and maxrel
  //
  Int_t cont=0;
  Int_t view = 0;
  Int_t plane = 0;
  Int_t strip = 0;
  Int_t indx = 0;
  Float_t vfpl[96];
  Int_t stfpl[96];
  memset(vfpl, 0, 96*sizeof(Float_t));
  memset(stfpl, 0, 96*sizeof(Int_t));
  Float_t mip = 0.;
  for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
    //
    mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
    //
    if ( !usepl18x && view==0 && plane==18 ) mip = 0.;
    //
    //
    // put in vfpl vector the energy release on the first plane
    //
    if ( strip != -1 && view == 1 && plane == 0 ) {
      stfpl[indx] = strip;
      vfpl[indx] = mip;
      indx++;
    };
    //
  };
  //
  if ( debug ) printf(" find energy released along the strip of maximum on the first plane and on the two neighbour strips \n");
  //
  // find energy released along the strip of maximum on the first plane and on the two neighbour strips
  //
  if ( indx > 0 ){
    Int_t mindx = (Int_t)TMath::LocMax(indx,vfpl);
    for (Int_t ii=0; ii<indx; ii++){
      if ( stfpl[ii] == stfpl[mindx] ) stdedx1 += vfpl[ii];
      if ( (mindx-1)>=0 && stfpl[ii] == (stfpl[mindx]-1) ) stdedx1 += vfpl[ii];
      if ( (mindx+1)<96 && stfpl[ii] == (stfpl[mindx]+1) ) stdedx1 += vfpl[ii];
      //      if ( (mindx-1)>=0 && stfpl[ii] == stfpl[mindx-1] ) stdedx1 += vfpl[ii];
      //      if ( (mindx+1)<96 && stfpl[ii] == stfpl[mindx+1] ) stdedx1 += vfpl[ii];
    };
  maxrel = vfpl[mindx];
  } else {
    stdedx1 = 0.;
    maxrel = 0.;
  };
  // cout<<stdedx1<<"      "<<maxrel<<"\n";
  //
  if ( debug ) printf(" if ( !usetrack ) return: usetrack %i ntr %i \n",usetrack,ntr);
  // 
  //
  //  if ( !usetrack ) return;
  //
  PamTrack *ptrack = 0;
  CaloTrkVar *track = 0;
  //
  if ( usetrack ){
    if ( ntr >= 0 ){
        ptrack = L2->GetTrack(ntr,trkAlg);
      if ( ptrack ) track = ptrack->GetCaloTrack();
    } else {
      track = L2->GetCaloStoredTrack(ntr);
    };
    //
    if ( !track && ntr >= 0 ){
      printf(" ERROR: cannot find any track! \n");      
      cout << "ERROR: trk.algorythm --> "<<trkAlg<<endl;
      printf(" ERROR: CaloNuclei variables not completely filled \n");
      return;   
    };
  } else {
    if ( ntr >= 0 ){
      if ( debug ) printf(" ERROR: you asked not to use a track but you are looking for track number %i !\n",ntr);
      if ( debug ) printf(" ERROR: CaloNuclei variables not completely filled \n");
      return;       
    };
  };
  //
  //  Float_t defethr = 6. * 0.90;
  Float_t defethr = 6.25; // = (sqrt(9) - 0.5) ** 2.;
  //
  // Calculate dedx1 and dedx3
  //
  for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
    //
    mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
    //
    if ( !usepl18x && view==0 && plane==18 ) mip = 0.;
    //
    if ( ntr >= 0 ){
      //
      if ( strip != -1 && 
           view == 1   && 
           plane == 0  &&
           ( strip == (track->tibar[0][1]-1) || strip == (track->tibar[0][1]-2) || strip == (track->tibar[0][1]) ) 
           && true ){      
        dedx1 += mip;
      };
      if ( strip != -1 && 
           (( view == 1 && ( plane == 0 || plane == 1 ) ) ||
            ( view == 0 && plane == 0 )) &&
           (( view == 0 && ( strip == track->tibar[0][0] || strip == (track->tibar[0][0]-1) || strip == (track->tibar[0][0]-2) )) ||
            ( view == 1 && ( strip == track->tibar[0][1] || strip == (track->tibar[0][1]-1) || strip == (track->tibar[0][1]-2) )) ||
            ( view == 1 && ( strip == track->tibar[1][1] || strip == (track->tibar[1][1]-1) || strip == (track->tibar[1][1]-2) ))) &&
           true ){
        dedx3 += mip;
      };
    } else {
      //
      if ( strip != -1 && 
           view == 1   && 
           plane == 0  &&
           ( strip == (L2->GetCaloLevel2()->cibar[0][1]-1) || strip == (L2->GetCaloLevel2()->cibar[0][1]-2) || strip == (L2->GetCaloLevel2()->cibar[0][1]) ) 
           && true ){      
        dedx1 += mip;
      };
      if ( strip != -1 && 
           (( view == 1 && ( plane == 0 || plane == 1 ) ) ||
            ( view == 0 && plane == 0 )) &&
           (( view == 0 && ( strip == L2->GetCaloLevel2()->cibar[0][0] || strip == (L2->GetCaloLevel2()->cibar[0][0]-1) || strip == (L2->GetCaloLevel2()->cibar[0][0]-2) )) ||
            ( view == 1 && ( strip == L2->GetCaloLevel2()->cibar[0][1] || strip == (L2->GetCaloLevel2()->cibar[0][1]-1) || strip == (L2->GetCaloLevel2()->cibar[0][1]-2) )) ||
            ( view == 1 && ( strip == L2->GetCaloLevel2()->cibar[1][1] || strip == (L2->GetCaloLevel2()->cibar[1][1]-1) || strip == (L2->GetCaloLevel2()->cibar[1][1]-2) ))) &&
           true ){
        dedx3 += mip;
      };
    };
    //    
  };
  //
  dedx3 /= 3.;
  //  Float_t mesethr = dedx1 * 0.90;
  Float_t mesethr = 0.;
  if ( dedx1 > 0. ) mesethr = (sqrt(dedx1) - 0.50)*(sqrt(dedx1) - 0.50);
  Bool_t aldone = false;
  //
 retry:
  // 
  if ( debug ) printf("retry\n");
  //
  interplane = 0;
  //
  ethr = TMath::Max(defethr,mesethr);
  //
  // Find the interaction plane "interplane"
  //
  Int_t gapth = 3;
  Int_t nhit[2] = {0,0};
  Int_t splane[2] = {-1,-1};
  Int_t sview[2] = {-1,-1};
  Int_t interpl[2] = {-1,-1};
  Int_t interv[2] = {-1,-1};
  Bool_t wmulthit[2] = {false,false};
  Bool_t wgap[2] = {false,false};
  Int_t ii = 0;
  while ( ii<L2->GetCaloLevel1()->istrip ){
    //
    mip = L2->GetCaloLevel1()->DecodeEstrip(ii,view,plane,strip);    
    //
    if ( !usepl18x && view==0 && plane==18 ) mip = 0.;
    //
    //
    if ( ntr >= 0 ){
      if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] &&
           ( strip == (track->tibar[plane][view]-1) || strip == (track->tibar[plane][view]-2) || strip == (track->tibar[plane][view]) ) 
           && true ){      
        if ( debug ) printf(" inside loop: ii %i mip %f view %i plane %i strip %i tibar %i nhit %i splane %i sview %i \n",ii,mip,view,plane,strip,track->tibar[plane][view]-1,nhit[view],splane[view],sview[view]);
        interpl[view] = plane;
        interv[view] = view;
        if ( splane[view] != plane || sview[view] != view ){
          if ( nhit[view] > 1 ){
            wmulthit[view] = true;
            //    if ( splane[view] == -1 ) splane[view] = 0; //
            //    if ( sview[view] == -1 ) sview[view] = view; //
            interpl[view] = splane[view];
            interv[view] = sview[view];   
        };
          if ( plane > splane[view]+gapth ){
            wgap[view] = true;
            //    if ( splane[view] == -1 ) splane[view] = 0;//
            //    if ( sview[view] == -1 ) sview[view] = view; //
            interpl[view] = splane[view];
            interv[view] = sview[view];
          };
          splane[view] = plane;
          sview[view] = view;
          nhit[view] = 1;
        } else {
          nhit[view]++;
        };
      };
    } else {
      if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] &&
           ( strip == (L2->GetCaloLevel2()->cibar[plane][view]-1) || strip == (L2->GetCaloLevel2()->cibar[plane][view]-2) || strip == (L2->GetCaloLevel2()->cibar[plane][view]) ) 
           && true ){      
        if ( debug ) printf(" inside loop: ii %i mip %f view %i plane %i strip %i cibar %i nhit %i splane %i sview %i \n",ii,mip,view,plane,strip,L2->GetCaloLevel2()->cibar[plane][view]-1,nhit[view],splane[view],sview[view]);
        interpl[view] = plane;
        interv[view] = view;
        if ( splane[view] != plane || sview[view] != view ){
          if ( nhit[view] > 1 ){
            wmulthit[view] = true;
            //    if ( splane[view] == -1 ) splane[view] = 0; //
            //    if ( sview[view] == -1 ) sview[view] = view; //
            interpl[view] = splane[view];
            interv[view] = sview[view];   
        };
          if ( plane > splane[view]+gapth ){
            wgap[view] = true;
            //    if ( splane[view] == -1 ) splane[view] = 0;//
            //    if ( sview[view] == -1 ) sview[view] = view; //
            interpl[view] = splane[view];
            interv[view] = sview[view];
          };
          splane[view] = plane;
          sview[view] = view;
          nhit[view] = 1;
        } else {
          nhit[view]++;
        };
      };
    };
    //
    ii++;
    //    
  };
  //
  if (debug ) printf("conversion interpl %i interv %i multhit %i interplane %i \n",interpl[0],interv[0],multhit,interplane);
  Int_t winterplane[2] = {-1,-1};
  //
  for ( Int_t view = 0; view < 2; view++){
    //
    if ( nhit[view] > 1 && !wmulthit[view] && !wgap[view] ){
      wmulthit[view] = true;
      interpl[view] = splane[view];
      interv[view] = sview[view];
    };
    //
    if ( wmulthit[view] ) multhit = true;
    if ( wgap[view] ) gap = true;
    //
    // convert view and plane number of interaction plane into number of available dE/dx measurements before the interaction plane
    //    
    if ( interpl[view] >= 0 ) {
      if ( interv[view] == 0 ){
        winterplane[view] = (1 + interpl[view]) * 2;
      } else {
        winterplane[view] = (1 + interpl[view]) + (1 + interpl[view] - 1);
      };
      if ( wmulthit[view] ) winterplane[view]--;
    };
  };
  if ( winterplane[0] > 0 && winterplane[1] > 0 ){
    if ( multhit ){
      interplane = TMath::Min(winterplane[0],winterplane[1]);
    } else {
      interplane = TMath::Max(winterplane[0],winterplane[1]);
    };
  } else {
    if ( !winterplane[0] || !winterplane[1] ){
      interplane = 0;
    } else {
      interplane = TMath::Max(winterplane[0],winterplane[1]);
    };
  };
  //  
  if ( debug ) printf("2conversion interpl %i interv %i multhit %i interplane %i \n",interpl[1],interv[1],multhit,interplane);
  if ( debug ) printf("3conversion winterpl0 %i winterpl1 %i \n",winterplane[0],winterplane[1]);
  //
  Int_t ipl = 0;
  if ( interplane > 0 ){
    //
    // Calculate preq, postq, qpremean
    //
    cont++;
    ii = 0;
    Int_t ind = -1;
    Int_t qsplane = -1;
    Int_t qsview = -1;
    Int_t ind2 = -1;
    Int_t qsplane2 = -1;
    Int_t qsview2 = -1;
    Float_t qme[200];
    memset(qme,0,200*sizeof(Float_t));
    Float_t qme2[2112];
    memset(qme2,0,2112*sizeof(Float_t));
    //
    while ( ii<L2->GetCaloLevel1()->istrip ){
      //
      mip = L2->GetCaloLevel1()->DecodeEstrip(ii,view,plane,strip);    
      //
      if ( !usepl18x && view==0 && plane==18 ) mip = 0.;
      //
      //
      if ( strip != -1 ){
        if ( view == 0 ){
          ipl = (1 + plane) * 2;
        } else {
          ipl = (1 + plane) + (1 + plane - 1 );
        };
        if ( ipl > interplane ){
          postq += mip;
        } else {
          preq += mip;
          if ( ntr >= 0 ){
            if (  strip == (track->tibar[plane][view]-1) || strip == (track->tibar[plane][view]-2) || strip == (track->tibar[plane][view]) ){
              if ( qsplane != plane || qsview != view ){
                qsplane = plane;
                qsview = view;
                ind++;
                if ( debug && ind > 199 ) printf(" AAAGH!! \n");
                qme[ind] = 0.;
              };
              qme[ind] += mip;
            };
            for ( Int_t ns = 0; ns < R ; ns++){
              Int_t ms =  track->tibar[plane][view] - 1 - ns + (R - 1)/2;
              if ( strip == ms ){
                if ( qsplane2 != plane || qsview2 != view ){
                  qsplane2 = plane;
                  qsview2 = view;
                  ind2++;
                  if ( debug && ind2 > 2112 ) printf(" AAAGH!! \n");
                  qme2[ind2] = 0.;
                };
                qme2[ind2] += mip;
              };
            };  
          } else {
            if (  strip == (L2->GetCaloLevel2()->cibar[plane][view]-1) || strip == (L2->GetCaloLevel2()->cibar[plane][view]-2) || strip == (L2->GetCaloLevel2()->cibar[plane][view]) ){
              if ( qsplane != plane || qsview != view ){
                qsplane = plane;
                qsview = view;
                ind++;
                if ( debug && ind > 199 ) printf(" AAAGH!! \n");
                qme[ind] = 0.;
              };
              qme[ind] += mip;
            };
            for ( Int_t ns = 0; ns < R ; ns++){
              Int_t ms =  L2->GetCaloLevel2()->cibar[plane][view] - 1 - ns + (R - 1)/2;
              if ( strip == ms ){
                if ( qsplane2 != plane || qsview2 != view ){
                  qsplane2 = plane;
                  qsview2 = view;
                  ind2++;
                  if ( debug && ind2 > 2112 ) printf(" AAAGH!! \n");
                  qme2[ind2] = 0.;
                };
                qme2[ind2] += mip;
              };
            };  
          };
        };      
        //
      };
      //
      ii++;
      //    
    };


 //
    // here we must calculate qpremean, order vector qme, select 3 lowest measurements and caculate the mean...
    //    
    if ( debug ){
      for (Int_t l=0; l < interplane; l++){
        printf(" qme[%i] = %f \n",l,qme[l]);  
      };
    };
    //
    Long64_t work[200];
    ind = 0;
    Int_t l = 0;
    Int_t RN = 0;
    Float_t qm = 0.;
    Float_t qm2 = 0.;
    //
    Float_t qmt = ethr*0.8; // *0.9
    //
    Int_t uplim = TMath::Max(3,N);
    Int_t uplim2 = interplane-1;
    //
    while ( l < uplim && ind < interplane ){
      qm = TMath::KOrdStat((Long64_t)interplane,qme,(Long64_t)ind,work);
      if ( qm >= qmt ){ 
        if ( l < 3 ){
          qpremean += qm;
          RN++;
        };
        l++;
        if ( debug ) printf(" value no %i qm %f qmt %f \n",l,qm,qmt); 
      };
      ind++;
    };
    //
    qpremean /= (Float_t)RN;
    ind = 0;
    l = 0;
    RN = 0;
    while ( l < uplim && ind < interplane ){
      qm2 = TMath::KOrdStat((Long64_t)interplane,qme2,(Long64_t)ind,work);
      if ( qm2 >= qmt ){        
        if ( l < N ){
          qpremeanN += qm2;
          RN++;
        };
        l++;
        if ( debug ) printf(" qm2 value no %i qm %f qmt %f RN %i \n",l,qm2,qmt,RN); 
      };
      ind++;
    };
    ////////////////////////////////////
    //to calculate qNmin1///////////////
    ///////////////////////////////////
    //values under threshold
    qm2=0;
    ind = 0;
    l = 0;
    RN = 0;
    S2=0;
    while ( l < uplim2 && ind<interplane){
      qm2 = TMath::KOrdStat((Long64_t)interplane,qme2,(Long64_t)ind,work);
      if ( qm2 < qmt ) S2++;
      ind++;
    }
    qm2=0;
    ind = 0;
    l = 0;
    RN = 0;
    while ( l < uplim2 && ind < interplane ){
      qm2 = TMath::KOrdStat((Long64_t)interplane,qme2,(Long64_t)ind,work);
      if ( qm2 >= qmt ){        
        if ( l < (interplane - 1 - S2)){  
          qNmin1_w += qm2;
          RN++;
        };
        l++;
        if ( debug ) printf(" qm2 value no %i qm %f qmt %f RN %i \n",l,qm2,qmt,RN); 
      };
      ind++;
    };
    qpremeanN /= (Float_t)RN;
    qNmin1_w /= (Float_t)RN;
    UN = RN;
    ///////set qNmin1 definition///////////
    if (interplane==1 || interplane==2){  
      if (dedx1>0) qNmin1=dedx1;
      else if (stdedx1>0) qNmin1=stdedx1;              
    }
    else if (interplane > 2){
      qNmin1 =  qNmin1_w;
    }
    ////////////////////////////////////
    //////////////////////////////////
    //
    if ( debug ) printf(" charge is %f \n",sqrt(qpremean));
    //
    if ( mesethr != ethr && interplane >= 3 && !aldone ){
      Float_t mesethr2 = (sqrt(qpremean) - 0.50)*(sqrt(qpremean) - 0.50);
      if ( mesethr2 < mesethr*0.90 ){
        mesethr = (sqrt(dedx1) - 0.25)*(sqrt(dedx1) - 0.25);
      } else {
        mesethr = mesethr2;
      };
      aldone = true;
      if ( mesethr > defethr ){
        interplane = 0;
        preq = 0.;
        postq = 0.;
        qpremean = 0.; 
        qpremeanN = 0.; 
        qNmin1 = 0;
        multhit = false;
        gap = false;
        goto retry;
      };
    };
  };


//=======================================================================
//===========    charge determination stdedx1 vs. beta    ===============
//======================     Siegen method    ===========================
//=======================================================================

// Data from file Calo_Bands_New_7.dat
Float_t C0[9] = {0 , 1 , 2 , 3 , 4 , 5 , 6 , 8 , 90 };
Float_t B0[9] = {0 , -2.03769 , 7.61781 , 19.7098 , 60.5598 , 57.9226 , 14.8368 , -1358.83 , 8200 };
Float_t B1[9] = {0 , 0.0211274 , 9.32057e-010 , 4.47241e-07 , 1.44826e-06 , 2.6189e-05 , 0.00278178 , 55.5445 , 0 };
Float_t B2[9] = {0 , -3.91742 , -20.0359 , -16.3043 , -16.9471 , -14.4622 , -10.9594 , -2.38014 , 0 };
Float_t B3[9] = {0 , 11.1469 , -6.63105 , -27.8834 , -132.044 , -55.341 , 173.25 , 4115 , 0 };
Float_t B4[9] = {0 , -14.3465 , -0.485215 , 18.8122 , 117.533 , -14.0898 , -325.269 , -4388.89 , 0 };
Float_t B5[9] = {0 , 6.24281 , 3.96018 , 0 , -26.1881 , 42.9731 , 182.697 , 1661.01 , 0 };

Float_t x1[9],y1[9];
Int_t n1 = 9;

Float_t charge = 1000.;
Float_t beta = 100.;

//------- First try track dependent beta
    if( L2->GetNTracks(trkAlg)>=1 ){
        PamTrack *TRKtrack = L2->GetTrack(0,trkAlg);
    if (fabs(TRKtrack->GetToFTrack()->beta[12]) < 100.) beta = fabs(TRKtrack->GetToFTrack()->beta[12]);
                                                }
//------- If no beta found, try standalone beta
    if (beta == 100.) {
    ToFTrkVar *ttrack = L2->GetToFStoredTrack(-1);
    beta = fabs(ttrack->beta[12]);
                      }

    if (beta<2.) {  // it makes no sense to allow beta=5 or so...

    Float_t  mip=0;
    mip=stdedx1 ;

    if (mip>0)   {

    Float_t betahelp =   pow(beta, 1.8);
    Float_t ym = mip*betahelp;
    Float_t xb = beta;

    for ( Int_t jj=0; jj<9; jj++ ){
    x1[jj] = B0[jj]+B1[jj]*pow(xb,B2[jj])+B3[jj]*xb+B4[jj]*xb*xb+B5[jj]*xb*xb*xb;
    y1[jj] = C0[jj]*C0[jj] ;
                                  }

    TGraph *gr1 = new TGraph(n1,x1,y1);
    TSpline3 *spl1 = new TSpline3("grs",gr1);    // use a cubic spline
    Float_t chelp = spl1->Eval(ym);
    charge = TMath::Sqrt(chelp);
    gr1->Delete();
    spl1->Delete();

                    } // if (mip1>0)
                    } // beta < 100


    charge_siegen1 = charge;

//=======================  end charge Siegen  ===========================


// //=======================================================================
// //===========    charge determination Maximum release  vs. beta    ===============
// //======================     Rome method    ===========================
// //=======================================================================     

    Float_t D0[9] = {0, 1, 2, 3 , 4 , 5 , 6, 8, 90};
    Float_t E1[9] = {0, 500, 500, 923.553 , 659.842, 1113.97, 3037.25, 3034.84, 0};
    Float_t E2[9] = {0, 11.0, 7.5, 6.92574 , 5.08865, 5.29349, 6.41442, 5.52969, 0};
    Float_t E3[9] = {0, 1.2, 4, 9.7227  , 13.18, 23.5444, 38.2057, 63.6784, 80000};

    Float_t xx1[9],yy1[9];
    n1 = 9;
    
    charge = 1000.;
    mip=0;
    

    if (beta<2.) {  // it makes no sense to allow beta=5 or so...
      
      
      mip=maxrel;

      if (mip>0)   {
        Float_t ym = mip;
        Float_t xb = beta;
        
        for ( Int_t jj=0; jj<n1; jj++ ){
          xx1[jj] = E1[jj]*exp(-E2[jj]*xb)+E3[jj];
          yy1[jj] = D0[jj]*D0[jj] ;
        }
        
        TGraph *gr1 = new TGraph(n1,xx1,yy1);
        TSpline3 *spl1 = new TSpline3("grs",gr1);    // use a cubic spline
        Float_t chelp = spl1->Eval(ym);
        charge = TMath::Sqrt(chelp);
        gr1->Delete();
        spl1->Delete();

        
      } // if (mip1>0)
    } // beta < 100


    ZCalo_maxrel_b = charge;

    //=======================  end charge Rome: maxril vs beta  ===========================


// =======================================================================
// ===========    charge determination dedx  vs. beta    ===============
// ======================     Rome method    ===========================
// =======================================================================   

    Float_t F0[9] = {0., 1., 2., 3. ,4., 5. , 6., 8, 90};
    Float_t G1[9] = {0, 500, 500, 642.935 , 848.684, 1346.05, 3238.82, 3468.6, 0};
    Float_t G2[9] = {0, 11, 7.5, 6.2038 , 5.51723, 5.65265, 6.54089, 5.72723, 0};
    Float_t G3[9] = {0, 1.2, 4, 9.2421 , 13.9858, 25.3912, 39.6332, 64.5674, 80000};

  
    charge = 1000.;
    mip=0;


    if (beta<2.) { // it makes no sense to allow beta=5 or so...
      
    
      if( L2->GetNTracks(trkAlg)>=1 ){
        mip=dedx1;
      } 
      if (mip==0) mip=stdedx1;
      

      if (mip>0)   {
        
        Float_t ym = mip;
        Float_t xb = beta;
        
        for ( Int_t jj=0; jj<n1; jj++ ){
         xx1[jj] = G1[jj]*exp(-G2[jj]*xb)+G3[jj];
         yy1[jj] = F0[jj]*F0[jj] ;
        }
        
        TGraph *gr1 = new TGraph(n1,xx1,yy1);
        TSpline3 *spl1 = new TSpline3("grs",gr1);    // use a cubic spline
        Float_t chelp = spl1->Eval(ym);
        charge = TMath::Sqrt(chelp);
        gr1->Delete();
        spl1->Delete();
        
      } //if (mip1>0)
    } //beta < 100
    
    ZCalo_dedx_b = charge;

    //=======================  end charge Rome: dedx vs beta  ===========================
    
    
//=======================================================================
//===========    charge determination dedx  vs. defl    ===============
//======================     Rome method    ===========================
//=======================================================================     

     Float_t H0[9] = {0, 1, 2, 3 , 4 , 5 , 6, 8, 90 };
     Float_t I1[9] = {0, 3.5, 40, 56.1019, 101.673, 109.225, 150.599, 388.531, 0};
     Float_t I2[9] = {0, -1, -13.6, -12.5581, -22.5543, -15.9823, -28.2207, -93.6871, 0};
     Float_t I3[9] = {0, 1, 5.3, 11.6218, 19.664, 32.1817, 45.7527, 84.5992, 80000};

    
    charge = 1000.;
    mip=0;
    Float_t defl=0;
    

    if (beta<2.) { // it makes no sense to allow beta=5 or so...

      if( L2->GetNTracks(trkAlg)>=1 ){
          PamTrack *TRKtrack = L2->GetTrack(0,trkAlg);
        mip=dedx1;
        if (mip==0) mip=stdedx1;
        defl=TRKtrack->GetTrkTrack()->al[4];
        
        
        if (mip>0 && defl<0.7 && defl>0)   {
          
          Float_t ym = mip;
          Float_t xb = defl;
          
          for ( Int_t jj=0; jj<n1; jj++ ){
            xx1[jj] = I1[jj]*xb*xb+I2[jj]*xb+I3[jj];
            yy1[jj] = H0[jj]*H0[jj] ;
          }
          
          TGraph *gr1 = new TGraph(n1,xx1,yy1);
          TSpline3 *spl1 = new TSpline3("grs",gr1);    // use a cubic spline
          Float_t chelp = spl1->Eval(ym);
          charge = TMath::Sqrt(chelp);
          gr1->Delete();
          spl1->Delete();
          
        } // if (mip1>0 && defl<0.5 && defl>0)
      }//Ntrack>=1
    } //beta < 100

    ZCalo_dedx_defl = charge;
 
//=======================  end charge Rome: dedx vs defl  ===========================


//============================================================================================
//===========    charge determination Truncated mean (N-1 planes) vs. defl ===================
//================================     Rome method    ========================================
//============================================================================================   

    Float_t L0[9] = {0, 1, 2, 3 , 4 , 5 , 6, 8, 90};
    Float_t M1[9] = {0, 3.5, 27, 63.0145, 120.504, 173.663, 245.33, 236.517, 0};
    Float_t M2[9] = {0, -1, -10.6, -15.005, -31.0635, -39.4988, -60.5011, -46.3992, 0};
    Float_t M3[9] = {0, 1, 7, 12.5037, 22.8652, 35.2907, 51.4678, 86.4155, 80000}; 

    charge = 1000.;
    mip=0;
   

    if (beta<2.) { // it makes no sense to allow beta=5 or so...

      if( L2->GetNTracks(trkAlg)>=1 ){
        mip=qNmin1;
        
        if (mip>0 && defl<0.7 && defl>0)   {
          
          Float_t ym = mip;
          Float_t xb = defl;
          
          for ( Int_t jj=0; jj<n1; jj++ ){
            xx1[jj] = M1[jj]*xb*xb+M2[jj]*xb+M3[jj];
            yy1[jj] = L0[jj]*L0[jj] ;
          }
          
          TGraph *gr1 = new TGraph(n1,xx1,yy1);
          TSpline3 *spl1 = new TSpline3("grs",gr1);    // use a cubic spline
          Float_t chelp = spl1->Eval(ym);
          charge = TMath::Sqrt(chelp);
          gr1->Delete();
          spl1->Delete();
          
        } // if (mip1>0 && defl<0.5 && defl>0)
      }//Ntrack>=1
    } //beta < 100
    
    ZCalo_Nmin1_defl = charge;
 
//=======================  end charge Rome: Nmin1 vs defl  ===========================


  //
  if ( debug ) this->Print();
  if ( debug ) printf(" esci \n");
  //
};