08-mar-2006 NvE Time offset correction in IceF2k extended to allow also a user define...
[u/mrichter/AliRoot.git] / RALICE / icepack / iceconvert / IceF2k.cxx
CommitLineData
f67e2651 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 IceF2k
20// Conversion of Amanda F2K data into IceEvent physics event structures.
5481c137 21// This class is derived from AliJob providing a task-based processing
22// structure on an event-by-event basis.
23// The main object in the job environment is an IceEvent* pointer.
24// In case the user has provided sub-tasks, these will be executed
25// on an event-by-event basis after the IceEvent structure has been filled
26// with the F2K data and before the final structures are written out.
1c9018c6 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
97a09735 30// to investigate/analyse F2K data using the Ralice/IcePack analysis tools.
31//
32// Note : Sometimes the filtering/reco process which produced the F2K file
33// may have introduced a shift (i.e. offset) in the hit times w.r.t.
34// the actual trigger time. The aim of this is to obtain the hit times
35// centered more or less around zero.
36// In case of real data, this is recorded in the F2K data itself and
37// as such will be taken automatically into account by this IceF2k
38// processor such that all times will be provided again unshifted.
39// In other words, all times will be w.r.t. the actual trigger time
40// as recorded in the trigger data device named "Trigger" in the IceEvent
41// structure.
42// In case of simulated data this shift is not available in the F2K data.
43// The offset denoted in the F2K record is related to the time of the
44// primary interaction to put it well ahead of the detector trigger.
45// This primary interaction time, however, is irrelevant for the
46// reconstruction of the recorded hit patterns.
47// If a user had introduced a shift in producing the MC data,
48// very frequently (but not always) a value of -19000 is used.
49// For the IceF2k processing, the user can manually introduce a
50// time offset in case of MC data via the memberfunction SetMcToffset().
51// This user defined offset value will then be used to correct all
52// the hit times such that they will be provided again unshifted
53// w.r.t. the actual trigger time as recorded in the device named
54// "Trigger" in the IceEvent structure.
55// By default the MC time offset is set to 0 in the constructor
56// of this class.
f67e2651 57//
58// Usage example :
59// ---------------
60//
f5b75967 61// Note : This example creates automatically the ROOT output file, which
62// is the most user friendly way of running the conversion job.
63// In the subdirectory /macros the example macro icef2k.cc provides
64// an example of how to create a ROOT output file yourself and passing
65// this file via a pointer to IceF2k.
66//
f67e2651 67// gSystem->Load("ralice");
68// gSystem->Load("icepack");
69// gSystem->Load("iceconvert");
70//
5481c137 71// IceF2k q("IceF2k","F2K to IcePack data structure conversion");
f67e2651 72//
73// // Limit the number of entries for testing
5481c137 74// q.SetMaxEvents(10);
f67e2651 75//
76// // Print frequency to produce a short summary print every printfreq events
5481c137 77// q.SetPrintFreq(1);
f67e2651 78//
79// // Split level for the output structures
5481c137 80// q.SetSplitLevel(2);
f67e2651 81//
82// // Buffer size for the output structures
5481c137 83// q.SetBufferSize(32000);
84//
88b8d522 85// // The F2K input filename(s)
86// q.AddInputFile("run7825.f2k");
f67e2651 87//
5481c137 88// // Output file for the event structures
f5b75967 89// q.SetOutputFile("events.root");
5481c137 90//
91// ///////////////////////////////////////////////////////////////////
92// // Here the user can specify his/her sub-tasks to be executed
93// // on an event-by-event basis after the IceEvent structure
94// // has been filled and before the data is written out.
95// // Sub-tasks (i.e. a user classes derived from TTask) are entered
96// // as follows :
97// //
98// // MyXtalk task1("task1","Cross talk correction");
99// // MyClean task2("task2","Hit cleaning");
100// // q.Add(&task1);
101// // q.Add(&task2);
102// //
103// // The sub-tasks will be executed in the order as they are entered.
104// ///////////////////////////////////////////////////////////////////
105//
106// // Perform the conversion and execute subtasks (if any)
107// // on an event-by-event basis
108// q.ExecuteJob();
f67e2651 109//
110// // Select various objects to be added to the output file
111//
f5b75967 112// TFile* ofile=q.GetOutputFile();
113//
114// if (ofile)
115// {
116// ofile->cd(); // Switch to the output file directory
1c9018c6 117//
f5b75967 118// AliObjMatrix* omdb=q.GetOMdbase();
119// if (omdb) omdb->Write();
f67e2651 120//
f5b75967 121// AliDevice* fitdefs=q.GetFitdefs();
122// if (fitdefs) fitdefs->Write();
f67e2651 123//
f5b75967 124// TDatabasePDG* pdg=q.GetPDG();
125// if (pdg) pdg->Write();
f67e2651 126//
f5b75967 127// // Flush additional objects to the output file.
128// // The output file is not explicitly closed here
129// // to allow interactive investigation of the data tree
130// // when this macro is run in an interactive ROOT/CINT session.
131// ofile->Write();
132// }
f67e2651 133//
134//--- Author: Nick van Eijndhoven 11-mar-2005 Utrecht University
135//- Modified: NvE $Date$ Utrecht University
136///////////////////////////////////////////////////////////////////////////
137
138#include "IceF2k.h"
139#include "Riostream.h"
140
141ClassImp(IceF2k) // Class implementation to enable ROOT I/O
142
5481c137 143IceF2k::IceF2k(const char* name,const char* title) : AliJob(name,title)
f67e2651 144{
145// Default constructor.
5481c137 146// By default maxevent=-1, split=99, bsize=32000, printfreq=1.
f67e2651 147
5481c137 148 fSplit=99;
149 fBsize=32000;
150 fMaxevt=-1;
151 fPrintfreq=1;
88b8d522 152 fInfiles=0;
5481c137 153 fOutfile=0;
f67e2651 154
155 fPdg=0;
156 fOmdb=0;
157 fFitdefs=0;
69884331 158 fTrigdefs=0;
97a09735 159 fToffset=0;
160 fMctoffset=0;
f67e2651 161}
162///////////////////////////////////////////////////////////////////////////
163IceF2k::~IceF2k()
164{
165// Default destructor.
5481c137 166
88b8d522 167 if (fInfiles)
168 {
169 delete fInfiles;
170 fInfiles=0;
171 }
172
f67e2651 173 if (fPdg)
174 {
175 delete fPdg;
176 fPdg=0;
177 }
178
179 if (fOmdb)
180 {
181 delete fOmdb;
182 fOmdb=0;
183 }
184
185 if (fFitdefs)
186 {
187 delete fFitdefs;
188 fFitdefs=0;
189 }
69884331 190
191 if (fTrigdefs)
192 {
193 delete fTrigdefs;
194 fTrigdefs=0;
195 }
f67e2651 196}
197///////////////////////////////////////////////////////////////////////////
5481c137 198void IceF2k::SetMaxEvents(Int_t n)
199{
200// Set the maximum number of events to be processed.
201// n=-1 implies processing of the complete input file, which is the default
202// initialisation in the constructor.
203 fMaxevt=n;
204}
205///////////////////////////////////////////////////////////////////////////
206void IceF2k::SetPrintFreq(Int_t f)
207{
208// Set the printfrequency to produce info every f events.
209// f=1 is the default initialisation in the constructor.
64c21700 210 if (f>=0) fPrintfreq=f;
5481c137 211}
212///////////////////////////////////////////////////////////////////////////
213void IceF2k::SetSplitLevel(Int_t split)
214{
215// Set the split level for the ROOT data file.
216// split=99 is the default initialisation in the constructor.
217 if (split>=0) fSplit=split;
218}
219///////////////////////////////////////////////////////////////////////////
220void IceF2k::SetBufferSize(Int_t bsize)
221{
222// Set the buffer size for the ROOT data file.
223// bsize=32000 is the default initialisation in the constructor.
224 if (bsize>=0) fBsize=bsize;
225}
226///////////////////////////////////////////////////////////////////////////
97a09735 227void IceF2k::SetMcToffset(Float_t toffset)
228{
229// Set a user defined time offset for Monte Carlo data.
230// A very frequently (but not always) used value is -19000.
231// See the introductory docs of this class for further details.
232 fMctoffset=toffset;
233}
234///////////////////////////////////////////////////////////////////////////
5481c137 235void IceF2k::SetInputFile(TString name)
236{
237// Set the name of the F2K input file.
88b8d522 238// This function has become obsolete but is kept for backward compatibility.
239// The user is advised to use AddInputFile() instead, which allows processing
240// of multiple F2K input files.
241// This function will reset the list of all F2K input files and put the specified
242// filename at the first position.
243// Additional F2K input files can be specified via AddInputFile().
244
245 if (fInfiles) delete fInfiles;
246
247 fInfiles=new TObjArray();
248 fInfiles->SetOwner();
249
250 TObjString* s=new TObjString();
251 s->SetString(name);
252 fInfiles->Add(s);
253}
254///////////////////////////////////////////////////////////////////////////
255void IceF2k::AddInputFile(TString name)
256{
257// Add the name of this F2K input file to the list to be processed.
258
259 if (!fInfiles)
260 {
261 fInfiles=new TObjArray();
262 fInfiles->SetOwner();
263 }
264
265 TObjString* s=new TObjString();
266 s->SetString(name);
267 fInfiles->Add(s);
5481c137 268}
269///////////////////////////////////////////////////////////////////////////
270void IceF2k::SetOutputFile(TFile* ofile)
271{
272// Set the output file for the ROOT data.
f5b75967 273 if (fOutfile) delete fOutfile;
5481c137 274 fOutfile=ofile;
275}
276///////////////////////////////////////////////////////////////////////////
f5b75967 277void IceF2k::SetOutputFile(TString name)
278{
279// Create the output file for the ROOT data.
280 if (fOutfile) delete fOutfile;
281 fOutfile=new TFile(name.Data(),"RECREATE","F2K data in IceEvent structure");
282}
283///////////////////////////////////////////////////////////////////////////
284TFile* IceF2k::GetOutputFile()
285{
286// Provide pointer to the ROOT output file.
287 return fOutfile;
288}
289///////////////////////////////////////////////////////////////////////////
f67e2651 290TDatabasePDG* IceF2k::GetPDG()
291{
292// Provide pointer to the PDG database
293 return fPdg;
294}
295///////////////////////////////////////////////////////////////////////////
296AliObjMatrix* IceF2k::GetOMdbase()
297{
298// Provide pointer to the OM geometry, calib. etc... database
299 return fOmdb;
300}
301///////////////////////////////////////////////////////////////////////////
302AliDevice* IceF2k::GetFitdefs()
303{
304// Provide pointer to the fit definitions
305 return fFitdefs;
306}
307///////////////////////////////////////////////////////////////////////////
69884331 308AliDevice* IceF2k::GetTrigdefs()
309{
310// Provide pointer to the trigger definitions
311 return fTrigdefs;
312}
313///////////////////////////////////////////////////////////////////////////
5481c137 314void IceF2k::Exec(Option_t* opt)
f67e2651 315{
5481c137 316// Job to loop over the specified number of events and convert the
f67e2651 317// F2K data into the IceEvent structure.
5481c137 318// If maxevents<0 (default) all the entries of the input file
f67e2651 319// will be processed.
320// Every "printfreq" events a short event summary will be printed.
321// The default value is printfreq=1.
5481c137 322// The output will be written on a standard output tree named "T".
323//
324// Notes :
325// -------
326// 1) This class is derived from AliJob, allowing a task based processing.
327// After the conversion of an F2K event into an IceEvent structure,
328// the processing of all available sub-tasks (if any) is invoked.
329// This provides an event-by-event (sub)task processing before the
330// final data structures are written out.
331// 2) The main object in this job environment is an IceEvent* pointer.
332
88b8d522 333 if (!fInfiles)
5481c137 334 {
88b8d522 335 cout << " *IceF2k Exec* No data input file(s) specified." << endl;
5481c137 336 return;
337 }
f67e2651 338
88b8d522 339 Int_t ninfiles=fInfiles->GetEntries();
340 if (!ninfiles)
5481c137 341 {
88b8d522 342 cout << " *IceF2k Exec* No data input file(s) specified." << endl;
5481c137 343 return;
344 }
f67e2651 345
1c9018c6 346 TTree* otree=0;
347 if (fOutfile)
5481c137 348 {
1c9018c6 349 otree=new TTree("T","F2K Data converted to IceEvent structures");
350 otree->SetDirectory(fOutfile);
5481c137 351 }
f67e2651 352
5481c137 353 IceEvent* evt=new IceEvent();
f67e2651 354 evt->SetTrackCopy(1);
355 evt->SetDevCopy(1);
356
357 // Branch in the tree for the event structure
1c9018c6 358 if (otree) otree->Branch("IceEvent","IceEvent",&evt,fBsize,fSplit);
f67e2651 359
360 // Create the particle database and extend it with some F2000 specific definitions
361 if (!fPdg) fPdg=new TDatabasePDG();
362 Double_t me=fPdg->GetParticle(11)->Mass();
363 fPdg->AddParticle("brems" ,"brems" ,0,1,0,0,"none",10001001,0,0);
364 fPdg->AddParticle("deltae" ,"deltae" ,me,1,0,-3,"Lepton",10001002,0,0);
365 fPdg->AddParticle("pairprod","pairprod",0,1,0,0,"none",10001003,0,0);
366 fPdg->AddParticle("nucl_int","nucl_Int",0,1,0,0,"none",10001004,0,0);
367 fPdg->AddParticle("mu_pair" ,"mu_pair" ,0,1,0,0,"none",10001005,0,0);
368 fPdg->AddParticle("hadrons" ,"hadrons" ,0,1,0,0,"none",10001006,0,0);
369 fPdg->AddParticle("fiberlaser","fiberlaser",0,1,0,0,"none",10002100,0,0);
370 fPdg->AddParticle("n2laser" ,"n2laser" ,0,1,0,0,"none",10002101,0,0);
371 fPdg->AddParticle("yaglaser" ,"yaglaser" ,0,1,0,0,"none",10002201,0,0);
372 fPdg->AddParticle("z_primary","z_primary",0,1,0,0,"none",10003000,0,0);
373 fPdg->AddParticle("a_primary","a_primary",0,1,0,0,"none",10003500,0,0);
374
5481c137 375 // Initialise the job working environment
376 SetMainObject(evt);
1c9018c6 377 if (fOutfile)
378 {
379 AddObject(fOutfile);
380 AddObject(otree);
381 }
5481c137 382
88b8d522 383 TString inputfile;
384
5481c137 385 cout << " ***" << endl;
386 cout << " *** Start processing of job " << GetName() << " ***" << endl;
387 cout << " ***" << endl;
88b8d522 388 for (Int_t i=0; i<ninfiles; i++)
389 {
390 TObjString* sx=(TObjString*)fInfiles->At(i);
391 if (!sx) continue;
392 inputfile=sx->GetString();
393 cout << " F2K input file : " << inputfile.Data() << endl;
394 }
5481c137 395 cout << " Maximum number of events to be processed : " << fMaxevt << endl;
396 cout << " Print frequency : " << fPrintfreq << endl;
1c9018c6 397 if (fOutfile)
398 {
399 cout << " ROOT output file : " << fOutfile->GetName() << endl;
400 cout << " Output characteristics : splitlevel = " << fSplit << " buffersize = " << fBsize << endl;
401 }
5481c137 402
403 ListEnvironment();
88b8d522 404
5481c137 405 Int_t nevt=0;
88b8d522 406 for (Int_t ifile=0; ifile<ninfiles; ifile++)
f67e2651 407 {
88b8d522 408 TObjString* sx=(TObjString*)fInfiles->At(ifile);
409 if (!sx) continue;
f67e2651 410
88b8d522 411 inputfile=sx->GetString();
412 if (inputfile=="") continue;
f67e2651 413
88b8d522 414 // Open the input file in the default ascii format (autodetection) for reading
415 fInput=rdmc_mcopen(inputfile.Data(),"r",RDMC_DEFAULT_ASCII_F);
f67e2651 416
88b8d522 417 if (!fInput)
418 {
419 cout << " *IceF2k Exec* No input file found with name : " << inputfile.Data() << endl;
420 continue;
421 }
69884331 422
88b8d522 423 // Initialise the event structure
424 rdmc_init_mevt(&fEvent);
f67e2651 425
88b8d522 426 // Read the file header information
427 rdmc_rarr(fInput,&fHeader);
f67e2651 428
88b8d522 429 // Fill the database with geometry, calib. etc... parameters
430 // for all the devices
431 FillOMdbase();
f67e2651 432
88b8d522 433 // Set the fit definitions according to the F2000 header info
434 SetFitdefs();
5481c137 435
88b8d522 436 // Set the trigger definitions according to the F2000 header info
437 SetTrigdefs();
438
439 while (!rdmc_revt(fInput,&fHeader,&fEvent))
64c21700 440 {
88b8d522 441 if (fMaxevt>-1 && nevt>=fMaxevt) break;
442
443 // Reset the complete Event structure
444 evt->Reset();
445
446 evt->SetRunNumber(fEvent.nrun);
447 evt->SetEventNumber(fEvent.enr);
448 evt->SetMJD(fEvent.mjd,fEvent.secs,fEvent.nsecs);
449
97a09735 450 // Take trigger offset into account which might have been
451 // introduced during the filtering process.
452 // For simulated data this will be treated separately in PutMcTracks().
453 fToffset=fEvent.t_offset;
454
88b8d522 455 PutTrigger();
456
457 PutMcTracks();
458
459 PutRecoTracks();
f67e2651 460
88b8d522 461 PutHits();
f67e2651 462
88b8d522 463 // Invoke all available sub-tasks (if any)
464 CleanTasks();
465 ExecuteTasks(opt);
466
467 if (fPrintfreq)
468 {
469 if (!(nevt%fPrintfreq)) evt->HeaderData();
470 }
471
472 // Write the complete structure to the output Tree
473 if (otree) otree->Fill();
474
475 // Update event counter
476 nevt++;
477 }
478 if (fMaxevt>-1 && nevt>=fMaxevt) break;
5481c137 479 }
1c9018c6 480
f5b75967 481 // Flush possible memory resident data to the output file
482 if (fOutfile) fOutfile->Write();
483
1c9018c6 484 // Remove the IceEvent object from the environment
485 // and delete it as well
486 if (evt)
487 {
488 RemoveObject(evt);
489 delete evt;
490 }
f67e2651 491}
492///////////////////////////////////////////////////////////////////////////
493void IceF2k::FillOMdbase()
494{
495// Fill the database with geometry, calib. etc... parameters
496// for all the devices.
497
88b8d522 498 if (fHeader.nch<=0)
499 {
500 if (fOmdb)
501 {
502 delete fOmdb;
503 fOmdb=0;
504 }
505 return;
506 }
f67e2651 507
64c21700 508 Int_t adccal=fHeader.is_calib.adc;
509 Int_t tdccal=fHeader.is_calib.tdc;
510 Int_t totcal=fHeader.is_calib.tot;
64c21700 511
512 TF1 fadccal("fadccal","(x-[1])*[0]");
513 TF1 fadcdecal("fadcdecal","(x/[0])+[1]");
216d1d91 514 fadccal.SetParName(0,"BETA-ADC");
515 fadccal.SetParName(1,"PED-ADC");
516 fadcdecal.SetParName(0,"BETA-ADC");
517 fadcdecal.SetParName(1,"PED-ADC");
518
64c21700 519 TF1 ftdccal("ftdccal","(x*[0])-[1]-([0]-1.)*32767.-[2]/sqrt([3])");
520 TF1 ftdcdecal("ftdcdecal","(x+([0]-1.)*32767.+[1]+[2]/sqrt([3]))/[0]");
216d1d91 521 ftdccal.SetParName(0,"BETA-TDC");
522 ftdccal.SetParName(1,"T0");
523 ftdccal.SetParName(2,"ALPHA-TDC");
524 ftdccal.SetParName(3,"ADC-SLEW");
525 ftdcdecal.SetParName(0,"BETA-TDC");
526 ftdcdecal.SetParName(1,"T0");
527 ftdcdecal.SetParName(2,"ALPHA-TDC");
528 ftdcdecal.SetParName(3,"ADC-SLEW");
529
64c21700 530 TF1 ftotcal("ftotcal","x*[0]");
531 TF1 ftotdecal("ftotdecal","x/[0]");
216d1d91 532 ftotcal.SetParName(0,"BETA-TOT");
533 ftotdecal.SetParName(0,"BETA-TOT");
64c21700 534
f67e2651 535 if (fOmdb)
536 {
537 fOmdb->Reset();
538 }
539 else
540 {
541 fOmdb=new AliObjMatrix();
542 fOmdb->SetNameTitle("OMDBASE","The OM geometry, calib. etc... database");
543 fOmdb->SetOwner();
544 }
545
546 IceAOM* dev=0;
547 Double_t pos[3]={0,0,0};
548 for (Int_t i=0; i<fHeader.nch; i++)
549 {
550 dev=new IceAOM();
551 dev->SetUniqueID(i+1);
64c21700 552
553 dev->SetSlotName("ADC",1);
554 dev->SetSlotName("LE",2);
555 dev->SetSlotName("TOT",3);
556
557 dev->SetSlotName("TYPE",4);
558 dev->SetSlotName("ORIENT",5);
559 dev->SetSlotName("THRESH",6);
560 dev->SetSlotName("SENSIT",7);
561 dev->SetSlotName("BETA-TDC",8);
562 dev->SetSlotName("T0",9);
563 dev->SetSlotName("ALPHA-TDC",10);
564 dev->SetSlotName("PED-ADC",11);
565 dev->SetSlotName("BETA-ADC",12);
566 dev->SetSlotName("KAPPA-ADC",13);
567 dev->SetSlotName("PED-TOT",14);
568 dev->SetSlotName("BETA-TOT",15);
569 dev->SetSlotName("KAPPA-TOT",16);
f67e2651 570
571 pos[0]=fHeader.x[i];
572 pos[1]=fHeader.y[i];
573 pos[2]=fHeader.z[i];
574 dev->SetPosition(pos,"car");
64c21700 575
576 fadccal.SetParameter(0,fHeader.cal[i].beta_a);
577 fadccal.SetParameter(1,fHeader.cal[i].ped);
578 fadcdecal.SetParameter(0,fHeader.cal[i].beta_a);
579 if (!fHeader.cal[i].beta_a) fadcdecal.SetParameter(0,1);
580 fadcdecal.SetParameter(1,fHeader.cal[i].ped);
581
582 ftdccal.SetParameter(0,fHeader.cal[i].beta_t);
583 ftdccal.SetParameter(1,fHeader.cal[i].t_0);
584 ftdccal.SetParameter(2,fHeader.cal[i].alpha_t);
585 ftdccal.SetParameter(3,1.e20);
586 ftdcdecal.SetParameter(0,fHeader.cal[i].beta_t);
587 if (!fHeader.cal[i].beta_t) ftdcdecal.SetParameter(0,1);
588 ftdcdecal.SetParameter(1,fHeader.cal[i].t_0);
589 ftdcdecal.SetParameter(2,fHeader.cal[i].alpha_t);
590 ftdcdecal.SetParameter(3,1.e20);
591
592 ftotcal.SetParameter(0,fHeader.cal[i].beta_tot);
593 ftotdecal.SetParameter(0,fHeader.cal[i].beta_tot);
594 if (!fHeader.cal[i].beta_tot) ftotdecal.SetParameter(0,1);
595
596 if (adccal)
597 {
598 dev->SetDecalFunction(&fadcdecal,1);
599 }
600 else
601 {
602 dev->SetCalFunction(&fadccal,1);
603 }
604
605 if (tdccal)
606 {
607 dev->SetDecalFunction(&ftdcdecal,2);
608 }
609 else
610 {
611 dev->SetCalFunction(&ftdccal,2);
612 }
613
614 if (totcal)
615 {
616 dev->SetDecalFunction(&ftotdecal,3);
617 }
618 else
619 {
620 dev->SetCalFunction(&ftotcal,3);
621 }
622
623 dev->SetSignal(fHeader.type[i],4);
624 dev->SetSignal((Float_t)fHeader.costh[i],5);
625 dev->SetSignal(fHeader.thresh[i],6);
626 dev->SetSignal(fHeader.sensit[i],7);
627 dev->SetSignal(fHeader.cal[i].beta_t,8);
628 dev->SetSignal(fHeader.cal[i].t_0,9);
629 dev->SetSignal(fHeader.cal[i].alpha_t,10);
630 dev->SetSignal(fHeader.cal[i].ped,11);
631 dev->SetSignal(fHeader.cal[i].beta_a,12);
632 dev->SetSignal(fHeader.cal[i].kappa,13);
633 dev->SetSignal(fHeader.cal[i].ped_tot,14);
634 dev->SetSignal(fHeader.cal[i].beta_tot,15);
635 dev->SetSignal(fHeader.cal[i].kappa_tot,16);
636
f67e2651 637 fOmdb->EnterObject(i+1,1,dev);
638 }
639}
640///////////////////////////////////////////////////////////////////////////
641void IceF2k::SetFitdefs()
642{
643// Obtain the names of the variables for each fit procedure from the
644// f2000 header. Each different fit procedure is then stored as a separate
69884331 645// "hit" of an AliDevice object and the various fit variables are stored
646// as separate signal slots of the corresponding "hit".
f67e2651 647// Via the GetFitdefs() memberfunction this AliDevice object can be
648// retrieved and stored in the ROOT output file if wanted.
649// The name of the object is FitDefinitions and the stored data can be
650// inspected via the AliDevice::Data() memberfunction and looks as follows :
651//
652// *AliDevice::Data* Id :0 Name : FitDefinitions
653// Position Vector in car coordinates : 0 0 0
654// Err. in car coordinates : 0 0 0
655// The following 8 hits are registered :
656// *AliSignal::Data* Id :0
657// Position Vector in car coordinates : 0 0 0
658// Err. in car coordinates : 0 0 0
659// Owned by device : AliDevice Name : FitDefinitions
69884331 660// Slot : 1 Signal value : 0 name : id
661// Slot : 2 Signal value : 0 name : rchi2
662// Slot : 3 Signal value : 0 name : prob
663// Slot : 4 Signal value : 0 name : sigth
664// Slot : 5 Signal value : 0 name : covmin
665// Slot : 6 Signal value : 0 name : covmax
666// Slot : 7 Signal value : 0 name : cutflag
667// Slot : 8 Signal value : 0 name : chi2
f67e2651 668// *AliSignal::Data* Id :1
669// Position Vector in car coordinates : 0 0 0
670// Err. in car coordinates : 0 0 0
671// Owned by device : AliDevice Name : FitDefinitions
69884331 672// Slot : 1 Signal value : 0 name : id
673// Slot : 2 Signal value : 0 name : rchi2
674// Slot : 3 Signal value : 0 name : prob
f67e2651 675// etc....
676//
677// This memberfunction is based on the original idea/code by Adam Bouchta.
678
88b8d522 679 if (fHeader.n_fit<=0)
680 {
681 if (fFitdefs)
682 {
683 delete fFitdefs;
684 fFitdefs=0;
685 }
686 return;
687 }
f67e2651 688
689 if (fFitdefs)
690 {
691 fFitdefs->Reset(1);
692 }
693 else
694 {
695 fFitdefs=new AliDevice();
696 }
697
698 fFitdefs->SetName("FitDefinitions");
699 fFitdefs->SetHitCopy (1);
700
701 AliSignal s;
702 s.Reset();
703
704 for (Int_t i=0; i<fHeader.n_fit; i++)
705 {
706 s.SetUniqueID(fHeader.def_fit[i].id);
69884331 707 s.SetName(TString(fHeader.def_fit[i].tag));
f67e2651 708
709 for (Int_t j=0; j<fHeader.def_fit[i].nwords; j++)
710 {
711 s.SetSlotName(TString(fHeader.def_fit[i].words[j]),j+1);
69884331 712 s.SetSignal(0,j+1);
f67e2651 713 }
714
715 fFitdefs->AddHit(s);
716 s.Reset(1);
717 }
718}
719///////////////////////////////////////////////////////////////////////////
69884331 720void IceF2k::SetTrigdefs()
721{
722// Obtain the names of the variables for each trigger procedure from the
723// f2000 header. Each different trigger procedure is then stored as a separate
724// "hit" of an AliDevice object and the various trigger variables are stored
725// as separate signal slots of the corresponding "hit".
726// Via the GetFitdefs() memberfunction this AliDevice object can be
727// retrieved and stored in the ROOT output file if wanted.
728// The name of the object is TrigDefinitions and the stored data can be
729// inspected via the AliDevice::Data() memberfunction and looks as follows :
730//
731// *AliDevice::Data* Id : 0 Name : TrigDefinitions
732// Position Vector in car (rad) coordinates : 0 0 0
733// Err. in car (rad) coordinates : 0 0 0
734// The following 9 hits are registered :
735// *AliSignal::Data* Id : 1 Name : main
736// Position Vector in car (rad) coordinates : 0 0 0
737// Err. in car (rad) coordinates : 0 0 0
738// Owned by device : AliDevice Id : 0 Name : TrigDefinitions
739// Slot : 1 Signal value : 0 name : trig_pulse_le
740// Slot : 2 Signal value : 0 name : trig_pulse_tot
741// Slot : 3 Signal value : 0 name : regi_flag
742// *AliSignal::Data* Id : 2 Name : amaa
743// Position Vector in car (rad) coordinates : 0 0 0
744// Err. in car (rad) coordinates : 0 0 0
745// Owned by device : AliDevice Id : 0 Name : TrigDefinitions
746// Slot : 1 Signal value : 0 name : trig_pulse_le
747// Slot : 2 Signal value : 0 name : trig_pulse_tot
748// Slot : 3 Signal value : 0 name : regi_flag
749// *AliSignal::Data* Id : 3 Name : amab10
750// Position Vector in car (rad) coordinates : 0 0 0
751// Err. in car (rad) coordinates : 0 0 0
752// Owned by device : AliDevice Id : 0 Name : TrigDefinitions
753// Slot : 1 Signal value : 0 name : trig_pulse_le
754// Slot : 2 Signal value : 0 name : trig_pulse_tot
755// Slot : 3 Signal value : 0 name : regi_flag
756// etc....
757
88b8d522 758 if (fHeader.n_trigger<=0)
759 {
760 if (fTrigdefs)
761 {
762 delete fTrigdefs;
763 fTrigdefs=0;
764 }
765 return;
766 }
69884331 767
768 if (fTrigdefs)
769 {
770 fTrigdefs->Reset(1);
771 }
772 else
773 {
774 fTrigdefs=new AliDevice();
775 }
776
777 fTrigdefs->SetName("TrigDefinitions");
778 fTrigdefs->SetHitCopy (1);
779
780 AliSignal s;
781 s.Reset();
782
783 for (Int_t i=0; i<fHeader.n_trigger; i++)
784 {
785 s.SetUniqueID(fHeader.def_trig[i].id);
786 s.SetName(TString(fHeader.def_trig[i].tag));
787
788 for (Int_t j=0; j<fHeader.def_trig[i].nwords; j++)
789 {
790 s.SetSlotName(TString(fHeader.def_trig[i].words[j]),j+1);
791 s.SetSignal(0,j+1);
792 }
793
794 fTrigdefs->AddHit(s);
795 s.Reset(1);
796 }
797}
798///////////////////////////////////////////////////////////////////////////
5481c137 799void IceF2k::PutMcTracks()
f67e2651 800{
801// Get the MC tracks from the F2000 file into the IcePack structure.
802// Note : MC tracks are given negative track id's in the event structure.
803// This memberfunction is based on the original code by Adam Bouchta.
804
5481c137 805 IceEvent* evt=(IceEvent*)GetMainObject();
f67e2651 806 if (!evt || fEvent.ntrack<=0) return;
807
97a09735 808 // User defined trigger offset in case of simulated data.
809 // The offset in the F2K file is meant to put the primary interaction
810 // well ahead of the detector trigger.
811 // See the introductory docs of this IceF2k class for further details.
812 fToffset=fMctoffset;
813
f67e2651 814 // Loop over all the tracks and add them to the current event
815 AliTrack t;
816 Double_t vec[3];
817 AliPosition r;
818 Ali3Vector p;
819 Int_t tid=0;
820 Int_t idpdg=0;
821 Int_t idf2k=0;
822 for (Int_t i=0; i<fEvent.ntrack; i++)
823 {
824 t.Reset ();
825
826 // Beginpoint of the track
827 vec[0]=fEvent.gen[i].x;
828 vec[1]=fEvent.gen[i].y;
829 vec[2]=fEvent.gen[i].z;
830 r.SetPosition(vec,"car");
831 t.SetBeginPoint(r);
832
833 // Endpoint of the track
834 vec[0]+=fEvent.gen[i].length*fEvent.gen[i].px;
835 vec[1]+=fEvent.gen[i].length*fEvent.gen[i].py;
836 vec[2]+=fEvent.gen[i].length*fEvent.gen[i].pz;
837 r.SetPosition(vec,"car");
838 t.SetEndPoint(r);
839
840 // Momentum in GeV/c
841 vec[0]=fEvent.gen[i].e*fEvent.gen[i].px*1e-3;
842 vec[1]=fEvent.gen[i].e*fEvent.gen[i].py*1e-3;
843 vec[2]=fEvent.gen[i].e*fEvent.gen[i].pz*1e-3;
844 p.SetVector (vec,"car");
845 t.Set3Momentum(p);
846
847 // MC tracks are indicated by negative track id's
848 tid=fEvent.gen[i].tag;
849 t.SetId(-abs(tid));
850
851 idf2k=fEvent.gen[i].id;
852 idpdg=0;
853 if (idf2k>1000)
854 {
855 idpdg=idf2k+10000000;
856 }
857 else if (idf2k <= 48)
858 {
859 idpdg=fPdg->ConvertGeant3ToPdg(idf2k);
860 }
861 else
862 {
863 if (idf2k==201) idpdg=12;
864 if (idf2k==202) idpdg=14;
865 if (idf2k==203) idpdg=16;
866 if (idf2k==204) idpdg=-12;
867 if (idf2k==205) idpdg=-14;
868 if (idf2k==206) idpdg=-16;
869 }
870
871 t.SetParticleCode(idpdg);
872 t.SetName(fPdg->GetParticle(idpdg)->GetName());
873 t.SetTitle("MC track");
874 t.SetMass(fPdg->GetParticle(idpdg)->Mass());
875 t.SetCharge(fPdg->GetParticle(idpdg)->Charge()/3.);
876
877 evt->AddTrack(t);
878 }
879
880 // Create the pointers to each particle's parent particle.
881 Int_t txid=0;
882 Int_t parid=0;
883 for (Int_t itk=1; itk<=evt->GetNtracks (); itk++)
884 {
885 AliTrack* tx=evt->GetTrack(itk);
886
887 if (!tx) continue;
888
889 txid=tx->GetId();
890
891 parid=-1;
892 for (Int_t j=0; j<fEvent.ntrack; j++)
893 {
894 tid=fEvent.gen[j].tag;
895 if (-abs(tid) == txid) parid=fEvent.gen[j].parent;
896 }
897
898 if (parid<0) continue;
899
900 AliTrack* tpar=evt->GetIdTrack(-abs(parid));
901
902 if (!tpar) continue;
903
904 tx->SetParentTrack(tpar);
905 }
906}
907///////////////////////////////////////////////////////////////////////////
5481c137 908void IceF2k::PutRecoTracks()
f67e2651 909{
910// Get the reconstructed tracks from the F2000 file into the IcePack structure.
911// Note : Reco tracks are given positive track id's in the event structure.
912// This memberfunction is based on the original code by Adam Bouchta.
913
5481c137 914 IceEvent* evt=(IceEvent*)GetMainObject();
f67e2651 915 if (!evt || fEvent.nfit<=0) return;
916
917 // Loop over all the tracks and add them to the current event
918 AliTrack t;
919 Double_t vec[3];
920 AliPosition r;
921 Ali3Vector p;
922 Int_t tid=0;
923 Int_t idpdg=0;
924 Int_t idf2k=0;
925 for (Int_t i=0; i<fEvent.nfit; i++)
926 {
927 t.Reset ();
928
929 // Beginpoint of the track
930 vec[0]=fEvent.rec[i].x;
931 vec[1]=fEvent.rec[i].y;
932 vec[2]=fEvent.rec[i].z;
933 r.SetPosition(vec,"car");
934 t.SetBeginPoint(r);
935
936 // Endpoint of the track
937 vec[0]+=fEvent.rec[i].length*fEvent.rec[i].px;
938 vec[1]+=fEvent.rec[i].length*fEvent.rec[i].py;
939 vec[2]+=fEvent.rec[i].length*fEvent.rec[i].pz;
940 r.SetPosition(vec,"car");
941 t.SetEndPoint(r);
942
943 // Momentum in GeV/c
944 if (fEvent.rec[i].e > 0)
945 {
946 vec[0]=fEvent.rec[i].e*fEvent.rec[i].px*1e-3;
947 vec[1]=fEvent.rec[i].e*fEvent.rec[i].py*1e-3;
948 vec[2]=fEvent.rec[i].e*fEvent.rec[i].pz*1e-3;
949 }
950 else // Give the track a nominal momentum of 1 GeV/c
951 {
952 vec[0]=fEvent.rec[i].px;
953 vec[1]=fEvent.rec[i].py;
954 vec[2]=fEvent.rec[i].pz;
955 }
956 p.SetVector (vec,"car");
957 t.Set3Momentum(p);
958
959 // Use the fit number as track id
960 tid=fEvent.rec[i].tag;
961 t.SetId(abs(tid));
962
963 idf2k=fEvent.rec[i].id;
964 idpdg=0;
965 if (idf2k>1000)
966 {
967 idpdg=idf2k+10000000;
968 }
969 else if (idf2k <= 48)
970 {
971 idpdg=fPdg->ConvertGeant3ToPdg(idf2k);
972 }
973 else
974 {
975 if (idf2k==201) idpdg=12;
976 if (idf2k==202) idpdg=14;
977 if (idf2k==203) idpdg=16;
978 if (idf2k==204) idpdg=-12;
979 if (idf2k==205) idpdg=-14;
980 if (idf2k==206) idpdg=-16;
981 }
982
983 t.SetParticleCode(idpdg);
c29371c1 984 t.SetNameTitle("Sieglinde","RECO track");
f67e2651 985 t.SetMass(fPdg->GetParticle(idpdg)->Mass());
986 t.SetCharge(fPdg->GetParticle(idpdg)->Charge()/3.);
987
988 // Retrieve the various fit parameters for this track
989 AliSignal* fitdata=fFitdefs->GetIdHit(i);
990 for (Int_t jval=0; jval<fEvent.fresult[i].nval; jval++)
991 {
992 fitdata->SetSignal(fEvent.fresult[i].val[jval],jval+1);
993 }
994
995 // Store the various fit parameters for this track
996 t.SetFitDetails(fitdata);
997
998 // Store the various reco tracks as track hypotheses.
999 // A copy of the first reco track is entered as a new track instance
1000 // into the event and all reco tracks (incl. the first one) are
1001 // stored as hypotheses linked to this new reco track.
1002 if (i==0)
1003 {
1004 evt->AddTrack(t);
1005 AliTrack* tx=evt->GetTrack(evt->GetNtracks());
1006 Int_t nrec=evt->GetNtracks(1);
1007 tx->SetId(nrec+1);
1008 }
1009 AliTrack* tx=evt->GetTrack(evt->GetNtracks());
1010 if (tx) tx->AddTrackHypothesis(t);
1011 }
1012}
1013///////////////////////////////////////////////////////////////////////////
5481c137 1014void IceF2k::PutHits()
f67e2651 1015{
1016// Get the hit and waveform info from the F2000 file into the IcePack structure.
1017// This memberfunction is based on the original code by Adam Bouchta.
1018
5481c137 1019 IceEvent* evt=(IceEvent*)GetMainObject();
f67e2651 1020 if (!evt) return;
1021
1022 // Loop over all the hits and add them to the current event
1023 IceAOM om;
1024 AliSignal s;
1025 s.SetSlotName("ADC",1);
1026 s.SetSlotName("LE",2);
1027 s.SetSlotName("TOT",3);
1028 Int_t chan=0;
1029 Int_t maxchan=800;
1030 if (fOmdb) maxchan=fHeader.nch;
1031 IceAOM* omx=0;
1032 AliSignal* sx=0;
1033 Int_t tid=0;
1034 AliTrack* tx=0;
64c21700 1035 Float_t adc=0;
f67e2651 1036 for (Int_t i=0; i<fEvent.nhits; i++)
1037 {
1038 chan=fEvent.h[i].ch+1;
1039 if (chan>maxchan) continue; // Channels 9001, 9002 etc are trigger channels
1040
1041 // Get corresponding device from the current event structure
1042 omx=(IceAOM*)evt->GetIdDevice(chan);
1043 if (!omx)
1044 {
1045 if (fOmdb)
1046 {
1047 omx=(IceAOM*)fOmdb->GetObject(chan,1);
1048 evt->AddDevice(omx);
1049 }
1050 else
1051 {
1052 om.Reset(1);
1053 om.SetUniqueID(chan);
1054 evt->AddDevice(om);
1055 }
1056 omx=(IceAOM*)evt->GetIdDevice(chan);
1057 }
1058
1059 if (!omx) continue;
1060
1061 s.Reset();
1062 s.SetUniqueID(fEvent.h[i].id);
1063 s.SetSignal(fEvent.h[i].amp,1);
97a09735 1064 s.SetSignal((fEvent.h[i].t-fToffset),2);
f67e2651 1065 s.SetSignal(fEvent.h[i].tot,3);
1066
1067 omx->AddHit(s);
1068
1069 sx=omx->GetHit(omx->GetNhits());
1070 if (!sx) continue;
1071
64c21700 1072 // ADC dependent TDC (de)calibration function for this hit
1073 TF1* fcal=omx->GetCalFunction("LE");
1074 TF1* fdecal=omx->GetDecalFunction("LE");
1075 if (fcal) sx->SetCalFunction(fcal,2);
1076 if (fdecal) sx->SetDecalFunction(fdecal,2);
1077 fcal=sx->GetCalFunction(2);
1078 fdecal=sx->GetDecalFunction(2);
1079 adc=sx->GetSignal(1,-4);
1080 if (adc>0)
1081 {
1082 if (fcal) fcal->SetParameter(3,adc);
1083 if (fdecal) fdecal->SetParameter(3,adc);
1084 }
1085 else
1086 {
a4b77ddf 1087 if (fcal) fcal->SetParameter(3,1.e20);
1088 if (fdecal) fdecal->SetParameter(3,1.e20);
64c21700 1089 }
1090
f67e2651 1091 // Bi-directional link between this hit and the track that caused the ADC value.
1092 // This F2K info is probably only present for MC tracks.
1093 tid=fEvent.h[i].ma;
1094 if (tid > 0)
1095 {
1096 tx=evt->GetIdTrack(tid); // Reco tracks
1097 if (!tx) tx=evt->GetIdTrack(-tid); // MC tracks
1098 if (tx) sx->AddLink(tx);
1099 }
1100 else
1101 {
1102 if (tid == -2) s.SetNameTitle("N","Noise");
1103 if (tid == -3) s.SetNameTitle("A","Afterpulse");
1104 }
1105 }
1106
1107 // Loop over all the waveforms and add the histo(s) to the corresponding OM's
1108 TH1F histo;
1109 Int_t nbins=0;
1110 Float_t xlow=0;
1111 Float_t xup=0;
1112 TString hname;
1113 for (Int_t iwf=0; iwf<fEvent.nwf; iwf++)
1114 {
1115 chan=fEvent.wf[iwf].om;
1116 if (chan<=0 || chan>maxchan) continue; // Skip trigger channels
1117
1118 // Get corresponding device from the current event structure
1119 omx=(IceAOM*)evt->GetIdDevice(chan);
1120 if (!omx)
1121 {
1122 if (fOmdb)
1123 {
1124 omx=(IceAOM*)fOmdb->GetObject(chan,1);
1125 evt->AddDevice(omx);
1126 }
1127 else
1128 {
1129 om.Reset(1);
1130 om.SetUniqueID(chan);
1131 evt->AddDevice(om);
1132 }
1133 omx=(IceAOM*)evt->GetIdDevice(chan);
1134 }
1135
1136 if (!omx) continue;
1137
1138 omx->SetSlotName("BASELINE",omx->GetNnames()+1);
1139 omx->SetSignal(fEvent.wf[iwf].baseline,"BASELINE");
1140
1141 // Fill the waveform histogram
1142 hname="OM";
1143 hname+=chan;
1144 hname+="-WF";
1145 hname+=omx->GetNwaveforms()+1;
1146
1147 histo.Reset();
1148 histo.SetName(hname.Data());
1149 nbins=fEvent.wf[iwf].ndigi;
1150 xlow=fEvent.wf[iwf].t_start;
1151 xup=xlow+float(nbins)*fEvent.wf[iwf].t_bin;
1152 histo.SetBins(nbins,xlow,xup);
1153
1154 for (Int_t jbin=1; jbin<=fEvent.wf[iwf].ndigi; jbin++)
1155 {
0e50ad4d 1156 histo.SetBinContent(jbin,fEvent.wf[iwf].digi[jbin-1]);
f67e2651 1157 }
1158
1159 omx->SetWaveform(&histo,omx->GetNwaveforms()+1);
1160 }
1161
1162 // Set bi-directional links between hits and reco track hypotheses.
1163 // Note : Reco tracks are recognised by their positive id.
1164 Int_t hid=0;
1165 TObjArray* rectracks=evt->GetTracks(1);
1166 for (Int_t jtk=0; jtk<rectracks->GetEntries(); jtk++)
1167 {
1168 tx=(AliTrack*)rectracks->At(jtk);
1169 if (!tx) continue;
1170
1171 for (Int_t jhyp=1; jhyp<=tx->GetNhypotheses(); jhyp++)
1172 {
1173 AliTrack* hypx=tx->GetTrackHypothesis(jhyp);
1174 if (!hypx) continue;
1175
1176 // Loop over all combinations of F2K fits and used OM hits
1177 for (Int_t k=0; k<fEvent.nfit_uses; k++)
1178 {
1179 if (fEvent.fit_uses[k].useid != hypx->GetId()) continue;
1180 hid=fEvent.fit_uses[k].hitid;
1181 sx=evt->GetIdHit(hid,"IceAOM");
1182 if (sx) sx->AddLink(hypx);
1183 }
1184 }
1185 }
1186}
1187///////////////////////////////////////////////////////////////////////////
69884331 1188void IceF2k::PutTrigger()
1189{
1190// Get the trigger info from the F2000 file into the IcePack structure.
1191
1192 if (!fTrigdefs) return;
1193
1194 IceEvent* evt=(IceEvent*)GetMainObject();
1195 if (!evt || fEvent.ntrig<=0) return;
1196
1197 AliDevice trig;
1198 trig.SetNameTitle("Trigger","Amanda/IceCube event triggers");
1199 AliSignal s;
1200 TString trigname;
1201 TString slotname;
1202 Int_t id=0;
1203 Int_t nval=0;
1204 for (Int_t i=0; i<fEvent.ntrig; i++)
1205 {
1206 id=fEvent.ptrig[i].id;
1207 nval=fEvent.ptrig[i].nval;
1208 if (!nval) continue;
1209 AliSignal* tdef=fTrigdefs->GetIdHit(id+1);
1210 if (!tdef) continue;
1211 trigname=tdef->GetName();
1212 s.Reset(1);
1213 s.SetName(trigname);
1214 s.SetUniqueID(id+1);
1215 for (Int_t jval=0; jval<fEvent.ptrig[i].nval; jval++)
1216 {
1217 slotname=tdef->GetSlotName(jval+1);
1218 s.SetSlotName(slotname,jval+1);
1219 s.SetSignal(fEvent.ptrig[i].val[jval],jval+1);
1220 }
1221 trig.AddHit(s);
1222 }
1223
1224 // Store the trigger data into the IceEvent structure
1225 evt->AddDevice(trig);
1226}
1227///////////////////////////////////////////////////////////////////////////