TMarSyst.C

Go to the documentation of this file.

// Class TMarSyst
//
// Author     : F. Cassol Brunner
// Created    : 15/1/2004
// Last update: 
//          by: 
// Comment: Class to handle the study of systematic errors

#include "Marana/TMarSyst.h"
#include "Marana/TMarEvent.h"

TMarSyst::TMarSyst()
{
   Sign=0;
   Type=0;
   Name="/";
   Scale_now=NULL;

   MapPossibleScales["/"]=0;
   MapPossibleScales["E_EM_UP"]=+1;       // energy electron scale up
   MapPossibleScales["E_EM_DOWN"]=-1;     // energy electron scale down
   MapPossibleScales["THE_EM_UP"]=+2;     // theta electron scale up
   MapPossibleScales["THE_EM_DOWN"]=-2;   // theta electron scale down
   MapPossibleScales["E_JET_UP"]=+3;      // energy jet scale up
   MapPossibleScales["E_JET_DOWN"]=-3;    // energy jet scale down
   MapPossibleScales["THE_JET_UP"]=+4;    // theta jet scale up
   MapPossibleScales["THE_JET_DOWN"]=-4;  // theta jet scale down
   MapPossibleScales["PT_MU_UP"]=+5;      // pt mu scale up
   MapPossibleScales["PT_MU_DOWN"]=-5;    // pt mu scale down
   MapPossibleScales["THE_MU_UP"]=+6;     // theta mu scale up
   MapPossibleScales["THE_MU_DOWN"]=-6;   // theta mu scale down
   MapPossibleScales["THE_ALL_UP"]=+7;    // theta electron,jet,mu scale up
   MapPossibleScales["THE_ALL_DOWN"]=-7;  // theta electron,jet,mu scale down
   MapPossibleScales["E_ALL_UP"]=+8;      // energy electron,jet and pt mu scale up
   MapPossibleScales["E_ALL_DOWN"]=-8;    // energy electron,jet and pt mu scale down
   MapPossibleScales["ID_EM_UP"]=+9;     // identification electron scale up
   MapPossibleScales["ID_EM_DOWN"]=-9;   // identification electron scale down
   MapPossibleScales["ID_JET_UP"]=+10;     // identification jet scale up
   MapPossibleScales["ID_JET_DOWN"]=-10;   // identification jet scale down
   MapPossibleScales["ID_MU_UP"]=+11;     // identification muon scale up
   MapPossibleScales["ID_MU_DOWN"]=-11;   // identification muon scale down
   MapPossibleScales["ID_PHO_UP"]=+12;     // identification electron scale up
   MapPossibleScales["ID_PHO_DOWN"]=-12;   // identification electron scale down
   
   MapPossibleScales["TRIGGER_UP"]=+13;      // trigger up
   MapPossibleScales["TRIGGER_DOWN"]=-13;    // trigger down
   MapPossibleScales["USERSEL_UP"]=+14;      // related to user selection up
   MapPossibleScales["USERSEL_DOWN"]=-14;    // related to user selection down
   MapPossibleScales["HFSNOISE_UP"]=+15;      // related to user selection up
   MapPossibleScales["HFSNOISE_DOWN"]=-15;    // related to user selection down

//---- for separation between correlated/uncorrelated parts (mainly NC/CC analysis)
   MapPossibleScales["E_EM_UP_COR"]=+16;       // energy electron scale up
   MapPossibleScales["E_EM_DOWN_COR"]=-16;     // energy electron scale down
   MapPossibleScales["E_EM_UP_UNCOR"]=+17;       // energy electron scale up
   MapPossibleScales["E_EM_DOWN_UNCOR"]=-17;     // energy electron scale down
   MapPossibleScales["E_JET_UP_COR"]=+18;      // energy jet scale up
   MapPossibleScales["E_JET_DOWN_COR"]=-18;    // energy jet scale down
   MapPossibleScales["E_JET_UP_UNCOR"]=+19;      // energy jet scale up
   MapPossibleScales["E_JET_DOWN_UNCOR"]=-19;    // energy jet scale down

//--- tau identification
   MapPossibleScales["ID_TAU_UP"]=+20;     // identification muon scale up
   MapPossibleScales["ID_TAU_DOWN"]=-20;   // identification muon scale down

//--- Forward detectors
   MapPossibleScales["EFF_FMD_UP"]=+21;     // Fmd efficiency up
   MapPossibleScales["EFF_FMD_DOWN"]=-21;   // Fmd efficiency down
   MapPossibleScales["EFF_PRT_UP"]=+22;     // PRT efficiency up
   MapPossibleScales["EFF_PRT_DOWN"]=-22;   // PRT efficiency down
   MapPossibleScales["EFF_PLUG_UP"]=+23;     // Plug efficiency up
   MapPossibleScales["EFF_PLUG_DOWN"]=-23;   // Plug efficiency down

//--- MC weighting errors
   MapPossibleScales["W_XPOM_UP"]=+24;     // MC xpomeron weighting  up
   MapPossibleScales["W_XPOM_DOWN"]=-24;   // MC xpomeron weighting  down
   MapPossibleScales["W_BETA_UP"]=+25;     // MC Beta weighting  up
   MapPossibleScales["W_BETA_DOWN"]=-25;   // MC Beta weighting  down
   MapPossibleScales["W_T_UP"]=+26;     // MC T weighting  up
   MapPossibleScales["W_T_DOWN"]=-26;   // MC T weighting  down
   MapPossibleScales["W_Q2_UP"]=+27;     // MC Q2 weighting up
   MapPossibleScales["W_Q2_DOWN"]=-27;   // MC Q2 weighting down



//... default:
   ScalesToBeDone.push_back("/");
};

