myEvent.cxx

Go to the documentation of this file.

// Authors: M.Duranti - INFN di Perugia
#include "myEvent.h"
#include "debug.h"
#include <TClass.h>
#include <iostream>
#include <config.hxx>
#include "myChain.h"

using namespace std;

//--------------------------------------------------------------------

ClassImp(myEvent);

//--------------------------------------------------------------------

myEvent* myEvent::ptr=0;
section_t gsection;

Double_t myEvent::zlay[9];

unsigned int TimePrevSlow = 0;
unsigned int TimePrevISS = 0;

myEvent::myEvent():
  mm(nullptr),ms(nullptr),RTI(nullptr)
{
#ifdef PDEBUG
  printf("In myEvent::myEvent\n");
#endif
  init();
  PRINTDEBUG;
  if (ptr) {
    printf("%s::%s was already existing (%p)... Substituting it (with %p)!\n", ptr->ClassName(), ptr->GetName(), ptr, this);
  }
  else printf("Object myEvent (%p) created...\n", this);//only for debug
    ptr = this;
}


myEvent::myEvent(const myEvent& orig){

  kMC=orig.kMC;

  if(orig.mm)  mm = new myMC(*orig.mm);
  else mm=0;
  
  if(orig.ms) ms = new myStatus(*orig.ms);
  else ms=0;
  
  fHeader = myHeader(orig.fHeader); 
  
  us1=0;

  vmp.clear();
  for (int ii=0; ii<(int)(vmp.size()); ii++) {
    vmp.push_back(orig.vmp[ii]);
  }
  
  vmt.clear();
  for (int ii=0; ii<(int)(vmt.size()); ii++) {
    vmt.push_back(orig.vmt[ii]);
  }
  vme.clear();
  for (int ii=0; ii<(int)(vme.size()); ii++) {
    vme.push_back(orig.vme[ii]);
  }
  vmu.clear();
  for (int ii=0; ii<(int)(vmu.size()); ii++) {
    vmu.push_back(orig.vmu[ii]);
  }
  vmb.clear();
  for (int ii=0; ii<(int)(vmb.size()); ii++) {
    vmb.push_back(orig.vmb[ii]);
  }
  vmbh.clear();
  for (int ii=0; ii<(int)(vmbh.size()); ii++) {
    vmbh.push_back(orig.vmbh[ii]);
  }
  vmr.clear();
  for (int ii=0; ii<(int)(vmr.size()); ii++) {
    vmr.push_back(orig.vmr[ii]);
  }
  
  vmup.clear();
  for (int ii=0; ii<(int)(vmup.size()); ii++) {
    vmup.push_back(orig.vmup[ii]);
  }
  /*
  vmukp.clear();
  for (int ii=0; ii<(int)(vmukp.size()); ii++) {
    vmukp.push_back(orig.vmukp[ii]);
  }
  */
  vmuktp = orig.vmuktp;
  vmukup = orig.vmukup;

  vmuqtp.clear();
  for (int ii=0; ii<(int)(vmuqtp.size()); ii++) {
    vmuqtp.push_back(orig.vmuqtp[ii]);
  }
  vmtp.clear();
  for (int ii=0; ii<(int)(vmtp.size()); ii++) {
    vmtp.push_back(orig.vmtp[ii]);
  }
  vmep.clear();
  for (int ii=0; ii<(int)(vmep.size()); ii++) {
    vmep.push_back(orig.vmep[ii]);
  }
  vmrp.clear();
  for (int ii=0; ii<(int)(vmrp.size()); ii++) {
    vmrp.push_back(orig.vmrp[ii]);
  }
  vmbp.clear();
  for (int ii=0; ii<(int)(vmbp.size()); ii++) {
    vmbp.push_back(orig.vmbp[ii]);
  }
  vmbhp.clear();
  for (int ii=0; ii<(int)(vmbhp.size()); ii++) {
    vmbhp.push_back(orig.vmbhp[ii]);
  }
  
  // fNmyParticle   =orig.fNmyParticle;

  // fNmyTrTrack    =orig.fNmyTrTrack;     
  // fNmyEcalShower       =orig.fNmyEcalShower;        
  // fNmyTrdTrack        =orig.fNmyTrdTrack;         
  // fNmyBeta       =orig.fNmyBeta;        
  // fNmyBetaH      =orig.fNmyBetaH;       
  // fNmyRichRing       =orig.fNmyRichRing;        
  
  // fNmyTrTrackPlus=orig.fNmyTrTrackPlus; 
  // fNmyEcalShowerPlus   =orig.fNmyEcalShowerPlus;    
  // fNmyTrdKFromTrTrack    =orig.fNmyTrdKFromTrTrack;     
  // fNmyTrdQtFromTrTrack    =orig.fNmyTrdQtFromTrTrack;     
  // fNmyBetaPlus   =orig.fNmyBetaPlus;    
  // fNmyBetaHPlus  =orig.fNmyBetaHPlus;   
  // fNmyRichRingPlus   =orig.fNmyRichRingPlus;    

  if (orig.RTI) RTI= new myRTI(*orig.RTI);
  else RTI=0;
  
  LiveTime=orig.           LiveTime;
  DeltaT=orig.             DeltaT;                 
  Lvl1Time=orig.           Lvl1Time;               
  Lvl1TimePrevious=orig.   Lvl1TimePrevious;       
  DeltaTMine=orig.         DeltaTMine;             
  TofFlag1=orig.           TofFlag1;               
  TofFlag2=orig.           TofFlag2;               
  JMembPatt=orig.          JMembPatt;              
  PhysBPatt=orig.          PhysBPatt;              
  PhysBPattRestored=orig.  PhysBPattRestored;      
  PhysBPattPhys=orig.      PhysBPattPhys;          
  for (int ii=0;ii<4;ii++) TofPatt1[ii]=orig.TofPatt1[ii];            
  for (int ii=0;ii<4;ii++) TofPatt2[ii]=orig.TofPatt2[ii];            
  EcalFlag=orig.           EcalFlag;               
  for (int ii=0;ii<19;ii++) TrigRates[ii]=orig.TrigRates[ii];          
  nAnti=orig.              nAnti;                  
  AntiPatt=orig.           AntiPatt;               
  ReasonToBeBadRun=orig.   ReasonToBeBadRun;       
  _IsBadRun=orig.          _IsBadRun;               
  RunAnalCat=orig.         RunAnalCat;             
  Run=orig.                Run;                    
  Event=orig.              Event;                  
  fHeader=orig.            fHeader;
  DeltaTJMDC=orig.         DeltaTJMDC;             
  TimeJMDCFirst=orig.      TimeJMDCFirst;          
  TimeJMDCLast=orig.       TimeJMDCLast;           
  Duration=orig.           Duration;               
  for (int ii=0;ii<4;ii++) nTofClusterL[ii]=orig.          nTofClusterL[ii];        
  nTofCluster=orig.        nTofCluster;            
  SensorA=orig.            SensorA;                
  SensorC=orig.            SensorC;                
  SensorK=orig.            SensorK;                
  SensorL=orig.            SensorL;                
  SensorM=orig.            SensorM;                
  Plane1N_M70=orig.        Plane1N_M70;            
  Plane1N_M71=orig.        Plane1N_M71;            
  Plane1N_M72=orig.        Plane1N_M72;            
  Plane1N_M73=orig.        Plane1N_M73;            
  Sensor8=orig.            Sensor8;                
  for (int ii=0;ii<3;ii++) GeoMagSphericalCoo[ii]=orig.GeoMagSphericalCoo[ii]; 
  LatPred=orig.            LatPred;                
  LonPred=orig.            LonPred;                
  AltPred=orig.            AltPred;                
  LatCGM=orig.             LatCGM;                 
  LonCGM=orig.             LonCGM;                 
  RadCGM=orig.             RadCGM;                 
  MaxGeoCutOff_p=orig.     MaxGeoCutOff_p;         
  MaxGeoCutOff_n=orig.     MaxGeoCutOff_n;         
  Field=orig.              Field;                  
  Lenght=orig.             Lenght;                 
  Tdr=orig.                Tdr;                    
  Udr=orig.                Udr;                    
  Sdr=orig.                Sdr;                    
  Rdr=orig.                Rdr;                    
  Edr=orig.                Edr;                    
  L1dr=orig.               L1dr;                   
  L3dr=orig.               L3dr;
  L3Error=orig.            L3Error;
  L3VEvent=orig.           L3VEvent;
  L3TimeD=orig.            L3TimeD;
  for (int ii=0;ii<4;ii++) JINJStatus[ii] =orig.JINJStatus[ii];          
  for (int ii=0;ii<24;ii++) JError[ii]    =orig.JError[ii];                 

  return;
}


myEvent::~myEvent(){
#ifdef PDEBUG
  printf("In myEvent::~myEvent\n");
#endif
  Clear();
  PRINTDEBUG;
  ptr=0;
  PRINTDEBUG;
}

void myEvent::init() {
#ifdef PDEBUG
  printf("In myEvent::init\n");
#endif

  Clear();

  return;
}

myEvent* myEvent::gethead() {
#ifdef PDEBUG
  printf("In myEvent::gethead\n");
#endif

  PRINTDEBUG;

  if (!ptr) {
    ptr = new myEvent();
    printf("Object myEvent (%p) created as singleton...\n", ptr);
  }
#ifdef PDEBUG
  else printf("%s existed (%p)!\n", ptr->ClassName(), ptr);//only for debug
  printf("Event=%d\n", ptr->Event);//only for debug
#endif
  return ptr;
}

void myEvent::printptr() {
#ifdef PDEBUG
  printf("In myEvent::printptr\n");
#endif

  PRINTDEBUG;

  printf("myEvent pointer is %p...\n", ptr);

  return;
}
 
