IceCal2Root.cxx

Go to the documentation of this file.

 
#include "IceCal2Root.h"
#include "Rstrstream.h"
 
ClassImp(IceCal2Root); // Class implementation to enable ROOT I/O
 
IceCal2Root::IceCal2Root(const char* name,const char* title) : NcJob(name,title)
{
 fAmacalFileName="";
 fTWRDaqFileName="";
 fRootFileName="";
 fOutfile=0;
 
 fPdg=0;
 fMuDaqdb=0;
 fTWRDaqdb=0;
}

IceCal2Root::~IceCal2Root()
{
 
 if (fPdg)
 {
  delete fPdg;
  fPdg=0;
 }
 
 if (fMuDaqdb)
 {
  delete fMuDaqdb;
  fMuDaqdb=0;
 }
 
 if (fTWRDaqdb)
 {
  delete fTWRDaqdb;
  fTWRDaqdb=0;
 }
}

void IceCal2Root::SetAmacalibFile(TString name)
{
 fAmacalFileName=name;
}

void IceCal2Root::SetTWRDaqFile(TString name)
{
 fTWRDaqFileName=name;
}

void IceCal2Root::SetOutputFile(TString name)
{
 fRootFileName=name;
}

TDatabasePDG* IceCal2Root::GetPDG()
{
 return fPdg;
}

NcObjMatrix* IceCal2Root::GetOMdbase(TString name)
{
 
 if (name=="MuDaq") return fMuDaqdb;
 if (name=="TWRDaq") return fTWRDaqdb;
 return 0;
}

void IceCal2Root::Exec(Option_t* opt)
{
 
 if (fOutfile)
 {
  delete fOutfile;
  fOutfile=0;
 }
 if (fRootFileName != "")
 {
  fOutfile=new TFile(fRootFileName.Data(),"RECREATE","Calibration data in IcePack structure");
 }
 
 // Create the particle database and extend it with some F2000 specific definitions
 if (fPdg) delete fPdg;
 fPdg=new TDatabasePDG();
 fPdg->SetNameTitle("PDG-DBASE","The extended PDG particle database");
 Double_t me=fPdg->GetParticle(11)->Mass();
 fPdg->AddParticle("brems"   ,"brems"   ,0,1,0,0,"none",10001001,0,0);
 fPdg->AddParticle("deltae"  ,"deltae"  ,me,1,0,-3,"Lepton",10001002,0,0);
 fPdg->AddParticle("pairprod","pairprod",0,1,0,0,"none",10001003,0,0);
 fPdg->AddParticle("nucl_int","nucl_Int",0,1,0,0,"none",10001004,0,0);
 fPdg->AddParticle("mu_pair" ,"mu_pair" ,0,1,0,0,"none",10001005,0,0);
 fPdg->AddParticle("hadrons" ,"hadrons" ,0,1,0,0,"none",10001006,0,0);
 fPdg->AddParticle("fiberlaser","fiberlaser",0,1,0,0,"none",10002100,0,0);
 fPdg->AddParticle("n2laser"   ,"n2laser"   ,0,1,0,0,"none",10002101,0,0);
 fPdg->AddParticle("yaglaser"  ,"yaglaser"  ,0,1,0,0,"none",10002201,0,0);
 fPdg->AddParticle("z_primary","z_primary",0,1,0,0,"none",10003000,0,0);
 fPdg->AddParticle("a_primary","a_primary",0,1,0,0,"none",10003500,0,0);
 
 // Initialise the job working environment
 AddObject(fPdg);
 if (fOutfile) AddObject(fOutfile);
 
 cout << " ***" << endl;
 cout << " *** Start processing of job " << GetName() << " ***" << endl;
 cout << " ***" << endl;
 cout << " Amacalib MuDaq input file : " << fAmacalFileName.Data() << endl;
 cout << " TWRDaq input file : " << fTWRDaqFileName.Data() << endl;
 if (fOutfile) cout << " ROOT output file : " << fOutfile->GetName() << endl;
 
 GetMuDaqData();
 if (fMuDaqdb) AddObject(fMuDaqdb);
 
 GetTWRDaqData();
 if (fTWRDaqdb) AddObject(fTWRDaqdb);
 
 ListEnvironment();
 
 // Invoke all available sub-tasks (if any)
 CleanTasks();
 ExecuteTasks(opt);
 
 // Write the datastructures to the output file
 if (fOutfile)
 {
  fOutfile->cd();
  if (fMuDaqdb) fMuDaqdb->Write();
  if (fTWRDaqdb) fTWRDaqdb->Write();
  if (fPdg) fPdg->Write();
 }
 
 // Flush remaining memory resident data to the output file
 if (fOutfile) fOutfile->Write();
}

