BetterMCTables.cxx

Go to the documentation of this file.

#include "CalDetDST/BetterMCTables.h"

#include <cmath>
#include <iostream>

#include "DatabaseInterface/DbiCascader.h"
#include "DatabaseInterface/DbiWriter.h"
#include "Plex/PlexStripEndId.h"
#include "Calibrator/Calibrator.h"
#include "Calibrator/CalMIPCalibration.h"
#include "Validity/VldTimeStamp.h"
#include "Validity/VldContext.h"
#include "Validity/VldRange.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "UgliGeometry/UgliStripHandle.h"

using std::cout;
using std::endl;

float calc_acor(float lc, float lg, float c1, float c2, float l1, float l2, float l3){

   float a=(1.0/(c1+c2))*(c1*std::exp(-lg/l1) + c2*std::exp(-lg/l2));
   float b=std::exp(-lc/l3);
   return (a*b);
}

void load_mc_mipcal(float scale){

// fake up a mipcal table for mc

//  Ask the singleton DbiTableProxyRegistry for the DbiCascader.
  DbiCascader& cascader=DbiTableProxyRegistry::Instance().GetCascader();

//  Define the table.
  string tableDescr ="(SEQNO int,  SEIDKEY int,  STRIPENDID int,  SCALE float)";
  //  Ask the cascader to find a database that will accept it.
  Int_t dbNoTemp = cascader.CreateTemporaryTable("CALMIPCALIBRATION",tableDescr);

  if (dbNoTemp<0){
    cout << "No  database to will accept temporary tables. " << endl;
    return;
  }
  else{
    cout<<"Database number "<<dbNoTemp<<endl;
  }
  VldTimeStamp vts(2001,1,1,0,0,0);
  VldTimeStamp vte(2010,1,1,0,0,0);

  VldRange vr(Detector::kCalDet,SimFlag::kMC,vts,vte,"by hand");
  VldTimeStamp create;
/*
  cout<<"Validity time stamp create: "<<endl;
  create.Print();
  cout<<endl;
  cout<<"Validity time stamp start: "<<endl;
  vts.Print();
  cout<<endl;
  cout<<"Validity time stamp end: "<<endl;
  vte.Print();
  cout<<endl;
*/

  DbiWriter<CalMIPCalibration> tcal(vr, -1, 0, create);

//  cout<<"Starting to loop over planes and strips"<<endl;
  for(int plane=0;plane<60;plane++){
    for(int strip=0;strip<24;strip++){
      PlexStripEndId seidn(Detector::kCalDet,plane,strip,
                           StripEnd::kNegative);
      PlexStripEndId seidp(Detector::kCalDet,plane,strip,
                           StripEnd::kPositive);
      Int_t seidkeyn = seidn.BuildPlnStripEndKey();
      Int_t seidkeyp = seidp.BuildPlnStripEndKey();
      CalMIPCalibration tcn(seidkeyn, StripEnd::kNegative,scale);
      tcal<<tcn;
      CalMIPCalibration tcp(seidkeyp, StripEnd::kPositive,scale);
      tcal<<tcp;
    }
  }
  if(!tcal.CanOutput()){
    cout<<"MakeTimeCal: Writer can't output!"<<endl;
  }
  
  cout<<"Done."<<endl;
  tcal.Close();

}


