// define tracking algorithm
  const char* alg = "NUC"; // or "EXT" "STD" "EXTF" "NUC" "NUCEXT" "NUCEXTF"

//....snip ....

// define PamLevel2  and ToFNuclei
  PamLevel2*  event  = new PamLevel2(dir,filelist,"+AUTO");       // << create pamela event
  ToFNuclei *tofn = new ToFNuclei(event,alg);


// ...snip....

// the n-sigma cuts used for ToFNuclei for the B/C analysis
Float_t nsigma_s12=1.5;  //
Float_t nsigma_s2=2.0; //
Float_t nsigma_s3=1.4; //

 Float_t *charge_layer_trk_raw = NULL;
 Float_t zs[9];  // charge in the 6 layers plus 3 mean charge

 char htitle[50],info[50];

  //==============================================================
  // Here I  read nuclei calibration file
  //==============================================================
  cout<<"reading  nuclei_meansigma  calibration file "<<endl;

sprintf(info,"nuclei_meansigma_10th_1.dat");

cout<<"reading meansigma file "<<info<<endl;

ifstream fin1(info);

  // 5 species:  Li, Be, B, C, O

  TF1 *mean0[9][5];
  TF1 *mean1[9][5];
  TF1 *sigma[9][5];

  for(Int_t i=0;i<9;i++) {
    for(Int_t j=0;j<5;j++) {
      sprintf(htitle, "mean0_%d_%d",j+5,i);
      mean0[i][j] = new TF1(htitle,"pol6",-0.5,1.3);
      sprintf(htitle, "mean1_%d_%d",j+5,i);
      mean1[i][j] = new TF1(htitle,"pol1",1.2,3.);
      sprintf(htitle, "sigma_%d_%d",j+5,i);
      sigma[i][j] = new TF1(htitle,"pol1",-0.5,3.);
    }
  }

  Float_t p0,p1,p2,p3,p4,p5,p6,p0a,p1a;
  Int_t izin,ilayin;


//  S11...S32,S1,S2,S3
// Z= 3,4,5,6,8

for(int i=0;i<9;i++) {
for(int j=0;j<5;j++) {
fin1>>ilayin>>izin>>p0>>p1>>p2>>p3>>p4>>p5>>p6>>p0a>>p1a;
cout<<ilayin<<" "<<izin<<" "<<p0<<" "<<p1<<" "<<p2<<" "<<p3<<" "<<p4<<" "<<p5<<" "<<p6<<" "<<p0a<<" "<<p1a<<endl;
mean0[i][j]->SetParameters(p0,p1,p2,p3,p4,p5,p6);
mean1[i][j]->SetParameters(p0a,p1a);
fin1>>ilayin>>izin>>p0>>p1;
cout<<ilayin<<" "<<izin<<" "<<p0<<" "<<p1<<endl;
sigma[i][j]->SetParameters(p0,p1);
                     }
                     }

fin1.close();