istream& operator>>(istream& in, TMarSyst& Syst){
   Int_t num=0, out=1;
   char nom[10];
   
   string tag;
 
   
   in >> nom;
   out= sscanf(nom,"%d",&num);
   
   if(!out){
     Error("TMarSyst::istream","Wrong input card format (systematics:). Exit");
     exit(1);
   }
     
   for (Int_t i=0;i<num;i++){
      in >> tag;
      
      if(Syst.MapPossibleScales.find(tag)!=Syst.MapPossibleScales.end())Syst.ScalesToBeDone.push_back(tag);
      else {
        Error("TMarSyst::istream","Warning Scale '%s' not known. Exit",tag.c_str());
        exit(1);
      }     
     
   }
   
   

   return in;
};
TMarSyst::~TMarSyst(){

   ScalesToBeDone.clear();
   MapPossibleScales.clear();
  
}
ostream& operator<<(ostream& out, TMarSyst& Syst){

   for(list<string>::iterator i=Syst.ScalesToBeDone.begin();i!=Syst.ScalesToBeDone.end();i++)
   {
      out << " " << (*i);
   }
   out << endl;

   return out;
};


TMarSyst& TMarSyst::operator= (const TMarSyst& other) 
{
   Sign=other.Sign;
   Type=other.Type;
   Name=other.Name;
   MapPossibleScales=other.MapPossibleScales;
   ScalesToBeDone=other.ScalesToBeDone;
   return *this;

};
Int_t TMarSyst::Next()
{

    if(ScalesToBeDone.size()<=1)return 0;

    if(Scale_now==NULL||Scale_now==ScalesToBeDone.end())Scale_now=ScalesToBeDone.begin();

//... set the type and sign of the scale
    
    Name=(*Scale_now);
    
    int scale=MapPossibleScales[Name];
    Type=abs(scale);      
    if(Type!=0)Sign=Type/(scale); 
    else Sign=0;

//... prepare the pointer for the next scale
    Scale_now++;
     
    return Type;
}

void TMarSyst::ScaleEnergyPartEm(vector<TMarBody>& PartList, const Double_t& Zvtx,const Int_t RunYear
                                ,const Int_t Correlation){

//... syst. on energy scale (see. thesis of B. Heinemann)
//... Move scale from 0.7-3% depending on zCluster (theta)
//... If up-correction use factors 1.007..1.03
//... Correlation is used to separe it into one correlated and uncorrelated part
//... 0 = all, 1 = correlated, 2 = uncorrelated

  Float_t Corr90E    = 0.007;    // 0.7% in energy scale  z<20cm
  Float_t Corr40E    = 0.015;     // 1.5% in energy scale   20<z<100
  Float_t Corr0E     = 0.03;     // 3% in energy scale    z>100
  Float_t CorrESpacal = 0.01;    // 1% in energy scale  for spacal
  Float_t CorrE;

//--- systematics for Hera-2 data
  if(RunYear==TMarEvent::Y0304||RunYear==TMarEvent::Y05) {
        Corr90E    = 0.015;      // 1.5% in energy scale  z<20cm
        Corr40E    = 0.015;       // 1.5% in energy scale   20<z<100
        Corr0E     = 0.03;       // 3% in energy scale  z>100
  }
  if(Correlation==2) { //--- uncorrelated error
        Corr90E    = 0.005;      
        Corr40E    = 0.005;      
        Corr0E     = 0.005;
  }
  
  
//... list of modified electrons
  
  Float_t zCluster=0;
  Double_t NewE=0;
  Double_t NewPt=0;

//... loop over all PartEm
  for(vector<TMarBody>::iterator body=PartList.begin();body!=PartList.end();body++) {

    
//... find the cluster position
    zCluster=100.5/atan(body->GetTheta())+Zvtx;

//... define the correction to be applied if electron in LAr
    if (zCluster<20)
      CorrE = Corr90E;
    else if (zCluster>20 && zCluster<100)
      CorrE = Corr40E;
    else 
      CorrE = Corr0E;

//... For spacal electron: up to now crude theta angle definition
    if(body->GetTheta()*MYR2D>153.) CorrE = CorrESpacal;      

//... correct energy and Pt
    
    NewE = body->GetE()*(1+(Sign*CorrE));
    NewPt = NewE*TMath::Sin(body->GetTheta());
    
    body->SetE(NewE);
    body->SetPt(NewPt);
  }

  return ;
}

