20-sep-2005 NvE Dummy values 0 set for KAPPA-ADC and KAPPA-TOT for OMs in IceCal2Root
[u/mrichter/AliRoot.git] / RALICE / icepack / iceconvert / IceCal2Root.cxx
CommitLineData
216d1d91 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
155ClassImp(IceCal2Root) // Class implementation to enable ROOT I/O
156
157IceCal2Root::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///////////////////////////////////////////////////////////////////////////
168IceCal2Root::~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///////////////////////////////////////////////////////////////////////////
185void IceCal2Root::SetAmacalibFile(TString name)
186{
187// Set the name of the Amacalib input file.
188 fAmacalFileName=name;
189}
190///////////////////////////////////////////////////////////////////////////
191void IceCal2Root::SetOutputFile(TString name)
192{
193// Set the name of the ROOT output file.
194 fRootFileName=name;
195}
196///////////////////////////////////////////////////////////////////////////
197TDatabasePDG* IceCal2Root::GetPDG()
198{
199// Provide pointer to the PDG database
200 return fPdg;
201}
202///////////////////////////////////////////////////////////////////////////
203AliObjMatrix* IceCal2Root::GetOMdbase()
204{
205// Provide pointer to the OM geometry, calib. etc... database
206 return fOmdb;
207}
208///////////////////////////////////////////////////////////////////////////
209void 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///////////////////////////////////////////////////////////////////////////
308void 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);
a4b77ddf 435 omx->SetSignal(0,16);
216d1d91 436
437 fcal=omx->GetCalFunction(2);
438 fdecal=omx->GetDecalFunction(2);
439 if (fcal)
440 {
441 fcal->SetParameter(0,beta);
442 fcal->SetParameter(1,ped);
443 fcal->SetParameter(2,alpha);
444 fcal->SetParameter(3,1.e20);
445 }
446 if (fdecal)
447 {
448 fdecal->SetParameter(0,beta);
449 if (!beta) fdecal->SetParameter(0,1);
450 fdecal->SetParameter(1,ped);
451 fdecal->SetParameter(2,alpha);
452 fdecal->SetParameter(3,1.e20);
453 }
454
455 fcal=omx->GetCalFunction(3);
456 fdecal=omx->GetDecalFunction(3);
457 if (fcal)
458 {
459 fcal->SetParameter(0,beta);
460 }
461 if (fdecal)
462 {
463 fdecal->SetParameter(0,beta);
464 }
465 }
466 else if (s == "A") // Read the Amplitude calibration constants
467 {
468 fInput >> jmod >> ped >> beta >> totped >> pol;
469 omx=(IceAOM*)fOmdb->GetObject(jmod,1);
470 if (!omx)
471 {
472 omx=new IceAOM(om);
473 omx->SetUniqueID(jmod);
474 fOmdb->EnterObject(jmod,1,omx);
475 }
476
477 omx->SetCalFunction(&fadccal,1);
478 omx->SetDecalFunction(&fadcdecal,1);
479
480 // Flag amplitude slots of bad OMs as dead and don't provide amplitude (de)calib functions
481 if (ped<-1e5 || beta<=0)
482 {
483 omx->SetDead(1);
484 omx->SetCalFunction(0,1);
485 omx->SetDecalFunction(0,1);
486 }
487 if (totped<-1e5)
488 {
489 omx->SetDead(3);
490 omx->SetCalFunction(0,3);
491 omx->SetDecalFunction(0,3);
492 }
493
494 omx->SetSignal(ped,11);
495 omx->SetSignal(beta,12);
a4b77ddf 496 omx->SetSignal(0,13);
216d1d91 497 omx->SetSignal(totped,14);
498
499 fcal=omx->GetCalFunction(1);
500 fdecal=omx->GetDecalFunction(1);
501 if (fcal)
502 {
503 fcal->SetParameter(0,beta);
504 fcal->SetParameter(1,ped);
505 }
506 if (fdecal)
507 {
508 fdecal->SetParameter(0,beta);
509 if (!beta) fdecal->SetParameter(0,1);
510 fdecal->SetParameter(1,ped);
511 }
512 }
513 else if (s == "K") // Read the cross talk probability constants
514 {
515 fInput >> jtrans >> jrec >> c >> b >> dlemin >> dlemax;
516 omx=(IceAOM*)fOmdb->GetObject(jtrans,1);
517 if (!omx)
518 {
519 omx=new IceAOM(om);
520 omx->SetUniqueID(jtrans);
521 fOmdb->EnterObject(jtrans,1,omx);
522 }
523
524 TF1* fx=new TF1(fxtalkp);
525 fx->SetParameter(0,c);
526 if (b)
527 {
528 fx->SetParameter(1,b);
529 }
530 else
531 {
532 fx->SetParameter(1,1);
533 }
534 fx->SetParameter(2,dlemin);
535 fx->SetParameter(3,dlemax);
536 fOmdb->EnterObject(jtrans,jrec+1,fx);
537 }
538 else // Skip this line
539 {
540 fInput.ignore(99999,'\n');
541 }
542 }
543}
544///////////////////////////////////////////////////////////////////////////