void myEvent::Clear(Option_t*){
#ifdef PDEBUG
  printf("In myEvent::Clear\n");
#endif
  PRINTDEBUG;

  kMC=false;
  if (mm) delete mm;
  mm=0;

  PRINTDEBUG;

  fHeader.Clear();

  PRINTDEBUG;

  if (ms) delete ms;
  ms=0;
    
  PRINTDEBUG;

  us1=0;

  PRINTDEBUG;

  vmp.clear();

  vmt.clear();
  vme.clear();
  vmu.clear();
  vmb.clear();
  vmbh.clear();
  vmr.clear();

  vmtp.clear();
  vmep.clear();
  vmup.clear();
  vmuktp.clear();
  vmukup.clear();
  vmuqtp.clear();
  vmbp.clear();
  vmbhp.clear();
  vmrp.clear();

  PRINTDEBUG;

  // fNmyParticle=0;

  // fNmyTrTrack=0;
  // fNmyEcalShower=0;
  // fNmyTrdTrack=0;
  // fNmyBeta=0;
  // fNmyBetaH=0;
  // fNmyRichRing=0;

  // fNmyTrTrackPlus=0;
  // fNmyEcalShowerPlus=0;
  // fNmyTrdKFromTrTrack=0;
  // fNmyTrdQtFromTrTrack=0;
  // fNmyBetaPlus=0;
  // fNmyBetaHPlus=0;
  // fNmyRichRingPlus=0;

  PRINTDEBUG;

  _nParticle = -999999;

  _nTrdTrack = -999999;
  _nTrdCluster = -999999;
  _nTrdRawHit = -999999;
  _nTrdSegment =-999999;

  PRINTDEBUG;

  _nEcalShower = -999999;
  
  PRINTDEBUG;

  _nBeta = -999999;
  _nBetaH = -999999;

  PRINTDEBUG;

  _nRichRing = -999999;
  _nRichHit  = -999999;

  _nTrTrack=-999999;
  _nTrCluster=-999999;
  _nTrRecHit=-999999;
  _nTrRawCluster = -999999;

  _UTCTime=0;

  fill_n(_nTrRawClusterL, 9, 0);

  PRINTDEBUG;

  LiveTime = -999999;
  DeltaT = 0xFFFF; // in musec
  Lvl1Time=-999999;
  //  Lvl1TimePrevious=-999999;//non resettarlo!
  DeltaTMine = -999999; // in musec
  TofFlag1 = -999999;
  TofFlag2 = -999999;
  JMembPatt = -999999;
  PhysBPatt = -999999;
  PhysBPattRestored = -999999;
  PhysBPattPhys = -999999;
  fill_n(TofPatt1,  4, 0);
  fill_n(TofPatt2,  4, 0);
  EcalFlag= -999999;
  fill_n(TrigRates, 19,0);

  PRINTDEBUG;

  nAnti = -999999;
  AntiPatt = 0;  

  PRINTDEBUG;

  ReasonToBeBadRun = "";
  _IsBadRun = -999999;
  RunAnalCat = -999999;

  Run = 0;
  Event = 99999999;

  //  TimeJMDCFirst=0;//non resettarlo!
  //  TimeJMDCLast=0;//non resettarlo!
  //  Duration=-999999;//non resettarlo!
  
  PRINTDEBUG;

  fill_n(nTofClusterL, 4, 0);
  nTofCluster=-999999;

  PRINTDEBUG;

  SensorA=-999999;
  SensorC=-999999;
  SensorK=-999999;
  SensorL=-999999;
  SensorM=-999999;
  Plane1N_M70=-999999;
  Plane1N_M71=-999999;
  Plane1N_M72=-999999;
  Plane1N_M73=-999999;
  Sensor8=-999999;

  PRINTDEBUG;

  fill_n(GeoMagSphericalCoo, 3, 0);

  LatPred=-999999;
  LonPred=-999999;
  AltPred=-999999;
  
  Field=-999999;
  
  LatCGM=-999999;
  LonCGM=-999999;
  RadCGM=-999999;

  // not cleared since evaluated only once per second
  // MaxGeoCutoff_p=-999999;
  // MaxGeoCutOff_n=-999999;
  
  PRINTDEBUG;

  if(RTI) delete RTI;
  RTI=0;

  PRINTDEBUG;
  
  Lenght=999999;
  Tdr=999999;
  Udr=999999;
  Sdr=999999;
  Rdr=999999;
  Edr=999999;
  L1dr=999999;
  L3dr=999999;
  L3Error=999999;
  L3VEvent=999999;
  L3TimeD=999999;
  fill_n(JINJStatus, 4, 0);
  fill_n(JError, 24, 0);

  PRINTDEBUG;

}

int myEvent::TofHSelection(){
  
  int beta_index=-1;
  int n_betah=0;
  if ( nBetaH() == 0) return false; // Skip this event                                                                                                                                                                
  n_betah= nBetaH();

  Int_t    btchk = 0;
  Double_t betat = 0;
  Int_t    max_UseHit=0;
  Int_t    maxPattern=0;
  myBetaH *bth=0;

  for (Int_t i = 0; i < n_betah; i++) {
    myBetaH *bb= pmyBetaH(i);
    if (!bb) continue;
    Double_t bt = bb->GetBeta();
    cout<<" i "<<i
        <<" bt "<<bt
        <<" GetUseHit "<<bb->GetUseHit()
        <<" GetBetaPattern "<<bb->GetBetaPattern() 
        <<endl;
    
    if (bt < betat) continue;

    Int_t chk = 0;
    if (bb->iTrdTrack()>=0  ) chk++;
    if (bb->iEcalShower()>=0) chk++;

    if (bb->GetUseHit()   > max_UseHit) chk++;
    if (bb->GetBetaPattern() == 4444) chk++;

    if(bb->GetUseHit()>= max_UseHit){
      max_UseHit=bb->GetUseHit();
      chk++;
    }
    int pattern=(bb->GetBetaPattern()%10)
      +((bb->GetBetaPattern()/10)%10)
      +((bb->GetBetaPattern()/100)%10)
      +((bb->GetBetaPattern()/1000)%10);
    
    if(pattern>= maxPattern){
      maxPattern=pattern;
      chk++;
    }
    if (chk < btchk) continue;  // Get the highest quality BetaH     
      
    bth   = bb;
    betat = bt;
    btchk = chk;
    beta_index=i;
    cout<<" i "<<i
        <<" bt "<<bt
        <<" btchk "<<btchk
        <<" beta_index "<<beta_index
        <<endl;
  }
  if (!bth) return -1; // Skip this event

  return beta_index;
}

//---------------------------------------------

myParticle* myEvent::pmyParticleByMySelf(int index){
  
  if (index>=0 && vmp.size()>0) {
    for (int ii=0; ii<(int)(vmp.size()); ii++) {
      if (vmp[ii] && ((myObject*)vmp[ii])->iMySelf()==index) {
        return vmp[ii];
      }
    }
  }
  
  return NULL;
}

myParticle* myEvent::pmyParticle(int index){
  
  if (index>=0 && (int)(vmp.size())>index)
    return vmp[index];
  
  return NULL;
}

//---------------------------------------------

myTrdTrack* myEvent::pmyTrdTrackByMySelf(int index){

  if (index>=0 && vmu.size()>0) {
    for (int ii=0; ii<(int)(vmu.size()); ii++) {
      if (vmu[ii] && ((myObject*)vmu[ii])->iMySelf()==index) {
        return vmu[ii];
      }
    }
  }
  
  return NULL;
}

myTrdTrack* myEvent::pmyTrdTrack(int index){
  
  if (index>=0 && (int)(vmu.size())>index)
    return vmu[index];
  
  return NULL;
}

myTrdTrackPlus* myEvent::pmyTrdTrackPlusByMySelf(int index){

  if (index>=0 && vmup.size()>0) {
    for (int ii=0; ii<(int)(vmup.size()); ii++) {
      if (vmup[ii] && ((myObject*)vmup[ii])->iMySelf()==index) {
        return vmup[ii];
      }
    }
  }
  
  return NULL;
}

myTrdTrackPlus* myEvent::pmyTrdTrackPlus(int index){
  
  if (index>=0 && (int)(vmup.size())>index)
    return vmup[index];
  
  return NULL;
}


myTrTrack* myEvent::pmyTrTrackByMySelf(int index){
  
  if (index>=0 && vmt.size()>0) {
    for (int ii=0; ii<(int)(vmt.size()); ii++) {
      if (vmt[ii] && ((myObject*)vmt[ii])->iMySelf()==index) {
        return vmt[ii];
      }
    }
  }
  
  return NULL;
}

myTrTrack* myEvent::pmyTrTrack(int index){
  
  if (index>=0 && (int)(vmt.size())>index)
    return vmt[index];
  
  return NULL;
}

myEcalShower* myEvent::pmyEcalShowerByMySelf(int index){

  if (index>=0 && vme.size()>0) {
    for (int ii=0; ii<(int)(vme.size()); ii++) {
      if (vme[ii] && ((myObject*)vme[ii])->iMySelf()==index) {
        return vme[ii];
      }
    }
  }
  
  return NULL;
}

myEcalShower* myEvent::pmyEcalShower(int index){
  
  if (index>=0 && (int)(vme.size())>index)
    return vme[index];
  
  return NULL;
}

myBeta* myEvent::pmyBetaByMySelf(int index){

  if (index>=0 && vmb.size()>0) {
    for (int ii=0; ii<(int)(vmb.size()); ii++) {
      if (vmb[ii] && ((myObject*)vmb[ii])->iMySelf()==index) {
        return vmb[ii];
      }
    }
  }
  
  return NULL;
}

myBeta* myEvent::pmyBeta(int index){
  
  if (index>=0 && (int)(vmb.size())>index)
    return vmb[index];
  
  return NULL;
}

myBetaH* myEvent::pmyBetaHByMySelf(int index){

  if (index>=0 && vmbh.size()>0) {
    for (int ii=0; ii<(int)(vmbh.size()); ii++) {
      if (vmbh[ii] && ((myObject*)vmbh[ii])->iMySelf()==index) {
        return vmbh[ii];
      }
    }
  }
  
  return NULL;
}

myBetaH* myEvent::pmyBetaH(int index){
  
  if (index>=0 && (int)(vmbh.size())>index)
    return vmbh[index];
  
  return NULL;
}

myRichRing* myEvent::pmyRichRingByMySelf(int index){

  if (index>=0 && vmr.size()>0) {
    for (int ii=0; ii<(int)(vmr.size()); ii++) {
      if (vmr[ii] && ((myObject*)vmr[ii])->iMySelf()==index) {
        return vmr[ii];
      }
    }
  }
  
  return NULL;
}

