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