void load_better_mipcal(float scale, VldTimeStamp vtm){

// fake up a mipcal table for mc

//  Ask the singleton DbiTableProxyRegistry for the DbiCascader.
  DbiCascader& cascader=DbiTableProxyRegistry::Instance().GetCascader();

//  Define the table.
  string tableDescr ="(SEQNO int,  SEIDKEY int,  STRIPENDID int,  SCALE float)";
  //  Ask the cascader to find a database that will accept it.
  Int_t dbNoTemp = cascader.CreateTemporaryTable("CALMIPCALIBRATION",tableDescr);

  if (dbNoTemp<0){
    cout << "No  database to will accept temporary tables. " << endl;
    return;
  }
  else{
    cout<<"Database number "<<dbNoTemp<<endl;
  }

  VldTimeStamp vts(2001,1,1,0,0,0);
  VldTimeStamp vte(2010,1,1,0,0,0);


  VldRange vr(Detector::kCalDet,SimFlag::kMC,vts,vte,"by hand");
  VldTimeStamp create;

  DbiWriter<CalMIPCalibration> tcal(vr, -1, 0, create);
  
  // the time that the mc run corresponds to
  VldContext vcm(Detector::kCalDet, SimFlag::kMC, vtm);
  
  cout<<"better_mc_tables! \nstart: "<<vts
      <<"\nend: "<<vte
      <<"\n vtm: "<<vtm
      <<"\ncontext: "<<vcm<<endl;

  // we will need some geometry information to calculate attenuations
  UgliGeomHandle ugh(vcm);
  
  // get a digit calibrator
  // we will use it to account for gains
  Calibrator::Instance().Reset(vcm);
  const float about_1pe=60;

//  cout<<"Starting to loop over planes and strips"<<endl;
  for(int plane=0;plane<60;plane++){
    for(int strip=0;strip<24;strip++){
      PlexStripEndId seidn(Detector::kCalDet,plane,strip,
                           StripEnd::kNegative);
      PlexStripEndId seidp(Detector::kCalDet,plane,strip,
                           StripEnd::kPositive);
      const Int_t seidkeyn = seidn.BuildPlnStripEndKey();
      const Int_t seidkeyp = seidp.BuildPlnStripEndKey();
      UgliStripHandle ugs_n = ugh.GetStripHandle(seidn);
      UgliStripHandle ugs_p = ugh.GetStripHandle(seidp);
//      cout<<plane<<"  "<<strip<<endl;
      const bool nvalid = ugs_n.IsValid();
      const bool pvalid = ugs_p.IsValid(); 
      if(!nvalid) cout<<"neg not valid!"<<endl;
      if(!pvalid) cout<<"pos not valid!"<<endl;
      if(!(pvalid & nvalid)) continue;

      const float clear_n = ugs_n.ClearFiber(seidn.GetEnd());
      const float green_n = ugs_n.GetHalfLength()
         +ugs_n.WlsPigtail(seidn.GetEnd());

      const float clear_p = ugs_p.ClearFiber(seidp.GetEnd());
      const float green_p = ugs_p.GetHalfLength()
         +ugs_p.WlsPigtail(seidp.GetEnd());
      // fiber attenuations from mid strip to phototube

float calc_acor(float lc, float lg, float c1=0.3333, float c2=0.6667, 
                float l1=1.052, float l2=7.078, float l3=11.4);


      const float an = calc_acor(clear_n, green_n, 
                                 0.3333, 0.6667, 1.052, 7.078, 11.4);
      const float ap = calc_acor(clear_p, green_p, 
                                 0.3333, 0.6667, 1.052, 7.078, 11.4);
      // channel gain
//      cout<<"here"<<endl;
//      cout<<dc.GetPhotoElectrons(about_1pe,seidn)<<endl;
//      cout<<dc.GetPhotoElectrons(about_1pe,seidp)<<endl;
      const float gn = about_1pe/(Calibrator::Instance().GetPhotoElectrons(about_1pe,seidn));
      const float gp = about_1pe/(Calibrator::Instance().GetPhotoElectrons(about_1pe,seidp));
      // normalization for attenuation
      // calculate the average (of both sides) attenuation
      // for a central strip with 4.65m green on one side
      // and 6m clear on the other with 0.65 green
      // this is just a convention but it corresponds to what we
      // usually are looking at for beam particles in 2002
      // scale then is average pe/mip for central strips
      const float norm=(calc_acor(0, 4.65)+calc_acor(6,0.65))/2.0;
      // constants
      const float sn=scale*gn*an/norm;
      const float sp=scale*gp*ap/norm;
      
      cout<<plane<<"  "<<strip<<"  "<<sn<<"  "<<sp<<endl;

      // eventual table values are ~ sigcor/mip
      CalMIPCalibration tcn(seidkeyn, StripEnd::kNegative,sn);
      tcal<<tcn;
      CalMIPCalibration tcp(seidkeyp, StripEnd::kPositive,sp);
      tcal<<tcp;
    }
  }
  if(!tcal.CanOutput()){
    cout<<"MakeTimeCal: Writer can't output!"<<endl;
  }

  cout<<"Done."<<endl;
  tcal.Close();
  cout<<"vtm: "<<vtm<<endl;
}

---