myRichRing* myEvent::pmyRichRing(int index){
  
  if (index>=0 && (int)(vmr.size())>index)
    return vmr[index];
  
  return NULL;
}

//--------------------

myTrdK* myEvent::pmyTrdKFromTrTrackByMySelf(int index){

  if (index>=0 && vmuktp.size()>0) {
    for (int ii=0; ii<(int)(vmuktp.size()); ii++) {
      if (vmuktp[ii] && ((myObject*)vmuktp[ii])->iMySelf()==index) {
        return vmuktp[ii];
      }
    }
  }
  
  return NULL;
}
myTrdK* myEvent::pmyTrdKFromTrdTrackByMySelf(int index){
  if (index>=0 && vmukup.size()>0) {
    for (size_t ii=0; ii<vmukup.size(); ii++) {
      myTrdK * k = vmukup[ii];
      if (k && k->iMySelf()==index) { return k; }
    }
  }
  return NULL;
}

myTrdK* myEvent::pmyTrdKFromTrTrack(int index){
  
  if (index>=0 && (int)(vmuktp.size())>index)
    return vmuktp[index];
  
  return NULL;
}

myTrdK* myEvent::pmyTrdKFromTrdTrack(int index){

  if (index>=0 && (int)(vmukup.size())>index)
    return vmukup[index];
  
  return NULL;
}

myTrdQtFromTrTrack* myEvent::pmyTrdQtFromTrTrackByMySelf(int index){

  if (index>=0 && vmuqtp.size()>0) {
    for (int ii=0; ii<(int)(vmuqtp.size()); ii++) {
      if (vmuqtp[ii] && ((myObject*)vmuqtp[ii])->iMySelf()==index) {
        return vmuqtp[ii];
      }
    }
  }
  
  return NULL;
}

myTrdQtFromTrTrack* myEvent::pmyTrdQtFromTrTrack(int index){
  
  if (index>=0 && (int)(vmuqtp.size())>index)
    return vmuqtp[index];
  
  return NULL;
}

myTrTrackPlus* myEvent::pmyTrTrackPlusByMySelf(int index){

  if (index>=0 && vmtp.size()>0) {
    for (int ii=0; ii<(int)(vmtp.size()); ii++) {
      if (vmtp[ii] && ((myObject*)vmtp[ii])->iMySelf()==index) {
        return vmtp[ii];
      }
    }
  }
  
  return NULL;
}

myTrTrackPlus* myEvent::pmyTrTrackPlus(int index){
  
  if (index>=0 && (int)(vmtp.size())>index)
    return vmtp[index];
  
  return NULL;
}

myEcalShowerPlus* myEvent::pmyEcalShowerPlusByMySelf(int index){

  if (index>=0 && vmep.size()>0) {
    for (int ii=0; ii<(int)(vmep.size()); ii++) {
      if (vmep[ii] && ((myObject*)vmep[ii])->iMySelf()==index) {
        return vmep[ii];
      }
    }
  }
  
  return NULL;
}

myEcalShowerPlus* myEvent::pmyEcalShowerPlus(int index){
  
  if (index>=0 && (int)(vmep.size())>index)
    return vmep[index];
  
  return NULL;
}

myBetaPlus* myEvent::pmyBetaPlusByMySelf(int index){

  if (index>=0 && vmbp.size()>0) {
    for (int ii=0; ii<(int)(vmbp.size()); ii++) {
      if (vmbp[ii] && ((myObject*)vmbp[ii])->iMySelf()==index) {
        return vmbp[ii];
      }
    }
  }
  
  return NULL;
}

myBetaPlus* myEvent::pmyBetaPlus(int index){
  
  if (index>=0 && (int)(vmbp.size())>index)
    return vmbp[index];
  
  return NULL;
}

myBetaHPlus* myEvent::pmyBetaHPlusByMySelf(int index){

  if (index>=0 && vmbhp.size()>0) {
    for (int ii=0; ii<(int)(vmbhp.size()); ii++) {
      if (vmbhp[ii] && ((myObject*)vmbhp[ii])->iMySelf()==index) {
        return vmbhp[ii];
      }
    }
  }
  
  return NULL;
}

myBetaHPlus* myEvent::pmyBetaHPlus(int index){
  
  if (index>=0 && (int)(vmbhp.size())>index)
    return vmbhp[index];
  
  return NULL;
}

myRichRingPlus* myEvent::pmyRichRingPlusByMySelf(int index){

  if (index>=0 && vmrp.size()>0) {
    for (int ii=0; ii<(int)(vmrp.size()); ii++) {
      if (vmrp[ii] && ((myObject*)vmrp[ii])->iMySelf()==index) {
        return vmrp[ii];
      }
    }
  }
  
  return NULL;
}

myRichRingPlus* myEvent::pmyRichRingPlus(int index){
  
  if (index>=0 && (int)(vmrp.size())>index)
    return vmrp[index];
  
  return NULL;
}

bool pnpoly(int npol, float *xp, float *yp, float x, float y){
  
  int i, j;
  bool c = false;
  
  for (i = 0, j = npol-1; i < npol; j = i++) {
    if (
        (((yp[i] <= y) && (y < yp[j])) || 
         ((yp[j] <= y) && (y < yp[i]))) && 
        (x < (xp[j] - xp[i]) * (y - yp[i])/(yp[j] - yp[i]) +  xp[i]))
      c = !c;
  }
  
  return c;
}

Bool_t myEvent::InSAA() {
  
  static float Phi_SAA[5] = {-83.0, -38.0, 0.0, 0.0, -83.0};
  static float Theta_SAA[5] = {-1.0, -1.0, -22.0, -60.0, -60.0};
  
  return pnpoly(5, Phi_SAA, Theta_SAA, Lon(), Lat());
}

Int_t myEvent::isInShadow(myPoint &pt, int ipart){

  myParticle* mp = pmyParticle(ipart);

  if (mp) {
    return mp->isInShadow(pt);
  }

  return -1;//no particle
}

#ifdef _WITHGBATCH_
int myEvent::GetRTI(AMSSetupR::RTI & a){
  
  return AMSEventR::getsetup()->getRTI(a,this->fHeader.Time[0]);
}

#endif

int myEvent::LoadRTI(unsigned int t1, unsigned int t2, const char *dir){

  if(!myChain::gethead()->HasSectionActivated("RTI") ){
    cerr<< " Section RTI not available or not activated "<<endl;
    return 2;
  }
  TTree* t3rti = myChain::gethead()->GetSection("RTI")->t3friend->GetTree();
  if (!t3rti){
    cerr<< " RTI TTree friend not available "<<endl;
    return 2;
  }
  if(t1>t2){
    cerr<< "LoadmyRTI-S-BegintimeNotLessThanEndTime "<<t1<<" "<<t2<<endl;
    return 2;
  }
  else if(t2-t1>864000){
    cerr<< "LoadmyRTI-S-EndBeginDifferenceTooBigMax864000 "<<t2-t1<<endl;
    if(t1<1305000000){
      t1=t2-864000;
    }
    else{
      t2=t1+864000;
    }
  }

  myRTI *a;
  t3rti->SetBranchAddress("fmyRTI",a);
  int bfound=0;
  int efound=0;
  for(int i_rti=0; i_rti<t3rti->GetEntries(); i_rti++){
    t3rti->GetEntry(i_rti);
    if(!a) continue;
    if(a->run==0) continue; //missing second
    unsigned int nt=a->utime;
    if(nt<t1){bfound=1;continue;}//beg found 
    if(nt>=t1&&nt<=t2){
      if(bfound!=2){
        fRTI.clear();bfound=2;
      }
      fRTI.insert(make_pair(nt,*a));
    }//found                             
    else if(nt>t2){efound=1; break;}//end  
  }
  int ret;
  if(bfound&&efound)ret=0;//All found
  else if(!bfound&&!efound)ret=2;//All not found
  else ret=1;//Only One Found
  cout<< "LoadmyRTI-I- "<<fRTI.size()<<" Entries Loaded "<<ret<<endl;

  return ret;
}                                               


int myEvent::GetRTI(myRTI & a){
  

  a = *(this->RTI);
  if(this->RTI) return 0;  
  else  return 1;
}


int myEvent::GetRTI(myRTI & a, unsigned int  xtime){

  static unsigned int stime[2]={1,1};
  static int vrti=-1;
  
  const int pt=3600*3;
  const int dt=3600*24;
  if(stime[0]==1||xtime>stime[1]||xtime<stime[0]){
    stime[0]=(xtime<=pt)?1:xtime-pt;
    stime[1]=xtime+dt;
    fRTI.clear();
    this->LoadRTI(stime[0],stime[1]);
  }
  //---Status                                                  
  myRTI b;
  a=b;
  map< unsigned int, myRTI >::iterator k=fRTI.lower_bound(xtime);
  if (fRTI.size()==0)return 2;
  if(k==fRTI.end())return 1;

  if(xtime==k->first){//find
    a=(k->second);
    return 0;
  }
  return 2;
}