void TMarSyst::ScaleThetaPartEm(vector<TMarBody>& PartList,const Int_t RunYear){


  Float_t CorrDiffThe= 0.003;   // 3 mrad in theta
  const Float_t CorrDiffTheSpacal= 0.001;   // 1 mrad in theta for spacal (if BDC, best is 0.5 mrad)

  if(RunYear==TMarEvent::Y0304||RunYear==TMarEvent::Y05) {
        CorrDiffThe= 0.005;     
  }


  Float_t NewThe=0;  
  Float_t NewEta=0;
  Float_t NewPt=0;

//... loop over all PartEm
  for(vector<TMarBody>::iterator body=PartList.begin();body!=PartList.end();body++) {

   Float_t ThisCorrDiffThe=CorrDiffThe;
   if(body->GetTheta()*MYR2D>153.) ThisCorrDiffThe = CorrDiffTheSpacal;
   
   if(body->GetTheta()*MYR2D>90.0)
        NewThe = acos(cos(body->GetTheta()-Sign*ThisCorrDiffThe));
    else
        NewThe = acos(cos(body->GetTheta()+Sign*ThisCorrDiffThe));

//     NewThe = acos(cos(body->GetTheta()+Sign*CorrDiffThe));
    
    NewEta = -log(tan(NewThe/2));
    NewPt = body->GetE()*TMath::Sin(NewThe);

//... correct Eta et Pt
    
    body->SetEta(NewEta);
    body->SetPt(NewPt);
  }

  return;
}

