RALICE/icepack/iceconvert/IceCal2Root.cxx

   1 /*******************************************************************************
   2  * Copyright(c) 2003, IceCube Experiment at the South Pole. All rights reserved.
   3  *
   4  * Author: The IceCube RALICE-based Offline Project.
   5  * Contributors are mentioned in the code where appropriate.
   6  *
   7  * Permission to use, copy, modify and distribute this software and its
   8  * documentation strictly for non-commercial purposes is hereby granted
   9  * without fee, provided that the above copyright notice appears in all
  10  * copies and that both the copyright notice and this permission notice
  11  * appear in the supporting documentation.
  12  * The authors make no claims about the suitability of this software for
  13  * any purpose. It is provided "as is" without express or implied warranty.
  14  *******************************************************************************/
  15
  16 // $Id$
  17
  18 ///////////////////////////////////////////////////////////////////////////
  19 // Class IceCal2Root
  20 // Conversion of Amanda ascii calibration data into an AliObjMatrix object
  21 // containing the complete OM position, calibration, Xtalk etc... database.
  22 // In addition a PDG particle database, extended with some specific Amanda
  23 // entries, is provided as well.
  24 // This class is derived from AliJob providing task-based processing.
  25 // The main object in the job environment is an AliObjMatrix* pointer
  26 // which contains the OM database.
  27 // Note that the data structures are only written out if an outputfile has
  28 // been specified via the SetOutputFile memberfunction.
  29 // In case no outputfile has been specified, this class provides a facility
  30 // to investigate/use the dbase data directly in subsequent (sub)tasks.
  31 //
  32 // The OM database information in the AliObjMatrix has the following structure :
  33 //
  34 // (j,1)    : Pointer to OM with identifier "j"
  35 // (j,k+1)  : Pointer to a TF1* being the probability function for Xtalk
  36 //            with OM "j" as transmitter and OM "k" as receiver.
  37 //
  38 // The geometry information is directly available from the OM pointer
  39 // in the form of its position and data words like "ORIENT" for orientation etc...
  40 // Just use the OM memberfunction Data() to obtain a full overview.
  41 //
  42 // Note : Position coordinates are indicated in meters and times are in nanoseconds,
  43 //        in accordance with the convention used previously for Amanda.
  44 //
  45 // From the OM pointer also the various (de)calibration functions for
  46 // ADC, LE and TOT can be obtained as TF1* pointers.
  47 // The actual values of the calibration constants are stored as parameters
  48 // of these (de)calibration functions and can be investigated via the
  49 // usual TF1::Print() or TF1::GetParameter() facilities.
  50 // The last two parameters of the Xtalk probability function have no effect
  51 // on the evaluated probability value. However, these two parameters provide
  52 // the minimum and maximum allowed LE differences between the transmitter
  53 // and receiver hits, respectively (as can be seen from the parameter names).
  54 //
  55 // The (de)calibration of signals and/or determination of the Xtalk probability
  56 // can be performed via the standard TF1::Eval(x) functionality, where "x"
  57 // represents the input argument of the function (e.g. an uncalibrated ADC value).
  58 //
  59 // In general the database is not directly accessed by the user in performing
  60 // physics analysis, since all the necessary information is contained in the
  61 // event data itself and available via the GetSignal() memberfunction of the hits.
  62 // However, specific tasks like e.g. calibration, Xtalk correction,
  63 // bad module removal, noise hit removal etc... might need explicit database access.
  64 // So, at the end of the example below some functionality is indicated for clarity.
  65 // The user may use exactly the same procedures to obtain explicit access to the calibration
  66 // functions etc... from the various OMs and/or hits within the actual event data which he/she
  67 // is analysing.
  68 //
  69 // The PDG particle database is a standard ROOT TDatabasePDG object
  70 // with the following extensions :
  71 //
  72 // Name        PDG code
  73 // ----        --------
  74 // brems       10001001
  75 // deltae      10001002
  76 // pairprod    10001003
  77 // nucl_int    10001004
  78 // mu_pair     10001005
  79 // hadrons     10001006
  80 // fiberlaser  10002100
  81 // n2laser     10002101
  82 // yaglaser    10002201
  83 // z_primary   10003000
  84 // a_primary   10003500
  85 //
  86 // Usage example :
  87 // ---------------
  88 //
  89 // gSystem->Load("ralice");
  90 // gSystem->Load("icepack");
  91 // gSystem->Load("iceconvert");
  92 //
  93 // IceCal2Root q("IceCal2Root","Amacalib to IcePack data structure conversion");
  94 //
  95 // // The Amacalib input filename
  96 // q.SetAmacalibFile("amacalib_amanda2_2003.txt");
  97 //
  98 // // Output file for the event structures
  99 // q.SetOutputFile("calib2003.root");
 100 //
 101 // ///////////////////////////////////////////////////////////////////
 102 // // Here the user can specify his/her sub-tasks to be executed
 103 // // after the database structures have been filled and before the
 104 // // data is written out.
 105 // // Sub-tasks (i.e. a user classes derived from TTask) are entered
 106 // // as follows :
 107 // //
 108 // //    MyTask1 task1("task1","Some modifications to specific OMs");
 109 // //    MyTask2 task2("task2","Removal of specific OMs");
 110 // //    MyTask3 task3("task3","Add private objects to the output file");
 111 // //    q.Add(&task1);
 112 // //    q.Add(&task2);
 113 // //    q.Add(&task3);
 114 // //
 115 // // The sub-tasks will be executed in the order as they are entered.
 116 // ///////////////////////////////////////////////////////////////////
 117 //
 118 // // Perform the conversion and execute subtasks (if any)
 119 // q.ExecuteJob();
 120 //
 121 // // Outline of dbase usage for (de)calibration and Xtalk
 122 //
 123 // AliObjMatrix* omdb=q.GetOMdbase();
 124 // IceAOM* om=(IceAOM*)omdb->GetObject(9,1); // Pointer to OM 9
 125 // om->Data(); // Overview of generic module parameters
 126 // TF1* fcal=0;   // Calibration function
 127 // TF1* fdecal=0; // De-calibration function
 128 // fcal=om->GetCalFunction("ADC");
 129 // Float_t adc=248; // Uncalibrated ADC
 130 // Float_t cadc=0;  // Calibrated ADC
 131 // if (fcal) cadc=fcal->Eval(adc);
 132 // fcal=om->GetCalFunction("TOT");
 133 // Float_t tot=1538; // Uncalibrated TOT
 134 // Float_t ctot=0;   // Calibrated TOT
 135 // if (fcal) ctot=fcal->Eval(tot);
 136 // fdecal=om->GetDecalFunction("LE");
 137 // Float_t le=21697; // Uncalibrated LE
 138 // Float_t cle=0;    // Calibrated LE
 139 // if (fcal) cle=fcal->Eval(le);
 140 //
 141 // // Xtalk probability between (trans) OM 90 and (rec) OM 113
 142 // // for a transmitter signal of uncalibrated amplitude "adc".
 143 // TF1* fxtalkp=(TF1*)omdb->GetObject(90,113+1);
 144 // Float_t prob=0;
 145 // adc=378;
 146 // if (fxtalkp) prob=fxtalkp->Eval(adc);
 147 //
 148 //--- Author: Nick van Eijndhoven 09-aug-2005 Utrecht University
 149 //- Modified: NvE $Date$ Utrecht University
 150 ///////////////////////////////////////////////////////////////////////////
 151
 152 #include "IceCal2Root.h"
 153 #include "Rstrstream.h"
 154
 155 ClassImp(IceCal2Root) // Class implementation to enable ROOT I/O
 156
 157 IceCal2Root::IceCal2Root(const char* name,const char* title) : AliJob(name,title)
 158 {
 159 // Default constructor.
 160  fAmacalFileName="";
 161  fRootFileName="";
 162  fOutfile=0;
 163
 164  fPdg=0;
 165  fOmdb=0;
 166 }
 167 ///////////////////////////////////////////////////////////////////////////
 168 IceCal2Root::~IceCal2Root()
 169 {
 170 // Default destructor.
 171
 172  if (fPdg)
 173  {
 174   delete fPdg;
 175   fPdg=0;
 176  }
 177
 178  if (fOmdb)
 179  {
 180   delete fOmdb;
 181   fOmdb=0;
 182  }
 183 }
 184 ///////////////////////////////////////////////////////////////////////////
 185 void IceCal2Root::SetAmacalibFile(TString name)
 186 {
 187 // Set the name of the Amacalib input file.
 188  fAmacalFileName=name;
 189 }
 190 ///////////////////////////////////////////////////////////////////////////
 191 void IceCal2Root::SetOutputFile(TString name)
 192 {
 193 // Set the name of the ROOT output file.
 194  fRootFileName=name;
 195 }
 196 ///////////////////////////////////////////////////////////////////////////
 197 TDatabasePDG* IceCal2Root::GetPDG()
 198 {
 199 // Provide pointer to the PDG database
 200  return fPdg;
 201 }
 202 ///////////////////////////////////////////////////////////////////////////
 203 AliObjMatrix* IceCal2Root::GetOMdbase()
 204 {
 205 // Provide pointer to the OM geometry, calib. etc... database
 206  return fOmdb;
 207 }
 208 ///////////////////////////////////////////////////////////////////////////
 209 void IceCal2Root::Exec(Option_t* opt)
 210 {
 211 // Job to convert the ascii database info into the IcePack structure.
 212 //
 213 // Notes :
 214 // -------
 215 // 1) This class is derived from AliJob, allowing a task based processing.
 216 //    After conversion of the ascii dbase data into the IcePack structure,
 217 //    the processing of all available sub-tasks (if any) is invoked.
 218 //    This provides a facility to investigate/use the dbase data in
 219 //    subsequent (sub)tasks processing before the final data structures
 220 //    are written out.
 221 // 2) The main object in this job environment is an AliObjMatrix* pointer
 222 //    which contains the OM database.
 223
 224  if (fAmacalFileName=="")
 225  {
 226   cout << " *IceCal2Root Exec* No amacalib input file specified." << endl;
 227   return;
 228  }
 229
 230  fInput.clear();
 231  fInput.open(fAmacalFileName.Data());
 232
 233  if (!fInput.good())
 234  {
 235   cout << " *IceCal2Root Exec* Bad input file : " << fAmacalFileName.Data() << endl;
 236   return;
 237  }
 238
 239  if (fOutfile)
 240  {
 241   delete fOutfile;
 242   fOutfile=0;
 243  }
 244  if (fRootFileName != "")
 245  {
 246   fOutfile=new TFile(fRootFileName.Data(),"RECREATE","Calibration data in IcePack structure");
 247  }
 248
 249  // The OM database object
 250  if (fOmdb)
 251  {
 252   fOmdb->Reset();
 253  }
 254  else
 255  {
 256   fOmdb=new AliObjMatrix();
 257   fOmdb->SetNameTitle("Cal-OMDBASE","The OM geometry, calib. etc... database");
 258   fOmdb->SetOwner();
 259  }
 260
 261  // Create the particle database and extend it with some F2000 specific definitions
 262  if (fPdg) delete fPdg;
 263  fPdg=new TDatabasePDG();
 264  fPdg->SetNameTitle("PDG-DBASE","The extended PDG particle database");
 265  Double_t me=fPdg->GetParticle(11)->Mass();
 266  fPdg->AddParticle("brems"   ,"brems"   ,0,1,0,0,"none",10001001,0,0);
 267  fPdg->AddParticle("deltae"  ,"deltae"  ,me,1,0,-3,"Lepton",10001002,0,0);
 268  fPdg->AddParticle("pairprod","pairprod",0,1,0,0,"none",10001003,0,0);
 269  fPdg->AddParticle("nucl_int","nucl_Int",0,1,0,0,"none",10001004,0,0);
 270  fPdg->AddParticle("mu_pair" ,"mu_pair" ,0,1,0,0,"none",10001005,0,0);
 271  fPdg->AddParticle("hadrons" ,"hadrons" ,0,1,0,0,"none",10001006,0,0);
 272  fPdg->AddParticle("fiberlaser","fiberlaser",0,1,0,0,"none",10002100,0,0);
 273  fPdg->AddParticle("n2laser"   ,"n2laser"   ,0,1,0,0,"none",10002101,0,0);
 274  fPdg->AddParticle("yaglaser"  ,"yaglaser"  ,0,1,0,0,"none",10002201,0,0);
 275  fPdg->AddParticle("z_primary","z_primary",0,1,0,0,"none",10003000,0,0);
 276  fPdg->AddParticle("a_primary","a_primary",0,1,0,0,"none",10003500,0,0);
 277
 278  // Initialise the job working environment
 279  SetMainObject(fOmdb);
 280  AddObject(fPdg);
 281  if (fOutfile) AddObject(fOutfile);
 282
 283  cout << " ***" << endl;
 284  cout << " *** Start processing of job " << GetName() << " ***" << endl;
 285  cout << " ***" << endl;
 286  cout << " Amacalib input file : " << fAmacalFileName.Data() << endl;
 287  if (fOutfile) cout << " ROOT output file : " << fOutfile->GetName() << endl;
 288
 289  ListEnvironment();
 290
 291  GetCalibData();
 292
 293  // Invoke all available sub-tasks (if any)
 294  ExecuteTasks(opt);
 295
 296  // Write the datastructures to the output file
 297  if (fOutfile)
 298  {
 299   fOutfile->cd();
 300   if (fOmdb) fOmdb->Write();
 301   if (fPdg) fPdg->Write();
 302  }
 303
 304  // Flush remaining memory resident data to the output file
 305  if (fOutfile) fOutfile->Write();
 306 }
 307 ///////////////////////////////////////////////////////////////////////////
 308 void IceCal2Root::GetCalibData()
 309 {
 310 // Obtain all the geometry, calibration and Xtalk data.
 311
 312  // Prescription of the various (de)calibration functions
 313  TF1 fadccal("fadccal","(x-[1])*[0]");
 314  TF1 fadcdecal("fadcdecal","(x/[0])+[1]");
 315  fadccal.SetParName(0,"BETA-ADC");
 316  fadccal.SetParName(1,"PED-ADC");
 317  fadcdecal.SetParName(0,"BETA-ADC");
 318  fadcdecal.SetParName(1,"PED-ADC");
 319
 320  TF1 ftdccal("ftdccal","(x*[0])-[1]-([0]-1.)*32767.-[2]/sqrt([3])");
 321  TF1 ftdcdecal("ftdcdecal","(x+([0]-1.)*32767.+[1]+[2]/sqrt([3]))/[0]");
 322  ftdccal.SetParName(0,"BETA-TDC");
 323  ftdccal.SetParName(1,"T0");
 324  ftdccal.SetParName(2,"ALPHA-TDC");
 325  ftdccal.SetParName(3,"ADC-SLEW");
 326  ftdcdecal.SetParName(0,"BETA-TDC");
 327  ftdcdecal.SetParName(1,"T0");
 328  ftdcdecal.SetParName(2,"ALPHA-TDC");
 329  ftdcdecal.SetParName(3,"ADC-SLEW");
 330
 331  TF1 ftotcal("ftotcal","x*[0]");
 332  TF1 ftotdecal("ftotdecal","x/[0]");
 333  ftotcal.SetParName(0,"BETA-TOT");
 334  ftotdecal.SetParName(0,"BETA-TOT");
 335
 336  // The cross talk probability function
 337  TF1 fxtalkp("fxtalkp","(1.+[2]-[2]+[3]-[3])/(1.+exp(([0]-x)/[1]))");
 338  fxtalkp.SetParName(0,"C");
 339  fxtalkp.SetParName(1,"B");
 340  fxtalkp.SetParName(2,"dLE-min");
 341  fxtalkp.SetParName(3,"dLE-max");
 342
 343  // The basic OM contents
 344  IceAOM om;
 345
 346  om.SetSlotName("ADC",1);
 347  om.SetSlotName("LE",2);
 348  om.SetSlotName("TOT",3);
 349
 350  om.SetSlotName("TYPE",4);
 351  om.SetSlotName("ORIENT",5);
 352  om.SetSlotName("THRESH",6);
 353  om.SetSlotName("SENSIT",7);
 354  om.SetSlotName("BETA-TDC",8);
 355  om.SetSlotName("T0",9);
 356  om.SetSlotName("ALPHA-TDC",10);
 357  om.SetSlotName("PED-ADC",11);
 358  om.SetSlotName("BETA-ADC",12);
 359  om.SetSlotName("KAPPA-ADC",13);
 360  om.SetSlotName("PED-TOT",14);
 361  om.SetSlotName("BETA-TOT",15);
 362  om.SetSlotName("KAPPA-TOT",16);
 363
 364  fInput.seekg(0); // Position at beginning of file
 365  fInput >> dec;   // Make sure all integers starting with 0 are taken in decimal format
 366
 367  TString s;
 368  Int_t jmod,type,serial,string,ix,iy,iz,ori;
 369  Float_t costh=0;
 370  Float_t thresh=0;
 371  Float_t sensit=1;
 372  Double_t pos[3]={0,0,0};
 373  Float_t ped,beta,alpha,kappa;
 374  Int_t pol;
 375  Float_t totped;
 376  Int_t jtrans,jrec;
 377  Float_t c,b,dlemin,dlemax;
 378  IceAOM* omx=0;
 379  TF1* fcal=0;
 380  TF1* fdecal=0;
 381  while (fInput >> s)
 382  {
 383   if (s == "P") // Read the Geom data
 384   {
 385    fInput >> jmod >> type >> serial >> string >> ix >> iy >> iz >> ori;
 386    omx=(IceAOM*)fOmdb->GetObject(jmod,1);
 387    if (!omx)
 388    {
 389     omx=new IceAOM(om);
 390     omx->SetUniqueID(jmod);
 391     fOmdb->EnterObject(jmod,1,omx);
 392    }
 393    pos[0]=double(ix)/1000.;
 394    pos[1]=double(iy)/1000.;
 395    pos[2]=double(iz)/1000.;
 396    omx->SetPosition(pos,"car");
 397    costh=1;
 398    if (ori==2) costh=-1;
 399    omx->SetSignal(type,4);
 400    omx->SetSignal(costh,5);
 401    omx->SetSignal(thresh,6);
 402    omx->SetSignal(sensit,7);
 403   }
 404   else if (s == "T") // Read the Time calibration constants
 405   {
 406    fInput >> jmod >> ped >> beta >> alpha >> pol;
 407    omx=(IceAOM*)fOmdb->GetObject(jmod,1);
 408    if (!omx)
 409    {
 410     omx=new IceAOM(om);
 411     omx->SetUniqueID(jmod);
 412     fOmdb->EnterObject(jmod,1,omx);
 413    }
 414
 415    omx->SetCalFunction(&ftdccal,2);
 416    omx->SetDecalFunction(&ftdcdecal,2);
 417    omx->SetCalFunction(&ftotcal,3);
 418    omx->SetDecalFunction(&ftotdecal,3);
 419
 420    // Flag time slots of bad OMs as dead and don't provide time (de)calib functions
 421    if (ped<-1e5 || beta<=0 || alpha<0)
 422    {
 423     omx->SetDead(2);
 424     omx->SetDead(3);
 425     omx->SetCalFunction(0,2);
 426     omx->SetDecalFunction(0,2);
 427     omx->SetCalFunction(0,3);
 428     omx->SetDecalFunction(0,3);
 429    }
 430
 431    omx->SetSignal(beta,8);
 432    omx->SetSignal(ped,9);
 433    omx->SetSignal(alpha,10);
 434    omx->SetSignal(beta,15);
 435
 436    fcal=omx->GetCalFunction(2);
 437    fdecal=omx->GetDecalFunction(2);
 438    if (fcal)
 439    {
 440     fcal->SetParameter(0,beta);
 441     fcal->SetParameter(1,ped);
 442     fcal->SetParameter(2,alpha);
 443     fcal->SetParameter(3,1.e20);
 444    }
 445    if (fdecal)
 446    {
 447     fdecal->SetParameter(0,beta);
 448     if (!beta) fdecal->SetParameter(0,1);
 449     fdecal->SetParameter(1,ped);
 450     fdecal->SetParameter(2,alpha);
 451     fdecal->SetParameter(3,1.e20);
 452    }
 453
 454    fcal=omx->GetCalFunction(3);
 455    fdecal=omx->GetDecalFunction(3);
 456    if (fcal)
 457    {
 458     fcal->SetParameter(0,beta);
 459    }
 460    if (fdecal)
 461    {
 462     fdecal->SetParameter(0,beta);
 463    }
 464   }
 465   else if (s == "A") // Read the Amplitude calibration constants
 466   {
 467    fInput >> jmod >> ped >> beta >> totped >> pol;
 468    omx=(IceAOM*)fOmdb->GetObject(jmod,1);
 469    if (!omx)
 470    {
 471     omx=new IceAOM(om);
 472     omx->SetUniqueID(jmod);
 473     fOmdb->EnterObject(jmod,1,omx);
 474    }
 475
 476    omx->SetCalFunction(&fadccal,1);
 477    omx->SetDecalFunction(&fadcdecal,1);
 478
 479    // Flag amplitude slots of bad OMs as dead and don't provide amplitude (de)calib functions
 480    if (ped<-1e5 || beta<=0)
 481    {
 482     omx->SetDead(1);
 483     omx->SetCalFunction(0,1);
 484     omx->SetDecalFunction(0,1);
 485    }
 486    if (totped<-1e5)
 487    {
 488     omx->SetDead(3);
 489     omx->SetCalFunction(0,3);
 490     omx->SetDecalFunction(0,3);
 491    }
 492
 493    omx->SetSignal(ped,11);
 494    omx->SetSignal(beta,12);
 495    omx->SetSignal(totped,14);
 496
 497    fcal=omx->GetCalFunction(1);
 498    fdecal=omx->GetDecalFunction(1);
 499    if (fcal)
 500    {
 501     fcal->SetParameter(0,beta);
 502     fcal->SetParameter(1,ped);
 503    }
 504    if (fdecal)
 505    {
 506     fdecal->SetParameter(0,beta);
 507     if (!beta) fdecal->SetParameter(0,1);
 508     fdecal->SetParameter(1,ped);
 509    }
 510   }
 511   else if (s == "K") // Read the cross talk probability constants
 512   {
 513    fInput >> jtrans >> jrec >> c >> b >> dlemin >> dlemax;
 514    omx=(IceAOM*)fOmdb->GetObject(jtrans,1);
 515    if (!omx)
 516    {
 517     omx=new IceAOM(om);
 518     omx->SetUniqueID(jtrans);
 519     fOmdb->EnterObject(jtrans,1,omx);
 520    }
 521
 522    TF1* fx=new TF1(fxtalkp);
 523    fx->SetParameter(0,c);
 524    if (b)
 525    {
 526     fx->SetParameter(1,b);
 527    }
 528    else
 529    {
 530     fx->SetParameter(1,1);
 531    }
 532    fx->SetParameter(2,dlemin);
 533    fx->SetParameter(3,dlemax);
 534    fOmdb->EnterObject(jtrans,jrec+1,fx);
 535   }
 536   else // Skip this line
 537   {
 538    fInput.ignore(99999,'\n');
 539   }
 540  }
 541 }
 542 ///////////////////////////////////////////////////////////////////////////