int myEvent::GetRTIdL1L9(int extlay,myPoint &nxyz, myPoint &dxyz,unsigned int xtime,unsigned int dt, int rti_flag){
  unsigned int bt= xtime/dt*dt;
  unsigned int et= bt+dt-1;

  static unsigned int rbt=0;
  static unsigned int ret=0;
  static unsigned int rst=0;//Eff  dt
  static myPoint rnxyz[2];
  static myPoint rdxyz[2];

  //----New External Layer If
  if(!(bt>=rbt&&et<=ret)){ //First Time arrive
    rst=0;
    for(int iexl=0;iexl<2;iexl++){rnxyz[iexl].setp(0,0,0);rdxyz[iexl].setp(0,0,0);}
    //----Sum BT
#ifdef _WITHGBATCH_
    AMSSetupR::RTI a;
#endif
    myRTI mya;
    for(unsigned int t=bt;t<=et;t++){
#ifdef _WITHGBATCH_
      if(rti_flag==0){
        if(GetRTI(mya,t)!=0) continue;
      }else{
        if(AMSEventR::GetRTI(a,t)!=0)continue;
      }
#else
      if(GetRTI(mya,t)!=0) continue;
#endif
      for(int iexl=0;iexl<2;iexl++){
        for(int ico=0;ico<3;ico++){
          double nev=0;
#ifdef _WITHGBATCH_
          if(rti_flag==0){
            nev=a.nl1l9[iexl][ico==0?0:1];
            rdxyz[iexl][ico]+=a.dl1l9[iexl][ico]*nev;
          }
          else{
            nev=mya.nl1l9[iexl][ico==0?0:1];
            rdxyz[iexl][ico]+=mya.dl1l9[iexl][ico]*nev;
          }
#else     
          nev=mya.nl1l9[iexl][ico==0?0:1];
          rdxyz[iexl][ico]+=mya.dl1l9[iexl][ico]*nev;
#endif
          rnxyz[iexl][ico]+=nev;//X Y Z Events number
        }
      }
      rst++;
    }
    //---Calculate Mean
    for(int iexl=0;iexl<2;iexl++){
      for(int ico=0;ico<3;ico++){
        rdxyz[iexl][ico]=(rnxyz[iexl][ico]>0)?rdxyz[iexl][ico]/rnxyz[iexl][ico]:0;
      }
    }
  }

  nxyz=rnxyz[extlay]; dxyz=rdxyz[extlay];
  rbt=bt;     ret=et;
  return rst;
}


//---------------------------------------------

//-------------------------------------------------------------------------        

#ifdef _WITHGBATCH_

//--------------------------------------------------------------------

ClassImp(myEventFiller);

//--------------------------------------------------------------------

static int TRDWeakForEPSeparation(AMSEventR* pev);
static int TrdRunCategory(AMSEventR* pev);
static bool BadRunCheck(AMSEventR *pev);
static int IsScienceRun(AMSEventR *pev);
static int RunAnalysisCategory(AMSEventR *pev);

/*
bool myEventFiller:: NoTRDPlus=false;
bool myEventFiller:: NoTRDKPlus=false;
bool myEventFiller:: NoTRDQtPlus=false;
bool myEventFiller:: NoEcalPlus=false;
bool myEventFiller:: NoBetaPlus=false;
bool myEventFiller:: NoBetaHPlus=false;
bool myEventFiller:: NoRichPlus=false;
bool myEventFiller:: NoTrackerPlus=false;
bool myEventFiller:: NoHeaderTree=false;
*/

myEventFiller::myEventFiller():myEvent(){
#ifdef PDEBUG
  printf("In myEventFiller::myEventFiller\n");
#endif
  PRINTDEBUG;
  init();
  PRINTDEBUG;
  ptr=this;
}

myEventFiller::~myEventFiller(){
#ifdef PDEBUG
  printf("In myEventFiller::~myEventFiller\n");
#endif
  PRINTDEBUG;
  Clear();
  PRINTDEBUG;
  ptr = 0;
  PRINTDEBUG;
}

myEventFiller* myEventFiller::gethead() {
#ifdef PDEBUG
  printf("In myEventFiller::gethead\n");
#endif

  PRINTDEBUG;

  if (!ptr) {
    ptr = new myEventFiller();
    printf("Object myEventFiller (%p) created as singleton...\n", ptr);
  }
  else printf("myEventFiller (%s) exists!\n", ptr->ClassName());//only for debug
  return (myEventFiller*)ptr;
}
 
void myEventFiller::init() {
#ifdef PDEBUG
  printf("In myEventFiller::init\n");
#endif

  myEvent::init();
  
  for (int ii=0; ii<9; ii++) {
    if(TkDBc::Head) zlay[ii]=TkDBc::Head->GetZlayer(ii+1);
    //    printf("%f\n", zlay[ii]);
  }
  
  Clear();

  return;
}

void myEventFiller::Clear(Option_t*){
#ifdef PDEBUG
  printf("In myEventFiller::Clear\n");
#endif
  PRINTDEBUG;

  myEvent::Clear();

  PRINTDEBUG;

  tmp_filler.SetOwner(kTRUE);
  tmp_filler.Clear();//since 'owner' this destroy all the myObject's hold by it

  PRINTDEBUG;

  mcevent = 0;

  pev = 0;

  head = 0;

  lvl1 = 0;

  tkrechitsevent = 0;

  setup = 0;

  slowcontrol = 0;

  daqev = 0;

#if !defined _B524_ && !defined _B550_ && !defined _B572_
  pdsperr = 0;
#endif

  PRINTDEBUG;
}