void IceCal2Root::GetMuDaqData()
{
 
 if (fAmacalFileName=="")
 {
  cout << " *IceCal2Root GetMuDaqData* No amacalib input file specified." << endl;
  return;
 }
 
 fInput.clear();
 fInput.open(fAmacalFileName.Data());
 
 if (!fInput.good())
 {
  cout << " *IceCal2Root GetMuDaqData* Bad input file : " << fAmacalFileName.Data() << endl;
  return;
 }
 
 // The MuDaq OM database object
 if (fMuDaqdb)
 {
  fMuDaqdb->Reset();
 }
 else
 {
  fMuDaqdb=new NcObjMatrix();
  fMuDaqdb->SetNameTitle("MuDaq-OMDBASE","The MuDaq OM geometry, calib. etc... database");
  fMuDaqdb->SetOwner();
 }
 
 // Prescription of the various (de)calibration functions
 TF1 fadccal("fadccal","(x-[1])*[0]");
 TF1 fadcdecal("fadcdecal","(x/[0])+[1]");
 fadccal.SetParName(0,"BETA-ADC");
 fadccal.SetParName(1,"PED-ADC");
 fadcdecal.SetParName(0,"BETA-ADC");
 fadcdecal.SetParName(1,"PED-ADC");
 
 TF1 ftdccal("ftdccal","(x*[0])-[1]-([0]-1.)*32767.-[2]/sqrt([3])");
 TF1 ftdcdecal("ftdcdecal","(x+([0]-1.)*32767.+[1]+[2]/sqrt([3]))/[0]");
 ftdccal.SetParName(0,"BETA-TDC");
 ftdccal.SetParName(1,"T0");
 ftdccal.SetParName(2,"ALPHA-TDC");
 ftdccal.SetParName(3,"ADC-SLEW");
 ftdcdecal.SetParName(0,"BETA-TDC");
 ftdcdecal.SetParName(1,"T0");
 ftdcdecal.SetParName(2,"ALPHA-TDC");
 ftdcdecal.SetParName(3,"ADC-SLEW");
 
 TF1 ftotcal("ftotcal","x*[0]");
 TF1 ftotdecal("ftotdecal","x/[0]");
 ftotcal.SetParName(0,"BETA-TOT");
 ftotdecal.SetParName(0,"BETA-TOT");
 
 // The cross talk probability function
 TF1 fxtalkp("fxtalkp","(1.+[2]-[2]+[3]-[3])/(1.+exp(([0]-x)/[1]))");
 fxtalkp.SetParName(0,"C");
 fxtalkp.SetParName(1,"B");
 fxtalkp.SetParName(2,"dLE-min");
 fxtalkp.SetParName(3,"dLE-max");
 
 // The basic OM contents
 IceAOM om;
 
 // Slots to hold the various (de)calibration functions
 om.AddNamedSlot("ADC");
 om.AddNamedSlot("LE");
 om.AddNamedSlot("TOT");
 // Slots with hardware parameters
 om.AddNamedSlot("TYPE");
 om.AddNamedSlot("ORIENT");
//@@@ om.AddNamedSlot("THRESH");
//@@@ om.AddNamedSlot("SENSIT");
//@@@ om.AddNamedSlot("READOUT"); // 0=unknown 1=electrical 2=optical 3=digital
 
 fInput.seekg(0); // Position at beginning of file
 fInput >> dec;   // Make sure all integers starting with 0 are taken in decimal format
 
 TString s;
 Int_t jmod,type,serial,string,ix,iy,iz,ori;
 Float_t costh=0;
//@@@ Float_t thresh=0;
//@@@ Float_t sensit=1;
//@@@ Float_t readout=0;
 Double_t pos[3]={0,0,0};
 Float_t ped,beta,alpha;
 Int_t pol;
 Float_t totped;
 Int_t jtrans,jrec;
 Float_t c,b,dlemin,dlemax;
 IceAOM* omx=0;
 TF1* fcal=0;
 TF1* fdecal=0;
 while (fInput >> s)
 {
  if (s == "P") // Read the Geom data
  {
   fInput >> jmod >> type >> serial >> string >> ix >> iy >> iz >> ori;
   omx=(IceAOM*)fMuDaqdb->GetObject(jmod,1);
   if (!omx)
   {
    omx=new IceAOM(om);
    omx->SetUniqueID(jmod);
    fMuDaqdb->EnterObject(jmod,1,omx);
   }
   pos[0]=double(ix)/1000.;
   pos[1]=double(iy)/1000.;
   pos[2]=double(iz)/1000.;
   omx->SetPosition(pos,"car");
   costh=1;
   if (ori==2) costh=-1;
   omx->SetSignal(type,"TYPE");
   omx->SetSignal(costh,"ORIENT");
//@@@   omx->SetSignal(thresh,"THRESH");
//@@@   omx->SetSignal(sensit,"SENSIT");
//@@@   omx->SetSignal(readout,"READOUT");
  }
  else if (s == "T") // Read the Time calibration constants
  {
   fInput >> jmod >> ped >> beta >> alpha >> pol;
   omx=(IceAOM*)fMuDaqdb->GetObject(jmod,1);
   if (!omx)
   {
    omx=new IceAOM(om);
    omx->SetUniqueID(jmod);
    fMuDaqdb->EnterObject(jmod,1,omx);
   }
 
   omx->SetCalFunction(&ftdccal,"LE");
   omx->SetDecalFunction(&ftdcdecal,"LE");
   omx->SetCalFunction(&ftotcal,"TOT");
   omx->SetDecalFunction(&ftotdecal,"TOT");
 
   // Flag time slots of bad OMs as dead and don't provide time (de)calib functions
   if (ped<-1e5 || beta<=0 || alpha<0)
   {
    omx->SetDead("LE");
    omx->SetDead("TOT");
    omx->SetCalFunction(0,"LE");
    omx->SetDecalFunction(0,"LE");
    omx->SetCalFunction(0,"TOT");
    omx->SetDecalFunction(0,"TOT");
   }
 
   fcal=omx->GetCalFunction("LE");
   fdecal=omx->GetDecalFunction("LE");
   if (fcal)
   {
    fcal->SetParameter(0,beta);
    fcal->SetParameter(1,ped);
    fcal->SetParameter(2,alpha);
    fcal->SetParameter(3,1.e20);
   }
   if (fdecal)
   {
    fdecal->SetParameter(0,beta);
    if (!beta) fdecal->SetParameter(0,1);
    fdecal->SetParameter(1,ped);
    fdecal->SetParameter(2,alpha);
    fdecal->SetParameter(3,1.e20);
   }
 
   fcal=omx->GetCalFunction("TOT");
   fdecal=omx->GetDecalFunction("TOT");
   if (fcal)
   {
    fcal->SetParameter(0,beta);
   }
   if (fdecal)
   {
    fdecal->SetParameter(0,beta);
   }
  }
  else if (s == "A") // Read the Amplitude calibration constants
  {
   fInput >> jmod >> ped >> beta >> totped >> pol;
   omx=(IceAOM*)fMuDaqdb->GetObject(jmod,1);
   if (!omx)
   {
    omx=new IceAOM(om);
    omx->SetUniqueID(jmod);
    fMuDaqdb->EnterObject(jmod,1,omx);
   }
 
   omx->SetCalFunction(&fadccal,"ADC");
   omx->SetDecalFunction(&fadcdecal,"ADC");
 
   // Flag amplitude slots of bad OMs as dead and don't provide amplitude (de)calib functions
   if (ped<-1e5 || beta<=0)
   {
    omx->SetDead("ADC");
    omx->SetCalFunction(0,"ADC");
    omx->SetDecalFunction(0,"ADC");
   }
   if (totped<-1e5)
   {
    omx->SetDead("TOT");
    omx->SetCalFunction(0,"TOT");
    omx->SetDecalFunction(0,"TOT");
   }
 
   fcal=omx->GetCalFunction("ADC");
   fdecal=omx->GetDecalFunction("ADC");
   if (fcal)
   {
    fcal->SetParameter(0,beta);
    fcal->SetParameter(1,ped);
   }
   if (fdecal)
   {
    fdecal->SetParameter(0,beta);
    if (!beta) fdecal->SetParameter(0,1);
    fdecal->SetParameter(1,ped);
   }
  }
  else if (s == "K") // Read the cross talk probability constants
  {
   fInput >> jtrans >> jrec >> c >> b >> dlemin >> dlemax;
   omx=(IceAOM*)fMuDaqdb->GetObject(jtrans,1);
   if (!omx)
   {
    omx=new IceAOM(om);
    omx->SetUniqueID(jtrans);
    fMuDaqdb->EnterObject(jtrans,1,omx);
   }
 
   TF1* fx=new TF1(fxtalkp);
   fx->SetParameter(0,c);
   if (b)
   {
    fx->SetParameter(1,b);
   }
   else
   {
    fx->SetParameter(1,1);
   }
   fx->SetParameter(2,dlemin);
   fx->SetParameter(3,dlemax);
   fMuDaqdb->EnterObject(jtrans,jrec+1,fx);
  }
  else // Skip this line
  {
   fInput.ignore(99999,'\n');
  } 
 }
 
 fInput.close();
}

