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