void myEventFiller::Fill(int fTree, Long64_t fEntry){
#ifdef PDEBUG
  printf("In myEventFiller::Fill\n");
#endif

  PRINTDEBUG;

  Clear();

  PRINTDEBUG;
  /* TODO not optimized, string finding every event... */

  myHeaderFiller mhfill;  
  mhfill.Fill();   
  fHeader = (myHeader) mhfill;

  PRINTDEBUG;

  Fill();

  PRINTDEBUG;

  mcevent = pev->pMCEventg(0);
  if (mcevent) {
    myMCFiller mmfill;
    mmfill.Fill(mcevent);
    mm = new myMC();
    *mm = mmfill;
    kMC = true;
    //    printf("MC block...\n");//only for debug
  }
  else kMC = false;

  PRINTDEBUG;

#ifdef PDEBUG
  printf("Loop over TrdTrack's (%d)\n", _nTrdTrack);//only for debug
#endif
  for (int ii=0; ii<_nTrdTrack; ii++) {
    TrdTrackR* trdtrack = pev->pTrdTrack(ii);
    PRINTDEBUG;
    if (trdtrack) {
      myTrdTrackFiller mu;
      mu.SetMySelf(ii);
      PRINTDEBUG;
      mu.Fill(trdtrack);
      PRINTDEBUG;
      myTrdTrack* mu2bewritten = new myTrdTrack(mu);
      tmp_filler.Add(mu2bewritten);
      vmu.push_back(mu2bewritten);
      vmukup.push_back(nullptr); // it's indicized with trd track
      // so far it has NO plus
    }
    PRINTDEBUG;
  }

  PRINTDEBUG;
  if (gsection.enabled("TRDPlus")) {
    if(pev->nTrdSegment()>0){
      myTrdTrackPlusFiller mup;
      mup.SetMySelf(0);
      PRINTDEBUG;
      mup.Fill(pev);
      PRINTDEBUG;
      myTrdTrackPlus* mup2bewritten = new myTrdTrackPlus(mup);
      tmp_filler.Add(mup2bewritten);
      vmup.push_back(mup2bewritten);
      PRINTDEBUG;
    }
  }
  
  PRINTDEBUG;

#ifdef PDEBUG
  printf("Loop over TrTrack's (%d)\n", _nTrTrack);
#endif
  for (int ii=0; ii<_nTrTrack; ii++) {
    TrTrackR* track = pev->pTrTrack(ii);
    PRINTDEBUG;
    if (track) {
      myTrTrackFiller mt;
      mt.SetMySelf(ii);
      PRINTDEBUG;
      mt.Fill(track);
      PRINTDEBUG;
      myTrTrack* mt2bewritten = new myTrTrack(mt);
      tmp_filler.Add(mt2bewritten);
      vmt.push_back(mt2bewritten);
      vmtp.push_back(NULL);
      vmuktp.push_back(NULL);//also the myTrdKFromTrTrack is 'indicized' as the myTrTrack
      vmuqtp.push_back(NULL);//also the myTrdQtFromTrTrack is 'indicized' as the myTrTrack
#ifdef PDEBUG
      printf("track=%d) iPart=%d, iMySelf=%d, i_part=%d, i_myself=%d\n", ii, mt.iPart(), mt.iMySelf(), mt.i_part, mt.i_myself);
#endif
    }
    PRINTDEBUG;
  }

  PRINTDEBUG;

#ifdef PDEBUG
  printf("Loop over EcalShower's (%d)\n", _nEcalShower);
#endif
  for (int ii=0; ii<_nEcalShower; ii++) {
    EcalShowerR* shower = pev->pEcalShower(ii);
    PRINTDEBUG;
    if (shower) {
      myEcalShowerFiller me;
      me.SetMySelf(ii);
      PRINTDEBUG;
      me.Fill(shower);
      PRINTDEBUG;
      myEcalShower* me2bewritten = new myEcalShower(me);
      tmp_filler.Add(me2bewritten);
      vme.push_back(me2bewritten);
      vmep.push_back(NULL);
    }
    PRINTDEBUG;
  }

  PRINTDEBUG;

#ifdef PDEBUG
  printf("Loop over RichRing's (%d)\n", _nRichRing);//only for debug
#endif
  RichRingR::setBetaCorrection(RichRingR::fullUniformityCorrection);
  for (int ii=0; ii<_nRichRing; ii++) {
    RichRingR* ring = pev->pRichRing(ii);
    PRINTDEBUG;
    if (ring) {
      myRichRingFiller mr;
      mr.SetMySelf(ii);
      PRINTDEBUG;
      mr.Fill(ring);
      PRINTDEBUG;
      myRichRing* mr2bewritten = new myRichRing(mr);
      tmp_filler.Add(mr2bewritten);
      vmr.push_back(mr2bewritten);
      vmrp.push_back(NULL);
    }
    PRINTDEBUG;
  }

  PRINTDEBUG;

#ifdef PDEBUG
  printf("Loop over Beta's (%d)\n", _nBeta);
#endif
  for (int ii=0; ii<_nBeta; ii++) {
    BetaR* beta = pev->pBeta(ii);
    PRINTDEBUG;
    if (beta) {
      myBetaFiller mb;
      mb.SetMySelf(ii);
      PRINTDEBUG;
      mb.Fill(beta);
      PRINTDEBUG;
      myBeta* mb2bewritten = new myBeta(mb);
      tmp_filler.Add(mb2bewritten);
      vmb.push_back(mb2bewritten);//questo crasha
      vmbp.push_back(NULL);
    }
    PRINTDEBUG;
  }

  PRINTDEBUG;

#ifdef PDEBUG
  printf("Loop over BetaH's (%d)\n", _nBetaH);
#endif
  for (int ii=0; ii<_nBetaH; ii++) {
    BetaHR* betah = pev->pBetaH(ii);
    PRINTDEBUG;
    if (betah) {
      myBetaHFiller mbh;
      mbh.SetMySelf(ii);
      PRINTDEBUG;
      mbh.Fill(betah);
      PRINTDEBUG;
      myBetaH* mbh2bewritten = new myBetaH(mbh);
      tmp_filler.Add(mbh2bewritten);
      vmbh.push_back(mbh2bewritten);
      vmbhp.push_back(NULL);
    }
    PRINTDEBUG;
  }

  PRINTDEBUG;

  int i_part_max_qbh=-1;
  float max_qbh=0.;

  int i_part_max_r=-1;
  float max_r=0.;
  int i_part_max_r_trd=-1;
  float max_r_trd=0.;
  int i_part_max_e=-1;
  float max_e=0.;
  //  int i_part_complete=-1;

  PRINTDEBUG;

#ifdef PDEBUG
  printf("Loop over Particle's (%d)\n", _nParticle);
#endif
  
  for (int ii=0; ii<_nParticle; ii++) {
    myParticleFiller mp;
    mp.Fill(pev->pParticle(ii));
    mp.SetPart(ii);
    mp.SetMySelf(ii);
    
    PRINTDEBUG;

    if (mp.trdtrack && mp.pmyTrdTrackByMySelf()) mp.pmyTrdTrackByMySelf()->SetPart(ii);
    PRINTDEBUG;
    if (mp.track && mp.pmyTrTrackByMySelf()) mp.pmyTrTrackByMySelf()->SetPart(ii);
    PRINTDEBUG;
    if (mp.shower && mp.pmyEcalShowerByMySelf()) mp.pmyEcalShowerByMySelf()->SetPart(ii);
    PRINTDEBUG;
    if (mp.ring && mp.pmyRichRingByMySelf()) mp.pmyRichRingByMySelf()->SetPart(ii);
    PRINTDEBUG;
    if (mp.beta && mp.pmyBetaByMySelf()) mp.pmyBetaByMySelf()->SetPart(ii);
    PRINTDEBUG;
    if (mp.betah && mp.pmyBetaHByMySelf()) mp.pmyBetaHByMySelf()->SetPart(ii);
    PRINTDEBUG;
#ifdef PDEBUG
    if (mp.track) printf("part=%d) mp.pmyTrTrack(), iPart=%d, iMySelf=%d, i_part=%d, i_myself=%d\n", ii, mp.pmyTrTrackByMySelf()->iPart(), mp.pmyTrTrackByMySelf()->iMySelf(), mp.pmyTrTrackByMySelf()->i_part, mp.pmyTrTrackByMySelf()->i_myself);
    if (mp.track) printf("mpart=%d) vmt[%d], iPart=%d, iMySelf=%d, i_part=%d, i_myself=%d\n", ii, mp.iTrTrack(), vmt[mp.iTrTrack()]->iPart(), vmt[mp.iTrTrack()]->iMySelf(), vmt[mp.iTrTrack()]->i_part, vmt[mp.iTrTrack()]->i_myself);
#endif
    
    PRINTDEBUG;

    bool NoEcalPlus = not gsection.enabled("EcalPlus"),
         NoHeaderTree = not gsection.enabled("HeaderTree");
    if (!(NoEcalPlus && NoHeaderTree) ) {
      if (mp.shower) {
        myEcalShowerPlusFiller* mep = new myEcalShowerPlusFiller();
        PRINTDEBUG;
        // cases:
        // NoEcalPlus | NoHeaderTree  --> Fill(EcalShowerR*, short)
        //     0             0        --> Fill(EcalShowerR*, 0), Fill is called normally, OK
        //     0             1        --> Fill(EcalShowerR*, 1), Fill is called normally, OK
        //     1             0        --> Fill(EcalShowerR*, 1), Fill is called "short" , OK
        //     1             1        --> Fill(EcalShowerR*, 1), Fill is called normally, KO
        //                                but we didn't arrive here since there's the if on top
        mep->Fill(mp.shower, NoEcalPlus&(!NoHeaderTree));
        PRINTDEBUG;
        mep->SetPart(ii);
        mep->SetMySelf(mp.pp->iEcalShower());
        myEcalShowerPlus* mep2bewritten = new myEcalShowerPlus(*mep);
        tmp_filler.Add(mep2bewritten);
        vmep[mep2bewritten->iMySelf()] = mep2bewritten;
        if(mp.shower->EnergyD>=max_e){
          max_e=mp.shower->EnergyD;
          i_part_max_e=ii;
        }
        if (mep) delete mep;
        mep=0;
      }
    }

    PRINTDEBUG;

    /* have to put here because Trd*plus need track */
    bool NoTrackerPlus = not gsection.enabled("TrackerPlus");
      if (!NoTrackerPlus) {
        if (mp.track) {
          myTrTrackPlusFiller* mtp = new myTrTrackPlusFiller();
          PRINTDEBUG;
          Double_t beta=999;
          if (mp.pmyBetaHByMySelf()) beta=mp.pmyBetaHByMySelf()->Beta;
          mtp->Fill(mp.track, ii, mp.shower, beta, kMC, true);
          PRINTDEBUG;
          mtp->FillCharge(mp.track, beta);
          PRINTDEBUG;
          mtp->SetPart(ii);
          mtp->SetMySelf(mp.pp->iTrTrack());
          myTrTrackPlus* mtp2bewritten = new myTrTrackPlus(*mtp);
          tmp_filler.Add(mtp2bewritten);
          vmtp[mtp2bewritten->iMySelf()] = mtp2bewritten;
          if(mtp->GetRigidity(0)>=max_r){
            max_r=mtp->GetRigidity(0);
            i_part_max_r=ii;
          }
          if (mtp) delete mtp;
          mtp=0;
        }
      }

      PRINTDEBUG;

      if( gsection.enabled("BetaHPlus") ){
        if ( mp.betah && mp.pmyBetaHByMySelf()) {
      myBetaHPlusFiller* mbhp = new myBetaHPlusFiller();
      PRINTDEBUG;
      mbhp->Fill(mp.betah);
      PRINTDEBUG;
      mbhp->SetPart(ii);
      mbhp->SetMySelf(mp.pp->iBetaH());
      myBetaHPlus* mbhp2bewritten = new myBetaHPlus(*mbhp);
      tmp_filler.Add(mbhp2bewritten);
      vmbhp[mbhp2bewritten->iMySelf()] = mbhp2bewritten;
      if(mbhp->Q>=max_qbh){
        max_qbh=mbhp->Q;
        i_part_max_qbh=ii;
      }
      if (mbhp) delete mbhp;
      mbhp=0;
    }
      }

    PRINTDEBUG;

    if( gsection.enabled("BetaPlus")){
    if (mp.beta && mp.pmyBetaByMySelf()) {
      myBetaPlusFiller* mbp = new myBetaPlusFiller();
      PRINTDEBUG;
      if (mp.pmyTrTrackByMySelf()) mbp->Fill(mp.beta, mp.track, mp.pmyTrTrackByMySelf()->kDef);
      else mbp->Fill(mp.beta, mp.track);
      PRINTDEBUG;
      mbp->SetPart(ii);
      mbp->SetMySelf(mp.pp->iBeta());
      myBetaPlus* mbp2bewritten = new myBetaPlus(*mbp);
      tmp_filler.Add(mbp2bewritten);
      vmbp[mbp2bewritten->iMySelf()] = mbp2bewritten;
      if (mbp) delete mbp;
      mbp=0;
    }
    }

    PRINTDEBUG;

    if( gsection.enabled("RichPlus")){
    if (mp.ring && mp.pmyRichRingByMySelf()) {
      myRichRingPlusFiller* mrp = new myRichRingPlusFiller();
      PRINTDEBUG;
      mrp->Fill(mp.ring, mp.pp);
      PRINTDEBUG;
      mrp->SetPart(ii);
      mrp->SetMySelf(mp.pp->iRichRing());
      myRichRingPlus* mrp2bewritten = new myRichRingPlus(*mrp);
      tmp_filler.Add(mrp2bewritten);
      vmrp[mrp2bewritten->iMySelf()] = mrp2bewritten;
      if (mrp) delete mrp;
      mrp=0;
    }
    }

    PRINTDEBUG;
    if (gsection.enabled("TRDKPlus")) {
      if (mp.track) {
        myTrdKFromTrTrackFiller* mukp = new myTrdKFromTrTrackFiller();
        PRINTDEBUG;
        float energy=0;
        if (mp.pmyEcalShowerByMySelf()) energy=mp.pmyEcalShowerByMySelf()->EnergyE;
        int kDef=0;
        if (mp.pmyTrTrackPlusByMySelf()) kDef=mp.pmyTrTrackPlusByMySelf()->kDef;
        if( kDef==0 ){
           if(!NoTrackerPlus)Error(__func__,"strange kDef for %p", mp.pmyTrTrackPlusByMySelf());
           kDef = mp.track->Gettrdefaultfit();
        }
        mukp->Fill(mp.iTrTrack(), mp.track, energy, kDef);
        PRINTDEBUG;
        mukp->SetPart(ii);
        mukp->SetMySelf(mp.pp->iTrTrack());//TrdKFromTrTrack is 'indicized' with the TrTrack used
        myTrdK* mukp2bewritten = new myTrdK(*mukp);
        tmp_filler.Add(mukp2bewritten);
        vmuktp[mukp2bewritten->iMySelf()] = mukp2bewritten;
        if(mp.pmyTrTrackByMySelf() && mp.pmyTrTrackByMySelf()->Rigidity>=max_r_trd){
        max_r_trd=mp.pmyTrTrackByMySelf()->Rigidity;
        i_part_max_r_trd=ii;
        }
        if (mukp) delete mukp;
        mukp=0;
      }
      if (mp.trdtrack ) {
        myTrdKFromTrdTrackFiller mukup;
        PRINTDEBUG;
        float rig = std::abs(pev->Particle(ii).Momentum / pev->Particle(ii).Charge);
        float energy=0;
        if (mp.pmyEcalShowerByMySelf()) energy=mp.pmyEcalShowerByMySelf()->EnergyE;
        mukup.Fill(mp.iTrdTrack(), mp.trdtrack, energy, rig);
        PRINTDEBUG;
        mukup.SetPart(ii);
        mukup.SetMySelf(mp.pp->iTrdTrack());//TrdKFromTrdTrack is 'indicized' with the TrdTrack used
        myTrdK* mukp2bewritten = new myTrdK(mukup);
        tmp_filler.Add(mukp2bewritten);
        vmukup[mukp2bewritten->iMySelf()] = mukp2bewritten;
      }
    }

    if (gsection.enabled("TRDQtPlus")) {
      if (mp.track) {
        myTrdQtFromTrTrackFiller* muqtp = new myTrdQtFromTrTrackFiller();
        PRINTDEBUG;
        float energy=0;
        if (mp.pmyEcalShowerByMySelf()) energy=mp.pmyEcalShowerByMySelf()->EnergyE;
        int kDef=0;
        if (mp.pmyTrTrackByMySelf()) kDef=mp.pmyTrTrackByMySelf()->kDef;
        muqtp->Fill(mp.iTrTrack(), mp.track, kDef, energy);
        PRINTDEBUG;
        muqtp->SetPart(ii);
        muqtp->SetMySelf(mp.pp->iTrTrack());//TrdQtFromTrTrack is 'indicized' with the TrTrack used
        myTrdQtFromTrTrack* muqtp2bewritten = new myTrdQtFromTrTrack(*muqtp);
        tmp_filler.Add(muqtp2bewritten);
        vmuqtp[muqtp2bewritten->iMySelf()] = muqtp2bewritten;
        if (muqtp) delete muqtp;
        muqtp=0;
      }
    }

    PRINTDEBUG;

    myParticle* mp2bewritten = new myParticle(mp);
    tmp_filler.Add(mp2bewritten);
    vmp.push_back(mp2bewritten);
    PRINTDEBUG;
  }

  PRINTDEBUG;
  if( gsection.enabled("StatusTree")){
#ifdef PDEBUG
  printf("Filling Status \n");//only for debug
#endif
  PRINTDEBUG;
  myStatusFiller msfill;
  msfill.fStatus=pev->fStatus;
  msfill.EntryNumber=fEntry;
  msfill.TreeNumber=fTree;
  msfill.Fill(i_part_max_r, i_part_max_qbh , i_part_max_e, i_part_max_r_trd);
  PRINTDEBUG;
  //FIXME: why not copy constructor?!
  ms = new myStatus();
  *ms = msfill;
  }
  PRINTDEBUG;

  return;
}