// ....... snip.....

  Long64_t nevents = event->GetEntries();
  cout<<" # events "<<nevents<<endl;

  for(Long64_t iev=0; iev<nevents; iev++) {
    event->GetEntry(iev);

    int   npt = event->GetNTracks(alg); 

    if( npt>0 ) { //if n.physical tracks>0

    PamTrack *pamtrack = event->GetTrack(0,alg); //the method selects the 0th physical track, discarding the image	

    ExtTrack* trktrack = pamtrack->GetExtTrack();
    ToFTrkVar* toftrack = pamtrack->GetToFTrack(); 

//========================================================================
//=======================  only  single tracks    ========================
//========================================================================

     if(
       npt ==1 && //the tracks is single
         trktrack->chi2 > 0 && //the fit converged
         trktrack->nstep < 100 &&    
    true){


      chi2=trktrack->chi2;
      defl = trktrack->GetDeflection();
      rig = 1./defl;
      beta = toftrack->CalcBeta(10.,10.,20.);       //beta12

            if(
                chi2 > 0    &&
               trktrack->GetNX() >= 4  &&
               trktrack->GetNY() >= 3  &&
                defl > 0 &&
                defl < 10 &&
                beta_mean>0.2 &&
                beta_mean<2.0 &&
                true
                ) {

//..........snip........
//=============================================================
//==========  ToFNuclei ===== =================================
//=============================================================

  charge_layer_trk_raw = tofn->Get_Charge_ToF_trk_layer_raw();
  for(Int_t j=0;j<6;j++) {
    zs[j]=charge_layer_trk_raw[j];  // S11,S12,...,S32
  }

// the three mean charges for <S1>, <S2>, <S3>
  zs[6] = 1000.0;
  if(zs[0]<1000 && zs[1]<1000)
    zs[6]=(zs[0]+zs[1])/2;
  else if(zs[0]<1000 && zs[1]>=1000)
    zs[6]=zs[0];
  else if(zs[0]>=1000 && zs[1]<1000)
    zs[6]=zs[1];

  zs[7] = 1000.0;
  if(zs[2]<1000 && zs[3]<1000)
    zs[7]=(zs[2]+zs[3])/2;
  else if(zs[2]<1000 && zs[3]>=1000)
    zs[7]=zs[2];
  else if(zs[2]>=1000 && zs[3]<1000)
    zs[7]=zs[3];

  zs[8] = 1000.0;
  if(zs[4]<1000 && zs[5]<1000)
    zs[8]=(zs[4]+zs[5])/2;
  else if(zs[4]<1000 && zs[5]>=1000)
    zs[8]=zs[4];
  else if(zs[4]>=1000 && zs[5]<1000)
    zs[8]=zs[5];

//================================================================
//============ Defining the cut for the raw selection=============
//============ Define charge using selection criteria  ===========
//================================================================

//===================  S12 & S2mean & S3mean  ====================

Int_t iztof1=0;
Int_t icount=0;

for (Int_t ij=0; ij<5; ij++) {  // Li, Be, B, C, O
  if(
     (rig<=pow(10,1.25) &&
      fabs(zs[1]-mean0[1][ij]->Eval(log10(rig)))< (nsigma_s12*sigma[1][ij]->Eval(log10(rig))) &&
      fabs(zs[7]-mean0[7][ij]->Eval(log10(rig)))< (nsigma_s2*sigma[7][ij]->Eval(log10(rig))) &&
      fabs(zs[8]-mean0[8][ij]->Eval(log10(rig)))< (nsigma_s3*sigma[8][ij]->Eval(log10(rig)))
      ) ||
     (rig>pow(10,1.25) &&
      fabs(zs[1]-mean1[1][ij]->Eval(log10(rig)))< (nsigma_s12*sigma[1][ij]->Eval(log10(rig))) &&
      fabs(zs[7]-mean1[7][ij]->Eval(log10(rig)))< (nsigma_s2*sigma[7][ij]->Eval(log10(rig))) &&
      fabs(zs[8]-mean1[8][ij]->Eval(log10(rig)))< (nsigma_s3*sigma[8][ij]->Eval(log10(rig)))
      )
     ){
   icount++;
   iztof1=ij+3;
   if (ij==4) iztof1 = 8;
        }
                      }

if (icount>1)  cout<<"overlap "<<endl;

//=====================  S12 & S2mean  =========================

Int_t iztof2=0;
icount=0;

for (Int_t ij=0; ij<5; ij++) {  // Li, Be, B, C, O
  if(
     (rig<=pow(10,1.25) &&
      fabs(zs[1]-mean0[1][ij]->Eval(log10(rig)))< (nsigma_s12*sigma[1][ij]->Eval(log10(rig))) &&
      fabs(zs[7]-mean0[7][ij]->Eval(log10(rig)))< (nsigma_s2*sigma[7][ij]->Eval(log10(rig)))
      ) ||
     (rig>pow(10,1.25) &&
      fabs(zs[1]-mean1[1][ij]->Eval(log10(rig)))< (nsigma_s12*sigma[1][ij]->Eval(log10(rig))) &&
      fabs(zs[7]-mean1[7][ij]->Eval(log10(rig)))< (nsigma_s2*sigma[7][ij]->Eval(log10(rig)))
      )
     ){
   icount++;
   iztof2=ij+3;
   if (ij==4) iztof2 = 8;
        }
                      }

// .............snip.......

     } // tracking cuts...
    } // fit converged
   }  // ntracks>1
  }   // event loop