void TMarSyst::ScaleEnergyPartJet(vector<TMarBody>& PartList,TLorentzVector& HadJetVec
                ,TLorentzVector& HadNoJetVec,const Int_t RunYear,const Int_t Correlation)
{
//... Correlation is used to separe it into one correlated and uncorrelated part
//... 0 = all, 1 = correlated, 2 = uncorrelated
//... Pth is used to check if jet calibration is applied or not 
//                      (ie: at low Q2: 5% systematic error is used)

//... scale factor
  Float_t CorrE    = 0.02;
  Float_t CorrESpacal    = 0.07;   //--- 7% if spacal hadrons
  if(Correlation==1) CorrE = 0.017; //--- correlated part
  if(Correlation==2) CorrE = 0.01; //--- uncorrelated part
  if((HadNoJetVec+HadJetVec).Pt()<8) CorrE = 0.05; //--- 5% applied at low Pt

//... initialize variables
  Float_t NewE=0;  
  Float_t NewPt=0;
  HadJetVec.SetPtEtaPhiE(0,0,0,0);

//... loop over all PartJet
  for(vector<TMarBody>::iterator body=PartList.begin();body!=PartList.end();body++) {
    
    Float_t ThisCorrE=CorrE;
    if(body->GetTheta()*MYR2D>153.) ThisCorrE=CorrESpacal;

//... correct energy and Pt
    
    NewE = body->GetE()*(1+(Sign*ThisCorrE));
    NewPt = NewE*TMath::Sin(body->GetTheta());
    
    body->SetE(NewE);
    body->SetPt(NewPt); 

    HadJetVec+=body->GetFourVector();
  }
//... scale no also the hfs not in jets
  if(HadNoJetVec.Theta()*MYR2D>153.) CorrE=CorrESpacal;
  HadNoJetVec=(1+(Sign*CorrE))*HadNoJetVec;


}
void TMarSyst::ScaleThetaPartJet(vector<TMarBody>& PartList,TLorentzVector& HadJetVec){

//... scale factor
  const Float_t CorrDiffTheC= 0.010; // Diff de 10 mrad central
  const Float_t CorrDiffTheF= 0.005; // Diff de 5 mrad forward

//... initialize variables
  Float_t CorrDiffThe = 0.0;
  Float_t NewThe=0;  
  Float_t NewPt=0; 
  Float_t NewEta=0; 
  Float_t OldThe=0; 
  HadJetVec.SetPtEtaPhiE(0,0,0,0);

//... loop over all PartJet
  for(vector<TMarBody>::iterator body=PartList.begin();body!=PartList.end();body++) {
 
    OldThe=body->GetTheta();

//... Define Error on thetajet
    if (OldThe*MYR2D>20.0){
      CorrDiffThe = CorrDiffTheC;
    }
    else{
      CorrDiffThe = CorrDiffTheF;
    }
    
    if(OldThe*MYR2D>90.0)
        NewThe = acos(cos(OldThe-Sign*CorrDiffThe));
    else
        NewThe = acos(cos(OldThe+Sign*CorrDiffThe));
//     NewThe = acos(cos(OldThe+Sign*CorrDiffThe));
     
    NewEta = -log(tan(NewThe/2));
    NewPt = body->GetE()*TMath::Sin(NewThe);

//... correct Eta et Pt
    
    body->SetEta(NewEta);
    body->SetPt(NewPt);

    HadJetVec+=body->GetFourVector();
  }

}
void TMarSyst::ScalePtPartMuon(vector<TMarBody>& PartList){

//... scale factor
  const Float_t CorrPt    = 0.05;

  Float_t NewE=0;  
  Float_t NewPt=0;


  for(vector<TMarBody>::iterator body=PartList.begin();body!=PartList.end();body++) {

//... correct energy and Pt
    
    NewPt = body->GetPt()*(1+(Sign*CorrPt));
    NewE = NewPt/TMath::Sin(body->GetTheta());
    

    body->SetPt(NewPt);
    body->SetE(NewE);
  }

}
void TMarSyst::ScaleThetaPartMuon(vector<TMarBody>& PartList,const Int_t RunYear){

//... scale factor
  Float_t CorrDiffThe= 0.003; // Diff de 3 mrad central
  if(RunYear==TMarEvent::Y0304||RunYear==TMarEvent::Y05) {
        CorrDiffThe= 0.005;     
  }

  Float_t NewThe=0;  
  Float_t NewPt=0; 
  Float_t NewEta=0; 
  Float_t OldThe=0; 

//... loop over all PartJet
  for(vector<TMarBody>::iterator body=PartList.begin();body!=PartList.end();body++) {
  
    OldThe=body->GetTheta();
    
    if(OldThe*MYR2D>90.0)
        NewThe = acos(cos(OldThe-Sign*CorrDiffThe));
    else
        NewThe = acos(cos(OldThe+Sign*CorrDiffThe));
//     NewThe = acos(cos(OldThe+Sign*CorrDiffThe));
     
    NewEta = -log(tan(NewThe/2));
    NewPt = body->GetE()*TMath::Sin(NewThe);

//... correct Eta et Pt
    
    body->SetEta(NewEta);
    body->SetPt(NewPt);

  }

}

void TMarSyst::ScaleNoiseHFS(TLorentzVector& HadNoJetVec)
{
//      Shift the part of the hfs with no jets by 4-vector of the removed 
//      noise to calculate systematic error comming from the noise substraction
        
//... scale factor
  const Float_t CorrNoise  = 0.10; //--- 25% of syst error on noise reduction, 10% now (Benji,8/07/04)

 static H1FloatPtr PtrNoisePt("HfsClusNoisePt");   
 static H1FloatPtr PtrNoiseTheta("HfsClusNoiseTheta");   
 static H1FloatPtr PtrNoisePhi("HfsClusNoisePhi");   
 static H1FloatPtr PtrNoiseE("HfsClusNoiseE");  


//--- Get the noise 4-vector from Hat
 TLorentzVector NoiseVec(0.,0.,0.,0.);
 NoiseVec.SetPtEtaPhiE(*PtrNoisePt,-TMath::Log(tan(*PtrNoiseTheta/2)),*PtrNoisePhi,*PtrNoiseE);

//... scale no also the hfs not in jets
  HadNoJetVec=HadNoJetVec+Sign*CorrNoise*NoiseVec;


}

---