void myEventFiller::Fill(){
#ifdef PDEBUG
  printf("In myEventFiller::Fill\n");
#endif
  PRINTDEBUG;
  
  pev = AMSEventR::Head();
  
  string reason;
  if (AMSEventR::getsetup())
    AMSEventR::getsetup()->IsBadRun(reason, pev->UTime(), pev->Run());
  ReasonToBeBadRun = reason;
  _IsBadRun = pev->IsBadRun("");
  RunAnalCat = RunAnalysisCategory(pev);

  Run = pev->Run();
  Event = pev->Event();
  
  head = &(pev->fHeader);


  _UTCTime = pev->UTCTime() ;

  PRINTDEBUG;
  
  if (!kMC) {
    static unsigned int TimeJMDCPrevious = fHeader.Time[0];//in this way also the first round the deltat is correct (0)
    static unsigned int MuTimeJMDCPrevious = fHeader.Time[1];//in this way also the first round the deltat is correct (0)
    DeltaTJMDC = fHeader.Time[0]+0.000001*fHeader.Time[1] - (TimeJMDCPrevious+0.000001*MuTimeJMDCPrevious);
    TimeJMDCPrevious = fHeader.Time[0];
    MuTimeJMDCPrevious = fHeader.Time[1];
    
    if (fHeader.Time[0]<TimeJMDCFirst) TimeJMDCFirst = fHeader.Time[0];
    if (fHeader.Time[0]>TimeJMDCLast) TimeJMDCLast = fHeader.Time[0];
    Duration = TimeJMDCLast - TimeJMDCFirst;//the total time in the "box"
  }
  
  PRINTDEBUG;
  
  lvl1 = pev->pLevel1(0);
  
  setup = pev->getsetup();
  
  daqev=pev->pDaqEvent(0);
  
  //-------------------------------------------------------------------------------------
  if (!kMC) {
    //  printf("(Lat, Lon, Alt) = (%f, %f, %f)\n", Lat(), Lon(), Alt()); //only for debug
    
    GM_GeoMagSphericalCoo(fHeader.Time[0], fHeader.RadS, fHeader.ThetaS, fHeader.PhiS, GeoMagSphericalCoo);

    //  printf("(LatMAG, LonMAG, RadMAG) = (%f, %f, %f)\n", LatMAG(), LonMAG(), RadMAG()); //only for debug 
    
#ifdef PREDICT
    static ISSCoo iss;
#endif
#ifdef CGM
    static cgmcoo cgm;
#endif
#ifdef FIELD
    static GeoMag geomag;
#endif
        
    if (InSAA()) {  
      static int counter=0;
      if (counter<100) {
        counter++;
        printf("%d-%d) in SAA (Lon,Lat)=(%f,%f)\n", Run, Event, Lon(), Lat());
      }
      else if (counter==100) {
        counter++;
        printf("%d-%d) in SAA (Lon,Lat)=(%f,%f). I said %d times, it's enough...\n", Run, Event, Lon(), Lat(), counter);
      }
    }
    
    if ( isISS() and fHeader.Time[0] != TimePrevISS) {//once per second it's enough...
      TimePrevISS = fHeader.Time[0];
      
#ifdef PREDICT
      iss.TrackISS((int)(fHeader.Time[0]), 0);
      //      iss.Print();
#endif
#ifdef CGM
      cgm.GetCGM(Lat, Lon, Alt, (int)(1970+(fHeader.Time[0]/(365*24*60*60))));
      //      cgm.GetCGM(iss.Lat, iss.Lon, iss.Alt);
      //      cgm.Print();
#endif
#ifdef FIELD
      geomag.Recalc(Lat, Lon, Alt, (int)(1970+(fHeader.Time[0]/(365*24*60*60))));
      //      geomag.Recalc(iss.Lat, iss.Lon, iss.Alt, 1970+(fHeader.Time[0]/(365*24*60*60)));
#endif
      
      Double_t cutoff[2];
      pev->GetMaxGeoCutoff(35.0, 1, cutoff);
      MaxGeoCutOff_n=cutoff[0];
      MaxGeoCutOff_p=cutoff[1];      

      if(gsection.enabled("RTI")){
      AMSSetupR::RTI::Version=2;
      AMSSetupR::RTI ev_RTI;
      int RTI_ko=pev->GetRTI(ev_RTI);
      if (!RTI_ko) {
        myRTIFiller rtifill;
        rtifill.Fill(&ev_RTI);
        if (RTI) delete RTI;
        RTI = new myRTI();
        *RTI = rtifill;
        //      printf("%d -- %d\n", Event, RTI->evno);//only for debug
      }else Error(__func__, "fail getting RTI for %d: %d", fHeader.Time[0], RTI_ko);
      }
    } // once per second
        
#ifdef PREDICT
    LonPred=(Double_t)((iss.Lon<180)?iss.Lon:(iss.Lon-360));
    LatPred=(Double_t)(iss.Lat);
    AltPred=(Double_t)(iss.Alt);
    //      printf("(LatPred, LonPred, AltPred) = (%f, %f, %f)\n", LatPred, LonPred, AltPred); //only for debug
#endif
#ifdef CGM
    //      LonCGM=(Double_t)((cgm.clor<180)?cgm.clor:(cgm.clor-360));
    LonCGM=M_PI*(Double_t)(cgm.clor)/180.0;
    LatCGM=M_PI*(Double_t)(cgm.clar)/180.0;
    RadCGM=((Double_t)(cgm.rbm));
    //      printf("(LatCGM, LonCGM, AltCGM) = (%f, %f, %f)\n", LatCGM, LonCGM, RadCGM); //only for debug
    //      printf("CGM radius (~ L) = %f -> %f km\n", cgm.rbm, EARTHRADIUS*RadCGM); //only for debug
#endif
#ifdef FIELD
    Field=geomag.F;
#endif
        
    if (setup) {
      if (fHeader.Time[0] != TimePrevSlow) {
        TimePrevSlow = fHeader.Time[0];
        slowcontrol = &(setup->fSlowControl);
#ifdef SLOWCONTROLTEMP
        std::vector<float> temp;
        slowcontrol->GetData("Sensor A", fHeader.Time[0], (float)(fHeader.Time[1]/1000000.0), temp);
        if (temp.size()>0) SensorA = temp[0];
        slowcontrol->GetData("Sensor C", fHeader.Time[0], (float)(fHeader.Time[1]/1000000.0), temp);
        if (temp.size()>0) SensorC = temp[0];
        slowcontrol->GetData("Sensor K", fHeader.Time[0], (float)(fHeader.Time[1]/1000000.0), temp);
        if (temp.size()>0) SensorK = temp[0];
        slowcontrol->GetData("Sensor L", fHeader.Time[0], (float)(fHeader.Time[1]/1000000.0), temp);
        if (temp.size()>0) SensorL = temp[0];
        slowcontrol->GetData("Sensor M", fHeader.Time[0], (float)(fHeader.Time[1]/1000000.0), temp);
        if (temp.size()>0) SensorM = temp[0];
        slowcontrol->GetData("M-7X 0 1N-1", fHeader.Time[0], (float)(fHeader.Time[1]/1000000.0), temp);
        if (temp.size()>0) Plane1N_M70 = temp[0];
        slowcontrol->GetData("M-7X 0 1N-2", fHeader.Time[0], (float)(fHeader.Time[1]/1000000.0), temp);
        if (temp.size()>0) Plane1N_M71 = temp[0];
        slowcontrol->GetData("M-7X 0 1N-3", fHeader.Time[0], (float)(fHeader.Time[1]/1000000.0), temp);
        if (temp.size()>0) Plane1N_M72 = temp[0];
        slowcontrol->GetData("M-7X 0 1N-4", fHeader.Time[0], (float)(fHeader.Time[1]/1000000.0), temp);
        if (temp.size()>0) Plane1N_M73 = temp[0];
        slowcontrol->GetData("Sensor 8", fHeader.Time[0], (float)(fHeader.Time[1]/1000000.0), temp);
        if (temp.size()>0) Sensor8 = temp[0];
#endif
      }
    }
    
    if (daqev) {
      Lenght=daqev->Length;
      Tdr=daqev->Tdr;
      Udr=daqev->Udr;
      Sdr=daqev->Sdr;
      Rdr=daqev->Rdr;
      Edr=daqev->Edr;
      L1dr=daqev->L1dr;
      L3dr=daqev->L3dr;
      L3Error=daqev->L3Error;
      L3VEvent=daqev->L3VEvent;
      L3TimeD=daqev->L3TimeD;
#if !defined _B524_ && !defined _B550_
      copy_n( daqev->JINJStatus, 4, JINJStatus);
      copy_n( daqev->JError,    24, JError);
#endif
    }
    
#if defined(DSPERRORS) && !defined _B524_ && !defined _B550_ && !defined _B572_
    pdsperr=head->pDSPError();
#endif
    
  }
  
  //------------------------------------------------------------------------------------------------------------------------
  
  PRINTDEBUG;
  
  //--------------------------LVL1 & TOF-----------------------------
  if (lvl1) {
    
    EcalFlag = lvl1->EcalFlag;
    copy_n( lvl1->TrigRates, 19,TrigRates);
    copy_n( lvl1->TofPatt1, 4  ,TofPatt1);
    copy_n( lvl1->TofPatt2, 4  ,TofPatt2);
    JMembPatt = lvl1->JMembPatt;
    PhysBPatt=lvl1->PhysBPatt;
    PhysBPattPhys=PhysBPatt&0x3E;//(physbpatt&0x3E) is the "physic" part of PhysBPatt (from bit 1 to bit 5)
    lvl1->RestorePhysBPat();
    PhysBPattRestored=lvl1->PhysBPatt;
    TofFlag1 = lvl1->TofFlag1;
    TofFlag2 = lvl1->TofFlag2;
    LiveTime = lvl1->LiveTime;
    DeltaT = lvl1->TrigTime[4]; // in musec
    Lvl1Time = 0.64*(lvl1->TrigTime[2] + pow(2.0, 32.0)*lvl1->TrigTime[3]);
    DeltaTMine = Lvl1Time - Lvl1TimePrevious;
    if (DeltaTMine<0) DeltaTMine = Lvl1Time - (Lvl1TimePrevious - lvl1->TrigTime[0]*0.64);
    //    printf("%d) DeltaT = %u - DeltaTMine = %f - Diff = %E\n", Event, DeltaT, DeltaTMine, DeltaT-DeltaTMine);
    //    printf("%u %u %u %u %u\n", lvl1->TrigTime[0], lvl1->TrigTime[1], lvl1->TrigTime[2], lvl1->TrigTime[3], lvl1->TrigTime[4]);
    //    printf("JMDCtime = %d,%06d, LVL1-time = %f\n", fHeader.Time[0], fHeader.Time[1], Lvl1Time);
    Lvl1TimePrevious = Lvl1Time;
    
    PRINTDEBUG;

    //--------------------------ANTI------------------------------    
    AntiPatt=lvl1->AntiPatt;
    nAnti=0;
    for (int ii=0; ii<=8; ii++){
      if((AntiPatt)&(1<<ii)) {
        nAnti++;//counting the number of anti fired
      }
    }

  }

  PRINTDEBUG;

  nTofCluster=pev->nTofCluster();
  for (int ii=0; ii<nTofCluster; ii++) {
    TofClusterR* tfr=pev->pTofCluster(ii);
    if (((tfr->Layer)-1)>3) printf("nTofCluster evaluation: layer number (%d) not possible!\n", tfr->Layer);
    else nTofClusterL[(tfr->Layer)-1]++;
  }

  PRINTDEBUG;
  
  //--------------------------TRD------------------------------
  _nTrdTrack = pev->nTrdTrack();
  _nTrdCluster = pev->nTrdCluster();
  _nTrdRawHit = pev->nTrdRawHit();
  _nTrdSegment = pev->nTrdSegment();

  PRINTDEBUG;

  if(pev->nMCEventg()>0){
    //---------------------TRACKER----------------------------------
    _nTrTrack = pev->NTrTrack();
    _nTrCluster = pev->NTrCluster();
    _nTrRecHit = pev->NTrRecHit();
    _nTrRawCluster = pev->NTrRawCluster();
    //--------------------------ECAL------------------------------
    _nEcalShower = pev->NEcalShower();
    
    //--------------------------TOF-------------------------------
    _nBetaH = pev->NBetaH();
    //   if(!nBetaH) TofRecH::ReBuild();   //reconstruct BetaHR object. This will also ReInit ParticleR and ChargeR AMSEventR
    _nBeta = pev->NBeta();
    
    //-------------------------RICH-------------------------------
    _nRichRing = pev->NRichRing();
    _nRichHit = pev->NRichHit();

    //-------------------------Particle---------------------------
    _nParticle = pev->NParticle();
  }
  else{
    //---------------------TRACKER----------------------------------
    _nTrTrack = pev->nTrTrack();
    _nTrCluster = pev->nTrCluster();
    _nTrRecHit = pev->nTrRecHit();
    _nTrRawCluster = pev->nTrRawCluster();
    //--------------------------ECAL------------------------------
    _nEcalShower = pev->nEcalShower();
    
    //--------------------------TOF-------------------------------
    _nBetaH = pev->nBetaH();
    //   if(!nBetaH) TofRecH::ReBuild();   //reconstruct BetaHR object. This will also ReInit ParticleR and ChargeR AMSEventR
    _nBeta = pev->nBeta();
    
    //-------------------------RICH-------------------------------
    _nRichRing = pev->nRichRing();
    _nRichHit = pev->nRichHit();

    //-------------------------Particle---------------------------
    _nParticle = pev->nParticle();
  }


  for (int ii=0; ii<_nTrRawCluster; ii++) {
    static TrClusterR* clu;
    clu = pev->pTrCluster(ii);
    if (clu) {
      static int layer;
      layer = clu->GetLayer();
      //      printf("Layer = %d\n", layer);//only for debug
      if ((layer-1)>8) printf("nTrRawCluster evaluation: layer number (%d) not possible!\n", layer);
      else _nTrRawClusterL[layer-1]++;
    }
  }
  tkrechitsevent = &(pev->TrRecHit());
  
  PRINTDEBUG;

  return;
}


