StFmsJetFilter.cxx

// @(#)root/eg:$Id: StFmsJetFilter.cxx,v 1.1 2017/03/01 14:58:23 jwebb Exp $
// Author: Victor Perev  17/03/2009

//______________________________________________________________________________
#include "stdlib.h"
#include "math.h"
#include "TMath.h"
#include "TRandom.h"
#include <iostream>

#include "StFmsJetFilter.h"
#include "StGenParticle.h"

using namespace std;


static StFmsJetFilter qwerty;




//______________________________________________________________________________
int StFmsJetFilter::RejectEG(const StGenParticleMaster &ptl) const
{
  //ptl.Print("************** In RejectEG ************** ");
  // Condition: number of tracks etaGate[0]<eta<= etaGate[1] must be bigger  etaGate[2]
  // const static double etaGate[3]={2.5,4.2,30};
  const static double etaGate[3]={2.4,4.4,30};
  //  const static double etaGate[3]={-4.2,-2.5,30};
  const StGenParticle *tk=0;
  int n = ptl.Size();
  

  Float_t se[10] = {0};
  Float_t seH[10] = {0};
  Float_t seE[10] = {0};
  
  Float_t se1[10] = {0};
  Float_t seH1[10] = {0};
  Float_t seE1[10] = {0};

  Float_t se2[10]  = {0};
  Float_t seH2[10] = {0};
  Float_t seE2[10] = {0};
  
  //  Float_t  e1=0;
  //  Float_t e2=0;
  //  Float_t  eta1=0;
  //  Float_t eta2=0;
  int kk=0;
  
  Float_t en, phi, deltaphi, tot_e;
  tot_e =0;
  Float_t M_pi=TMath::Pi();
  Float_t max_e = 0;
  //  Int_t partonPresent = 0;
  //  TRandom *r = new TRandom(); 
  
  for (int i=0;i<n;i++) {
    tk = ptl(i); if (!tk)       continue;
    /*
    if(i==6){
      //      e1 = tk->Energy();
      eta1 = tk->Eta();
    }
    
    if(i==7){
      //      e2 = tk->Energy();
      eta2 = tk->Eta();
    }
    */
    //    if(eta1>2.0) partonPresent = 1;
    //    if(eta2>2.0) partonPresent = 1;
    
    if (tk->Eta() < etaGate[0]) continue;
    if (tk->Eta() > etaGate[1]) continue;
    if (tk->GetStatusCode()!=1) continue;
    
    Int_t isNeutral = 0;
    if(tk->GetPdgCode() ==111 || tk->GetPdgCode() ==221 || tk->GetPdgCode() ==22 || tk->GetPdgCode() ==11) isNeutral =1;
    if(TMath::Abs(tk->GetPdgCode() ==13)|| TMath::Abs(tk->GetPdgCode() ==14))continue;
    
    en =     tk->Energy();
    tot_e += en;
    phi = tk->Phi();
    
    kk = 9;
    
    Float_t phi1 = phi;
    if(phi<0)  phi1 = phi+2*TMath::Pi();
 
    for(int ikk =0; ikk <8; ikk++){
      kk = 9;  
      Float_t ik = ikk;
      if((phi1 >  ik*(M_pi/4.))  && (phi1 < (ik+2.)*(M_pi/4.)))  kk = 7-ikk;
      if(ikk==7) {if((phi1 >  ik*(M_pi/4.))  && (phi1 < 1.*(M_pi/4)))  kk = 7-ikk;}
      
      se2[kk] = se2[kk] + en;    
      if(isNeutral ==1){
        seE2[kk] =  seE2[kk] + en;
      }    
      if(isNeutral ==0){
        seH2[kk] =  seH2[kk] + en;
      }
 }
    
    if(max_e <en)max_e = en;   
    kk = 9;
    
    if((phi< 1.0*(M_pi/4))    && (phi > -1.0*(M_pi/4))) kk = 0;
    if((phi<   0*(M_pi/4))    && (phi > -2.0*(M_pi/4))) kk = 1;
    if((phi< -1.0*(M_pi/4))   && (phi > -3.0*(M_pi/4))) kk = 2;
    if((phi< -2.0*(M_pi/4))   && (phi > -4.0*(M_pi/4))) kk = 3;
    
    if((phi < -3.0*(M_pi/4))  && (phi > 3.0*(M_pi/4))) kk = 4;
    if((phi <  4.0*(M_pi/4))  && (phi > 2.0*(M_pi/4))) kk = 5;
    if((phi <  3.0*(M_pi/4))  && (phi >  1.*(M_pi/4))) kk = 6;
    if((phi <  2.0*(M_pi/4))  && (phi >   0*(M_pi/4))) kk = 7;
    
    se1[kk] = se1[kk] + en;    
    if(isNeutral ==1){
      seE1[kk] =  seE1[kk] + en;
    }
    if(isNeutral ==0){
      seH1[kk] =  seH1[kk] + en;
    }
    
for (int k=0;k<5;k++) {
  deltaphi = (k*(-M_pi/4.)) - phi ;
  if(deltaphi >  M_pi) {deltaphi -= 2.0 * M_pi;}
  if(deltaphi < -1.*M_pi) {deltaphi += 2.0 * M_pi;}
  if (deltaphi< -M_pi/2) deltaphi=deltaphi+2*M_pi;
  deltaphi  =  TMath::Abs(deltaphi);
  
  if(deltaphi < M_pi/4.){
    se[k] = se[k] + en;    
    if(isNeutral ==1){
      seE[k] =  seE[k] + en;
    }    
    if(isNeutral ==0){
      seH[k] =  seH[k] + en;
    }
  }
 }

 
 for (int k=5;k<8;k++) {
    deltaphi = (8-k)*(M_pi/4.) - phi ;
    if(deltaphi >  M_pi) {deltaphi -= 2.0 * M_pi;}
    if(deltaphi < -1.*M_pi) {deltaphi += 2.0 * M_pi;}
    if (deltaphi< -M_pi/2) deltaphi=deltaphi+2*M_pi;
    deltaphi  =  TMath::Abs(deltaphi);
    if(deltaphi < M_pi/4.){
      se[k] = se[k] + en; 
      if(isNeutral ==1){
        seE[k] =  seE[k] + en;
      }
        if(isNeutral ==0){      
          seH[k] =  seH[k] + en;      
        }
    }
 }
 
  }
  
  
 int trig = 0;
 Float_t maxE  =0;
 // Float_t maxEn  =0;
 // Float_t maxHn  =0;
 Float_t maxE1  =0;
 // Float_t maxEn1  =0;
 //Float_t maxHn1  =0;
 
 Float_t maxE2  =0;
 //Float_t maxEn2  =0;
 //Float_t maxHn2  =0;

 Int_t max1 =-1;
 Int_t max2 =-1;
 Int_t max3 =-1;
 
 for (int k=0;k<8;k++) {
   // printf("M1 T,E,H : %4.2f  %4.2f  %4.2f \n",se[k],seE[k],seH[k]);
   //printf("M2 T,E,H : %4.2f  %4.2f  %4.2f \n",se1[k],seE1[k],seH1[k]);
   //printf("M3 T,E,H : %4.2f  %4.2f  %4.2f \n",se2[k],seE2[k],seH2[k]);
   //printf(" -------------------------- \n");
   // cout<<se[k]<<" ---  "<<(seE[k]+seH[k])<<endl;  
   se[k]  = seE[k]  + 0.75*seH[k];
   se1[k] = seE1[k] + 0.75*seH1[k];
   se2[k] = seE2[k] + 0.75*seH2[k];
   
}
 for (int k=0;k<8;k++) {
   // if(se[k]>0&& partonPresent ) trig =1;  // 30 GeV energy in Qudrant
   // if(max_e>60 ) trig =1;  // 30 GeV energy in Qudrant
   if(se[k]>maxE){
     maxE = se[k];
     //maxEn = seE[k];
     //maxHn = seH[k];
     max1 = k;  
   }
   
   if(se1[k]>maxE1){
     maxE1 = se1[k];
     //maxEn1 = seE1[k];
     //maxHn1 = seH1[k];
     max2 = k;  
   }

   if(se2[k]>maxE2){
     maxE2 = se2[k];
     //maxEn2 = seE2[k];
     //maxHn2 = seH2[k];
     max3 = k;   
   }
   
 }
 
 
 if(maxE>=30 ) trig =1;  // 30 GeV energy in Qudrant
 //if(maxE>0&& partonPresent==1 ) trig =1;  // 30 GeV energy in Qudrant
 if (trig==0){
   //   cout<<"NO TRIGGER :p6, p7, el had el+%50Had-quad ToT : "<<e1<<"  "<<e2<<" ("<<maxEn<<"  "<<maxHn<<") "<<maxE<<" --- "<<tot_e<<endl;
   return 1;
 }
 //  printf("  %4.2f  %4.2f  %4.2f  %4.2f  %4.2f  %4.2f ---- 75%  \n",e1,e2,maxEn,maxHn,maxE,tot_e );
 //printf("  %4.2f  %4.2f  %4.2f %d %d %d\n ",maxE, maxE1, maxE2, max1, max2, max3);
 //printf("Trig2 :   %4.2f  %4.2f  %4.2f  %4.2f  %4.2f  %4.2f  \n",e1,e2,maxEn1,maxHn1,maxE1,tot_e);
 //printf("Trig3 :   %4.2f  %4.2f  %4.2f  %4.2f  %4.2f  %4.2f  \n",e1,e2,maxEn2,maxHn2,maxE2,tot_e);
 //cout<<" "<<e1<<"  "<<e2<<" "<<maxEn1<<"  "<<maxHn1<<" "<<maxE1<<" "<<tot_e<<endl;
 //cout<<" "<<e1<<"  "<<e2<<" "<<maxEn2<<"  "<<maxHn2<<" "<<maxE2<<" "<<tot_e<<endl;

 printf("MAX        :    %4.2f  %4.2f  %4.2f %d %d %d\n ",maxE, maxE1, maxE2, max1, max2, max3);
 
 
 return 0;
}

//Float_t StFmsJetFilter::DPhi(Float_t phim, Float_t phi, Float_t &dphi){
Float_t StFmsJetFilter::DPhi(Float_t phim, Float_t phi){
  Float_t dphi;  
  Float_t dp = phim - phi ;
  Float_t M_pi = TMath::Pi();
  if(dp >  M_pi) {dp -= 2.0 * M_pi;}
  if(dp < -1.*M_pi) {dp += 2.0 * M_pi;}
  dphi=dp;
  if (dphi< -M_pi/2) dphi=dphi+2*M_pi;
  dphi  =  TMath::Abs(dphi);
  return dphi;
  
}

//______________________________________________________________________________
int StFmsJetFilter::RejectGT(const StGenParticleMaster &ptl) const
{
  //  ptl.Print("************** In RejectGT ************** ");
  return 0;
}
//______________________________________________________________________________
int StFmsJetFilter::RejectGE(const StGenParticleMaster &ptl) const
{
  //ptl.Print("************** In RejectGE ************** ");
  return 0;
}