10-aug-2005 NvE AliSignal::GetSignal extended with mode=8 to support dead flag of...
[u/mrichter/AliRoot.git] / 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 ///////////////////////////////////////////////////////////////////////////