// void myEventFiller::RebuildxLkTrk(myParticle& mp, myTrTrackPlus& mt, myEcalShower& me, myBetaPlus& mb){

//   TrkLHVar var1;
//   //    printf("myEventFiller::RebuildLkTrk: ipart = %d\n", ipart);//only for debug                 
//   var1.btp = mt.gbpj;   // bitpattern                              
//   var1.energy = me.EnergyE;
//   var1.rigidity_ch = mt.rigidity;  // algo 1                       
//   var1.chiy_ch = pow(10, mt.logchisqY);
//   var1.chix_ch = pow(10, mt.logchisqX);
//   var1.err_ch = mt.sigmarinv_SameWeight;   // same weigth  algo 21 
//   var1.rigidity_al = mt.rigidity_al;  // algo 2              
//   var1.rigidity_in = mt.rigidity_inn;  // pattern 3                
//   var1.rigidity_bt = mt.rigidity_k9;  // pattern 6                 
//   var1.rigidity_tp = mt.rigidity_k8;  // pattern 5                 
//   var1.rigidity_ibt = mt.rigidity_lh; // pattern 2                 
//   var1.rigidity_itp = mt.rigidity_uh; // pattern 1                 
//   var1.trkqin = mt.trackchargeinn;             // inner trk ch     
//   var1.trkqrms = mt.trackchargerms;           // rms trk char      
//   var1.minoise1 = mt.trnoisemindist[0].y();
//   var1.minoise2 = mt.trnoisemindist[1].y();
//   var1.minoise3 = mt.trnoisemindist[2].y();
//   var1.minoise4 = mt.trnoisemindist[3].y();
//   var1.nlayer = 0.0;
//   for (int layer=0; layer<=8; layer++) {
//     if((mt.gbpj)&(1<<layer))
//       var1.nlayer=var1.nlayer+1.0;  // found number hitted layer 
//   }
//   var1.tof_qd = mb.TofTrackChargePlaneLower;
//   var1.tof_qu = mb.TofTrackChargePlaneUpper;
//   var1.ntrk = 1.*(nTrTrack);                             
//   var1.nanti = 1.*(nanti);


//   // ECAL RELATED 
//   //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//   myPoint ecalp9;//=PointTrackEntryEcal;
//   myPoint ecalp0;//=PointTrackExitEcal;
//   myPoint ecalCofG=mt.PointTrackCogEcal;

//   var1.ecal_ecogx = me.CofG[0];
//   var1.ecal_ecogy = me.CofG[1];
//   var1.trk_ecogx = ecalCofG.x();
//   var1.trk_ecogy = ecalCofG.y();
//   myDir dird=mt.DirTrackExitEcal;
//   var1.cosy = fabs(dird.y());
  
//   Double_t dzn = (me.Entry[2]);
//   Double_t dzx = (me.Exit[2]);
//   if (dzn*dzn < dzx*dzx) {// standard direction                                                     
//     ecalp9.setp(me.Entry[0], me.Entry[1], me.Entry[2]);
//     ecalp0.setp(me.Exit[0], me.Exit[1], me.Exit[2]);
//   }
//   else {// contrary direction
//     ecalp0.setp(me.Entry[0], me.Entry[1], me.Entry[2]);
//     ecalp9.setp(me.Exit[0], me.Exit[1], me.Exit[2]);
//   }
//   Double_t ecp9x = ecalp9.x();
//   Double_t ecp9y = ecalp9.y();
//   Double_t ecp9z = ecalp9.z();
//   Double_t ecp0x = ecalCofG.x();
//   Double_t ecp0y = ecalCofG.y();
//   Double_t ecp0z = ecalCofG.z();
//   Double_t tecdx = dird.x();
//   Double_t tecdy = dird.y();
//   Double_t tecdz = dird.z();

