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