void IceCal2Root::GetTWRDaqData()
{
 
 if (fTWRDaqFileName=="")
 {
  cout << " *IceCal2Root GetTWRDaqData* No TWRDaq calibration data will be produced." << endl;
  return;
 }
 
 fInput.clear();
 fInput.open(fTWRDaqFileName.Data());
 
 if (!fInput.good())
 {
  cout << " *IceCal2Root GetTWRDaqData* Bad input file : " << fTWRDaqFileName.Data() << endl;
  return;
 }
 
 // The geometry info will be obtained from the MuDaq OM database
 if (!fMuDaqdb)
 {
  cout << " *IceCal2Root GetTWRDaqData* MuDaq OM geometry database is missing." << endl;
  return;
 }
 
 // The TWRDaq OM database object
 if (fTWRDaqdb)
 {
  fTWRDaqdb->Reset();
 }
 else
 {
  fTWRDaqdb=new NcObjMatrix();
  fTWRDaqdb->SetNameTitle("TWRDaq-OMDBASE","The TWRDaq OM geometry, calib. etc... database");
  fTWRDaqdb->SetOwner();
 }
 
 // Prescription of the various (de)calibration functions
 TF1 fadccal("fadccal","x*(5./4096.)/(50.*[0])");
 TF1 fadcdecal("fadcdecal","x*(50.*[0])/(5./4096.)");
 fadccal.SetParName(0,"nC/PE");
 fadcdecal.SetParName(0,"nC/PE");
 
 TF1 ftdccal("ftdccal","x-[0]");
 TF1 ftdcdecal("ftdcdecal","x+[0]");
 ftdccal.SetParName(0,"T0");
 ftdcdecal.SetParName(0,"T0");
 
 TF1 ftotcal("ftotcal","x*[0]");
 TF1 ftotdecal("ftotdecal","x/[0]");
 ftotcal.SetParName(0,"TOT-FACT");
 ftotdecal.SetParName(0,"TOT-FACT");
 
 fInput.seekg(0); // Position at beginning of file
 fInput >> dec;   // Make sure all integers starting with 0 are taken in decimal format
 
 Int_t jmod;
 Float_t t0;
 Float_t ncpe;
 IceAOM* omx=0;
 TF1* fcal=0;
 TF1* fdecal=0;
 while (fInput >> jmod >> t0 >> ncpe)
 {
  // Copy the Geom data from the MuDaq OM database
  IceAOM* omg=(IceAOM*)fMuDaqdb->GetObject(jmod,1);
  if (!omg) continue;
 
  omx=new IceAOM(*omg);
 
  // Reset all previous "dead" flags
  omx->SetAlive("ADC");
  omx->SetAlive("LE");
  omx->SetAlive("TOT");
 
  // Enter the TWRDaq (de)calibration functions
  omx->SetCalFunction(&fadccal,"ADC");
  omx->SetDecalFunction(&fadcdecal,"ADC");
  omx->SetCalFunction(&ftdccal,"LE");
  omx->SetDecalFunction(&ftdcdecal,"LE");
  omx->SetCalFunction(&ftotcal,"TOT");
  omx->SetDecalFunction(&ftotdecal,"TOT");
 
  // Flag slots of bad OMs as dead and don't provide time (de)calib functions
  if (ncpe<=0)
  {
   omx->SetDead("ADC");
   omx->SetCalFunction(0,"ADC");
   omx->SetDecalFunction(0,"ADC");
  }
  if (t0<-999)
  {
   omx->SetDead("LE");
   omx->SetDead("TOT");
   omx->SetCalFunction(0,"LE");
   omx->SetDecalFunction(0,"LE");
   omx->SetCalFunction(0,"TOT");
   omx->SetDecalFunction(0,"TOT");
  }
 
  // Set the TWRDaq (de)calibration function parameters for the good OMs
  fcal=omx->GetCalFunction("ADC");
  fdecal=omx->GetDecalFunction("ADC");
  if (fcal)
  {
   fcal->SetParameter(0,ncpe);
  }
  if (fdecal)
  {
   fdecal->SetParameter(0,ncpe);
  }
 
  fcal=omx->GetCalFunction("LE");
  fdecal=omx->GetDecalFunction("LE");
  if (fcal)
  {
   fcal->SetParameter(0,t0);
  }
  if (fdecal)
  {
   fdecal->SetParameter(0,t0);
  }
 
  fcal=omx->GetCalFunction("TOT");
  fdecal=omx->GetDecalFunction("TOT");
  if (fcal)
  {
   fcal->SetParameter(0,1);
  }
  if (fdecal)
  {
   fdecal->SetParameter(0,1);
  }
 
  fTWRDaqdb->EnterObject(jmod,1,omx);
 
 } // End of reading loop
 
 fInput.close();
}