//   Double_t norm1 = (ecp9x-ecp0x)*(ecp9x-ecp0x)+(ecp9y-ecp0y)*(ecp9y-ecp0y)+(ecp9z-ecp0z)*(ecp9z-ecp0z);
//   norm1 = sqrt(norm1);
//   Double_t norm2 = sqrt(tecdx*tecdx+tecdy*tecdy+tecdz*tecdz);
//   Double_t coset = (ecp9x-ecp0x)*tecdx+(ecp9y-ecp0y)*tecdy+(ecp9z-ecp0z)*tecdz;
//   coset = fabs(coset)/(norm1*norm2);
//   coset = fabs(1.-coset*coset);    // is sinsquare trk-ecal                                                                                                               
//   var1.sinq_trk_ecal = coset;
//   var1.sinq_trk_trd = 0;   // to be done and understood MC reliability of TRD track
//   //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
  
//   int nm = TrkLH::gethead()->fillvar(var1);
//   if (nm!=0) cout << " error var filling " << nm << endl;
//   //   cout << " var filled " << endl;                                                                                                                                    
//   mt.lktrk = TrkLH::gethead()->dotrklk(var1.energy, var1.pat);  // var1.pat e' il pattern                                                                                              
//   //  printf("lktrk = %f\n", lktrk);                                             
// }

#endif

#ifdef _WITHGBATCH_

int TRDWeakForEPSeparation(AMSEventR* pev){

  if (!pev) return 2;
  
  string reason = "";
  if (AMSEventR::getsetup())
    AMSEventR::getsetup()->IsBadRun(reason, pev->UTime(), pev->Run());
  else
    return 2;//cannot access to setup
  
  //  if (tag==1) cout << "Reason = " << reason << endl;//only for debug
  
  // if the run is only "weak" for EP separation (example: circulation of gas after refill) but not "bad" or "unusable" accept it
  if (reason.find("TRD_WeakForEPseparation") != string::npos && reason.find("TRD_BadForEPseparation") == string::npos && reason.find("TRD_UnusableForAnalysis") == string::npos)
    return 0;//weak for ep separation
  
  return 1;//other
}

int TrdRunCategory( AMSEventR* pev ){
  
  if (!pev) return -1;
  
  string reason = "";
  if (AMSEventR::getsetup())
    AMSEventR::getsetup()->IsBadRun(reason, pev->UTime(), pev->Run());
  else
    return -1;//cannot access to setup                                 
  
  bool weak =   ( reason.find("TRD_WeakForEPseparation") != string::npos );
  bool bad =    ( reason.find("TRD_BadForEPseparation")  != string::npos );
  bool unusab = ( reason.find("TRD_UnusableForAnalysis") != string::npos );
  
  if ( weak && !bad && !unusab )  return 1;// TRD catB
  if ( bad && !unusab )           return 2;// TRD catC
  if ( unusab )                   return 3;// TRD catD
  
  return 0;//TRD catA
}

bool BadRunCheck(AMSEventR *pev){
  
  if (!pev) return false;//return false to have the cutnotpassed evaluated more reasonably
  
  int tag = pev->IsBadRun("");
  
  if(tag==0) return false;        //good run
  else if(tag==1) return true;    //bad run
  else if(tag==2) {
    printf("Run:%d Ev:%d -> Cannot acces to root_setup. Skipping this event\n", (int)pev->Run(), (int)pev->Event());
    return true;
  }
  return true;    //default
}

int IsScienceRun(AMSEventR *pev){
  
  if(!pev) return -1;//return true to have the cutnotpassed evaluated more reasonably
  
  // Science runtag
  HeaderR* header = &(pev->fHeader);
  if (!header) return -1;//return true to have the cutnotpassed evaluated more reasonably
  
  if (header->RunType>=0xF000) return 1; //faster than: if ((header->RunType >> 12) != 0xf) return true;
  else {
    if (header->RunType>=0xED0D) return 2; //to include TrdRefill
    else if (header->RunType==0xD53C) return 2; //to include TrdRefill
    else if (header->RunType>=0xD510 && header->RunType<=0xD51A) return 2;//to include TrdRefill
  }
  return 0;
}

int RunAnalysisCategory(AMSEventR *pev){
  
  //************
  //This function returns the category of a run
  //return 0:   ''GOLDEN'' is not in the bad run list and has the SCI run flag
  //return 1:   ''SILVER'' is in the bad run list and has the SCI run flag OR without SCI run flag but TRD-usable
  //return -1:  ''NOT USABLE'' run is not SCI or not usable for TRD, so should not be used in the analysis
  //************
  
  //check if the run is "NotUsable" for TRD
  int trd_cat = TrdRunCategory(pev);
  bool trd_unusable = (trd_cat==3);
  
  //check if the run is a "science run" --- consider TrdRefills periods as science runs
  int runtag_cat = IsScienceRun(pev);
  bool sci = (runtag_cat==1 || runtag_cat==2);

  //check if the run is in the bad run list
  //this includes also trd refills runs (however the 'unusable' are checked with the bool below)
  bool bad = (BadRunCheck(pev) || runtag_cat==2);
  
  int tag=-1;//this will be returned if !sci 
  if (trd_unusable) return -1;
  else if (sci) {
    if (bad) tag=1;
    else tag=0; 
  }
  return tag;
}

#endif

int myEvent::RebuildTrigPatt(int &L1TrMemPatt,int &PhysTrPatt){
  //cout<<"Rebuild event...."<<endl;
  //returns: two new trig-patterns and int. flag 0/1->"WasOnlyUnbiasedTrig"/"WasPhysTrig"

  L1TrMemPatt=0;
  PhysTrPatt=0;
  int TrigPatt(0),PhysTrigPatt(0);
  int tofpattft[4]={0,0,0,0};
  int tofpattbz[4]={0,0,0,0};
  int AccPatt(0),NAccs(0);
  int TOFTrigFl1(-1),TOFTrigFl2(-1),ECTrigFl(-1);
  bool FTZ(0);

  int PhysTrSet[8]={0,0,0,0,0,0,0,0};
  PhysTrSet[0]=(1<<1);//ch_unb:FTCP0
  PhysTrSet[1]=(1<<4)+(1<<7);//singl_charged(not el):FTCT1&ACC
  PhysTrSet[2]=(1<<8)+(1<<9);//ions:ACC1&BZ
  PhysTrSet[3]=(1<<5);//sl_ions:FTZ
  PhysTrSet[4]=(1<<4)+(1<<10);//el:FTCT1&EC-F_and
  PhysTrSet[5]=(1<<12);//ph:EC-A_and
  PhysTrSet[6]=(1<<11);//unbEC:EC-A_or
  
  //------> get current trig. parameters:                                                                                                                                          
  for(int ii=0;ii<4;ii++){
    //    tofpattft[ii]=me->TofPatt1[ii];
    // tofpattbz[ii]=me->TofPatt2[ii];
    tofpattft[ii]=TofPatt1[ii];
    tofpattbz[ii]=TofPatt2[ii];
  }
  
  //  TrigPatt=me->JMembPatt;
  //PhysTrigPatt=me->PhysBPatt;
  TrigPatt=JMembPatt;
  PhysTrigPatt=PhysBPatt;
  FTZ=((TrigPatt&(1<<5))!=0);
  
  //AccPatt=me->AntiPatt;
  //TOFTrigFl1=me->TofFlag1;
  //TOFTrigFl2=me->TofFlag2;
  //ECTrigFl=me->EcalFlag;
  
  AccPatt=AntiPatt;
  TOFTrigFl1=TofFlag1;
  TOFTrigFl2=TofFlag2;
  ECTrigFl=EcalFlag;
  //                                N=0/1/2/3-> /undef/noLev1/Lev1or/Lev1and    

bool CentrOK=((tofpattft[2]&0x1FE01FEL)>0);//have any central counter in Lay-3 trig.patt

//-----------------------------------------                                                              

//------>FTC,FTCP0,FTCP1 rebuild:
//
 if(TOFTrigFl1<=4 && TOFTrigFl1>=0){
   L1TrMemPatt|=(1<<1);// set FTCP0(>=3of4)
   L1TrMemPatt|=1;// set FTC
 }

 if(TOFTrigFl1==0)L1TrMemPatt|=(1<<2);// set FTCP1(4of4)
  //-----------------------------------------
  //------>FTCT1 rebuild:
  //
  if(TOFTrigFl1==0 && CentrOK){//was tof 4of4
    L1TrMemPatt|=(1<<4);// set FTCT1
  }
  //-----------------------------------------
  //------> ACC0/1 rebuild:
  //
  for(int i=0;i<8;i++)if((AccPatt & (1<<i))>0)NAccs+=1;//count nsectors
  if(NAccs<1)L1TrMemPatt|=(1<<7);// set ACC0
  if(NAccs<5)L1TrMemPatt|=(1<<8);// set ACC1
  //-----------------------------------------
  //------> FTZ,BZ-bits rebuild(no need to exclude out.c in L3):
  //
  if(FTZ)L1TrMemPatt|=(1<<5);// set FTZ as it was
  if(TOFTrigFl2==0)L1TrMemPatt|=(1<<9);// set BZ "4of4"
  //-----------------------------------------
  //------> EC-bits rebuild:
  //
  if((ECTrigFl/10)==3){
    L1TrMemPatt|=(1<<10);// set EC-F_and
    L1TrMemPatt|=(1<<11);// set EC-F_or
  }
  else if((ECTrigFl/10==1) || (ECTrigFl/10==2))L1TrMemPatt|=(1<<11);// set EC-F_or
  if(ECTrigFl%10>=2){
    //    L1TrMemPatt|=(1<<13);// set EC-A_or
    if(ECTrigFl%10==3)L1TrMemPatt|=(1<<12);// set EC-A_and
  }
  if((L1TrMemPatt&(1<<11))>0)L1TrMemPatt|=(1<<6);// set FTE (EC-F_or is required)
  //-----------------------------------------
  //------>rebuild PhysTrigPatt:
  //
  for(int i=0;i<7;i++){
    if((L1TrMemPatt&PhysTrSet[i])==PhysTrSet[i])PhysTrPatt|=(1<<i);
  }
  //-----------------------------------------

  if((PhysTrPatt&0x3EL)>0)return 1;
  else return 0;
}