Corrected overloading of operator new (icc)
[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.
f67e2651 27//
28// Usage example :
29// ---------------
30//
31// gSystem->Load("ralice");
32// gSystem->Load("icepack");
33// gSystem->Load("iceconvert");
34//
5481c137 35// IceF2k q("IceF2k","F2K to IcePack data structure conversion");
f67e2651 36//
37// // Limit the number of entries for testing
5481c137 38// q.SetMaxEvents(10);
f67e2651 39//
40// // Print frequency to produce a short summary print every printfreq events
5481c137 41// q.SetPrintFreq(1);
f67e2651 42//
43// // Split level for the output structures
5481c137 44// q.SetSplitLevel(2);
f67e2651 45//
46// // Buffer size for the output structures
5481c137 47// q.SetBufferSize(32000);
48//
49// // The F2K input filename
50// q.SetInputFile("run7825.f2k");
f67e2651 51//
5481c137 52// // Output file for the event structures
53// TFile* ofile=new TFile("events.root","RECREATE","F2K data in IceEvent structure");
54// q.SetOutputFile(ofile);
55//
56// ///////////////////////////////////////////////////////////////////
57// // Here the user can specify his/her sub-tasks to be executed
58// // on an event-by-event basis after the IceEvent structure
59// // has been filled and before the data is written out.
60// // Sub-tasks (i.e. a user classes derived from TTask) are entered
61// // as follows :
62// //
63// // MyXtalk task1("task1","Cross talk correction");
64// // MyClean task2("task2","Hit cleaning");
65// // q.Add(&task1);
66// // q.Add(&task2);
67// //
68// // The sub-tasks will be executed in the order as they are entered.
69// ///////////////////////////////////////////////////////////////////
70//
71// // Perform the conversion and execute subtasks (if any)
72// // on an event-by-event basis
73// q.ExecuteJob();
f67e2651 74//
75// // Select various objects to be added to the output file
76//
77// AliObjMatrix* omdb=q.GetOMdbase();
78// if (omdb) omdb->Write();
79//
80// AliDevice* fitdefs=q.GetFitdefs();
81// if (fitdefs) fitdefs->Write();
82//
83// TDatabasePDG* pdg=q.GetPDG();
84// if (pdg) pdg->Write();
85//
86// // Close output file
87// ofile->Write();
88// ofile->Close();
89//
90//--- Author: Nick van Eijndhoven 11-mar-2005 Utrecht University
91//- Modified: NvE $Date$ Utrecht University
92///////////////////////////////////////////////////////////////////////////
93
94#include "IceF2k.h"
95#include "Riostream.h"
96
97ClassImp(IceF2k) // Class implementation to enable ROOT I/O
98
5481c137 99IceF2k::IceF2k(const char* name,const char* title) : AliJob(name,title)
f67e2651 100{
101// Default constructor.
5481c137 102// By default maxevent=-1, split=99, bsize=32000, printfreq=1.
f67e2651 103
5481c137 104 fSplit=99;
105 fBsize=32000;
106 fMaxevt=-1;
107 fPrintfreq=1;
108 fInfile="";
109 fOutfile=0;
f67e2651 110
111 fPdg=0;
112 fOmdb=0;
113 fFitdefs=0;
f67e2651 114}
115///////////////////////////////////////////////////////////////////////////
116IceF2k::~IceF2k()
117{
118// Default destructor.
5481c137 119 IceEvent* evt=(IceEvent*)GetMainObject();
120 if (evt) delete evt;
121
f67e2651 122 if (fPdg)
123 {
124 delete fPdg;
125 fPdg=0;
126 }
127
128 if (fOmdb)
129 {
130 delete fOmdb;
131 fOmdb=0;
132 }
133
134 if (fFitdefs)
135 {
136 delete fFitdefs;
137 fFitdefs=0;
138 }
139}
140///////////////////////////////////////////////////////////////////////////
5481c137 141void IceF2k::SetMaxEvents(Int_t n)
142{
143// Set the maximum number of events to be processed.
144// n=-1 implies processing of the complete input file, which is the default
145// initialisation in the constructor.
146 fMaxevt=n;
147}
148///////////////////////////////////////////////////////////////////////////
149void IceF2k::SetPrintFreq(Int_t f)
150{
151// Set the printfrequency to produce info every f events.
152// f=1 is the default initialisation in the constructor.
153 if (f>0) fPrintfreq=f;
154}
155///////////////////////////////////////////////////////////////////////////
156void IceF2k::SetSplitLevel(Int_t split)
157{
158// Set the split level for the ROOT data file.
159// split=99 is the default initialisation in the constructor.
160 if (split>=0) fSplit=split;
161}
162///////////////////////////////////////////////////////////////////////////
163void IceF2k::SetBufferSize(Int_t bsize)
164{
165// Set the buffer size for the ROOT data file.
166// bsize=32000 is the default initialisation in the constructor.
167 if (bsize>=0) fBsize=bsize;
168}
169///////////////////////////////////////////////////////////////////////////
170void IceF2k::SetInputFile(TString name)
171{
172// Set the name of the F2K input file.
173 fInfile=name;
174}
175///////////////////////////////////////////////////////////////////////////
176void IceF2k::SetOutputFile(TFile* ofile)
177{
178// Set the output file for the ROOT data.
179 fOutfile=ofile;
180}
181///////////////////////////////////////////////////////////////////////////
f67e2651 182TDatabasePDG* IceF2k::GetPDG()
183{
184// Provide pointer to the PDG database
185 return fPdg;
186}
187///////////////////////////////////////////////////////////////////////////
188AliObjMatrix* IceF2k::GetOMdbase()
189{
190// Provide pointer to the OM geometry, calib. etc... database
191 return fOmdb;
192}
193///////////////////////////////////////////////////////////////////////////
194AliDevice* IceF2k::GetFitdefs()
195{
196// Provide pointer to the fit definitions
197 return fFitdefs;
198}
199///////////////////////////////////////////////////////////////////////////
5481c137 200void IceF2k::Exec(Option_t* opt)
f67e2651 201{
5481c137 202// Job to loop over the specified number of events and convert the
f67e2651 203// F2K data into the IceEvent structure.
5481c137 204// If maxevents<0 (default) all the entries of the input file
f67e2651 205// will be processed.
206// Every "printfreq" events a short event summary will be printed.
207// The default value is printfreq=1.
5481c137 208// The output will be written on a standard output tree named "T".
209//
210// Notes :
211// -------
212// 1) This class is derived from AliJob, allowing a task based processing.
213// After the conversion of an F2K event into an IceEvent structure,
214// the processing of all available sub-tasks (if any) is invoked.
215// This provides an event-by-event (sub)task processing before the
216// final data structures are written out.
217// 2) The main object in this job environment is an IceEvent* pointer.
218
219 if (fInfile=="")
220 {
221 cout << " *IceF2k Exec* No data input file specified." << endl;
222 return;
223 }
f67e2651 224
5481c137 225 // Open the input file in the default ascii format (autodetection) for reading
226 fInput=rdmc_mcopen(fInfile.Data(),"r",RDMC_DEFAULT_ASCII_F);
f67e2651 227
5481c137 228 if (!fInput)
229 {
230 cout << " *IceF2k Exec* No input file found with name : " << fInfile.Data() << endl;
231 return;
232 }
f67e2651 233
5481c137 234 // Initialise the event structure
235 rdmc_init_mevt(&fEvent);
236
237 // Read the file header information
238 rdmc_rarr(fInput,&fHeader);
239
240 if (!fOutfile)
241 {
242 cout << " *IceF2k Exec* No ROOT output file specified." << endl;
243 return;
244 }
f67e2651 245
5481c137 246 TTree* otree=new TTree("T","F2K Data converted to IceEvent structures");
247
248 IceEvent* evt=new IceEvent();
f67e2651 249 evt->SetTrackCopy(1);
250 evt->SetDevCopy(1);
251
252 // Branch in the tree for the event structure
253 otree->Branch("IceEvent","IceEvent",&evt,fBsize,fSplit);
254
255 // Create the particle database and extend it with some F2000 specific definitions
256 if (!fPdg) fPdg=new TDatabasePDG();
257 Double_t me=fPdg->GetParticle(11)->Mass();
258 fPdg->AddParticle("brems" ,"brems" ,0,1,0,0,"none",10001001,0,0);
259 fPdg->AddParticle("deltae" ,"deltae" ,me,1,0,-3,"Lepton",10001002,0,0);
260 fPdg->AddParticle("pairprod","pairprod",0,1,0,0,"none",10001003,0,0);
261 fPdg->AddParticle("nucl_int","nucl_Int",0,1,0,0,"none",10001004,0,0);
262 fPdg->AddParticle("mu_pair" ,"mu_pair" ,0,1,0,0,"none",10001005,0,0);
263 fPdg->AddParticle("hadrons" ,"hadrons" ,0,1,0,0,"none",10001006,0,0);
264 fPdg->AddParticle("fiberlaser","fiberlaser",0,1,0,0,"none",10002100,0,0);
265 fPdg->AddParticle("n2laser" ,"n2laser" ,0,1,0,0,"none",10002101,0,0);
266 fPdg->AddParticle("yaglaser" ,"yaglaser" ,0,1,0,0,"none",10002201,0,0);
267 fPdg->AddParticle("z_primary","z_primary",0,1,0,0,"none",10003000,0,0);
268 fPdg->AddParticle("a_primary","a_primary",0,1,0,0,"none",10003500,0,0);
269
270 // Fill the database with geometry, calib. etc... parameters
271 // for all the devices
272 FillOMdbase();
273
274 // Set the fit definitions according to the F2000 header info
275 SetFitdefs();
276
5481c137 277 // Initialise the job working environment
278 SetMainObject(evt);
279 AddObject(fOutfile);
280 AddObject(otree);
281
282 cout << " ***" << endl;
283 cout << " *** Start processing of job " << GetName() << " ***" << endl;
284 cout << " ***" << endl;
285 cout << " F2K input file : " << fInfile.Data() << endl;
286 cout << " Maximum number of events to be processed : " << fMaxevt << endl;
287 cout << " Print frequency : " << fPrintfreq << endl;
288 cout << " ROOT output file : " << fOutfile->GetName() << endl;
289 cout << " Output characteristics : splitlevel = " << fSplit << " buffersize = " << fBsize << endl;
290
291 ListEnvironment();
292
293 Int_t nevt=0;
294 while (!rdmc_revt(fInput,&fHeader,&fEvent))
f67e2651 295 {
5481c137 296 if (fMaxevt>-1 && nevt>=fMaxevt) break;
f67e2651 297
298 // Reset the complete Event structure
299 evt->Reset();
300
301 evt->SetRunNumber(fEvent.nrun);
302 evt->SetEventNumber(fEvent.enr);
303 evt->SetMJD(fEvent.mjd,fEvent.secs,fEvent.nsecs);
304
5481c137 305 PutMcTracks();
f67e2651 306
5481c137 307 PutRecoTracks();
f67e2651 308
5481c137 309 PutHits();
f67e2651 310
5481c137 311 // Invoke all available sub-tasks (if any)
312 ExecuteTasks(opt);
313
314 if (!(nevt%fPrintfreq)) evt->HeaderData();
f67e2651 315
316 // Write the complete structure to the output Tree
317 otree->Fill();
f67e2651 318
5481c137 319 // Update event counter
320 nevt++;
321 }
f67e2651 322}
323///////////////////////////////////////////////////////////////////////////
324void IceF2k::FillOMdbase()
325{
326// Fill the database with geometry, calib. etc... parameters
327// for all the devices.
328
329 if (fHeader.nch<=0) return;
330
331 if (fOmdb)
332 {
333 fOmdb->Reset();
334 }
335 else
336 {
337 fOmdb=new AliObjMatrix();
338 fOmdb->SetNameTitle("OMDBASE","The OM geometry, calib. etc... database");
339 fOmdb->SetOwner();
340 }
341
342 IceAOM* dev=0;
343 Double_t pos[3]={0,0,0};
344 for (Int_t i=0; i<fHeader.nch; i++)
345 {
346 dev=new IceAOM();
347 dev->SetUniqueID(i+1);
348 dev->SetSlotName("TYPE",1);
349 dev->SetSlotName("ORIENT",2);
350 dev->SetSlotName("T0",3);
351 dev->SetSlotName("ALPHA",4);
352 dev->SetSlotName("KADC",5);
353 dev->SetSlotName("KTOT",6);
354 dev->SetSlotName("KTDC",7);
355
356 pos[0]=fHeader.x[i];
357 pos[1]=fHeader.y[i];
358 pos[2]=fHeader.z[i];
359 dev->SetPosition(pos,"car");
360 dev->SetSignal(fHeader.type[i],1);
361 dev->SetSignal((Float_t)fHeader.costh[i],2);
362 dev->SetSignal(fHeader.cal[i].t_0,3);
363 dev->SetSignal(fHeader.cal[i].alpha_t,4);
364 dev->SetSignal(fHeader.cal[i].beta_a,5);
365 dev->SetSignal(fHeader.cal[i].beta_tot,6);
366 dev->SetSignal(fHeader.cal[i].beta_t,7);
367 fOmdb->EnterObject(i+1,1,dev);
368 }
369}
370///////////////////////////////////////////////////////////////////////////
371void IceF2k::SetFitdefs()
372{
373// Obtain the names of the variables for each fit procedure from the
374// f2000 header. Each different fit procedure is then stored as a separate
375// hit of an AliDevice object and the various fit variables are stored
376// as separate signal slots of the corresponding hit.
377// Via the GetFitdefs() memberfunction this AliDevice object can be
378// retrieved and stored in the ROOT output file if wanted.
379// The name of the object is FitDefinitions and the stored data can be
380// inspected via the AliDevice::Data() memberfunction and looks as follows :
381//
382// *AliDevice::Data* Id :0 Name : FitDefinitions
383// Position Vector in car coordinates : 0 0 0
384// Err. in car coordinates : 0 0 0
385// The following 8 hits are registered :
386// *AliSignal::Data* Id :0
387// Position Vector in car coordinates : 0 0 0
388// Err. in car coordinates : 0 0 0
389// Owned by device : AliDevice Name : FitDefinitions
390// Slot : 1 Signal value : 1 name : id
391// Slot : 2 Signal value : 2 name : rchi2
392// Slot : 3 Signal value : 3 name : prob
393// Slot : 4 Signal value : 4 name : sigth
394// Slot : 5 Signal value : 5 name : covmin
395// Slot : 6 Signal value : 6 name : covmax
396// Slot : 7 Signal value : 7 name : cutflag
397// Slot : 8 Signal value : 8 name : chi2
398// *AliSignal::Data* Id :1
399// Position Vector in car coordinates : 0 0 0
400// Err. in car coordinates : 0 0 0
401// Owned by device : AliDevice Name : FitDefinitions
402// Slot : 1 Signal value : 1 name : id
403// Slot : 2 Signal value : 2 name : rchi2
404// Slot : 3 Signal value : 3 name : prob
405// etc....
406//
407// This memberfunction is based on the original idea/code by Adam Bouchta.
408
409 if (fHeader.n_fit<=0) return;
410
411 if (fFitdefs)
412 {
413 fFitdefs->Reset(1);
414 }
415 else
416 {
417 fFitdefs=new AliDevice();
418 }
419
420 fFitdefs->SetName("FitDefinitions");
421 fFitdefs->SetHitCopy (1);
422
423 AliSignal s;
424 s.Reset();
425
426 for (Int_t i=0; i<fHeader.n_fit; i++)
427 {
428 s.SetUniqueID(fHeader.def_fit[i].id);
429
430 for (Int_t j=0; j<fHeader.def_fit[i].nwords; j++)
431 {
432 s.SetSlotName(TString(fHeader.def_fit[i].words[j]),j+1);
433 s.SetSignal(j+1,j+1);
434 }
435
436 fFitdefs->AddHit(s);
437 s.Reset(1);
438 }
439}
440///////////////////////////////////////////////////////////////////////////
5481c137 441void IceF2k::PutMcTracks()
f67e2651 442{
443// Get the MC tracks from the F2000 file into the IcePack structure.
444// Note : MC tracks are given negative track id's in the event structure.
445// This memberfunction is based on the original code by Adam Bouchta.
446
5481c137 447 IceEvent* evt=(IceEvent*)GetMainObject();
f67e2651 448 if (!evt || fEvent.ntrack<=0) return;
449
450 // Loop over all the tracks and add them to the current event
451 AliTrack t;
452 Double_t vec[3];
453 AliPosition r;
454 Ali3Vector p;
455 Int_t tid=0;
456 Int_t idpdg=0;
457 Int_t idf2k=0;
458 for (Int_t i=0; i<fEvent.ntrack; i++)
459 {
460 t.Reset ();
461
462 // Beginpoint of the track
463 vec[0]=fEvent.gen[i].x;
464 vec[1]=fEvent.gen[i].y;
465 vec[2]=fEvent.gen[i].z;
466 r.SetPosition(vec,"car");
467 t.SetBeginPoint(r);
468
469 // Endpoint of the track
470 vec[0]+=fEvent.gen[i].length*fEvent.gen[i].px;
471 vec[1]+=fEvent.gen[i].length*fEvent.gen[i].py;
472 vec[2]+=fEvent.gen[i].length*fEvent.gen[i].pz;
473 r.SetPosition(vec,"car");
474 t.SetEndPoint(r);
475
476 // Momentum in GeV/c
477 vec[0]=fEvent.gen[i].e*fEvent.gen[i].px*1e-3;
478 vec[1]=fEvent.gen[i].e*fEvent.gen[i].py*1e-3;
479 vec[2]=fEvent.gen[i].e*fEvent.gen[i].pz*1e-3;
480 p.SetVector (vec,"car");
481 t.Set3Momentum(p);
482
483 // MC tracks are indicated by negative track id's
484 tid=fEvent.gen[i].tag;
485 t.SetId(-abs(tid));
486
487 idf2k=fEvent.gen[i].id;
488 idpdg=0;
489 if (idf2k>1000)
490 {
491 idpdg=idf2k+10000000;
492 }
493 else if (idf2k <= 48)
494 {
495 idpdg=fPdg->ConvertGeant3ToPdg(idf2k);
496 }
497 else
498 {
499 if (idf2k==201) idpdg=12;
500 if (idf2k==202) idpdg=14;
501 if (idf2k==203) idpdg=16;
502 if (idf2k==204) idpdg=-12;
503 if (idf2k==205) idpdg=-14;
504 if (idf2k==206) idpdg=-16;
505 }
506
507 t.SetParticleCode(idpdg);
508 t.SetName(fPdg->GetParticle(idpdg)->GetName());
509 t.SetTitle("MC track");
510 t.SetMass(fPdg->GetParticle(idpdg)->Mass());
511 t.SetCharge(fPdg->GetParticle(idpdg)->Charge()/3.);
512
513 evt->AddTrack(t);
514 }
515
516 // Create the pointers to each particle's parent particle.
517 Int_t txid=0;
518 Int_t parid=0;
519 for (Int_t itk=1; itk<=evt->GetNtracks (); itk++)
520 {
521 AliTrack* tx=evt->GetTrack(itk);
522
523 if (!tx) continue;
524
525 txid=tx->GetId();
526
527 parid=-1;
528 for (Int_t j=0; j<fEvent.ntrack; j++)
529 {
530 tid=fEvent.gen[j].tag;
531 if (-abs(tid) == txid) parid=fEvent.gen[j].parent;
532 }
533
534 if (parid<0) continue;
535
536 AliTrack* tpar=evt->GetIdTrack(-abs(parid));
537
538 if (!tpar) continue;
539
540 tx->SetParentTrack(tpar);
541 }
542}
543///////////////////////////////////////////////////////////////////////////
5481c137 544void IceF2k::PutRecoTracks()
f67e2651 545{
546// Get the reconstructed tracks from the F2000 file into the IcePack structure.
547// Note : Reco tracks are given positive track id's in the event structure.
548// This memberfunction is based on the original code by Adam Bouchta.
549
5481c137 550 IceEvent* evt=(IceEvent*)GetMainObject();
f67e2651 551 if (!evt || fEvent.nfit<=0) return;
552
553 // Loop over all the tracks and add them to the current event
554 AliTrack t;
555 Double_t vec[3];
556 AliPosition r;
557 Ali3Vector p;
558 Int_t tid=0;
559 Int_t idpdg=0;
560 Int_t idf2k=0;
561 for (Int_t i=0; i<fEvent.nfit; i++)
562 {
563 t.Reset ();
564
565 // Beginpoint of the track
566 vec[0]=fEvent.rec[i].x;
567 vec[1]=fEvent.rec[i].y;
568 vec[2]=fEvent.rec[i].z;
569 r.SetPosition(vec,"car");
570 t.SetBeginPoint(r);
571
572 // Endpoint of the track
573 vec[0]+=fEvent.rec[i].length*fEvent.rec[i].px;
574 vec[1]+=fEvent.rec[i].length*fEvent.rec[i].py;
575 vec[2]+=fEvent.rec[i].length*fEvent.rec[i].pz;
576 r.SetPosition(vec,"car");
577 t.SetEndPoint(r);
578
579 // Momentum in GeV/c
580 if (fEvent.rec[i].e > 0)
581 {
582 vec[0]=fEvent.rec[i].e*fEvent.rec[i].px*1e-3;
583 vec[1]=fEvent.rec[i].e*fEvent.rec[i].py*1e-3;
584 vec[2]=fEvent.rec[i].e*fEvent.rec[i].pz*1e-3;
585 }
586 else // Give the track a nominal momentum of 1 GeV/c
587 {
588 vec[0]=fEvent.rec[i].px;
589 vec[1]=fEvent.rec[i].py;
590 vec[2]=fEvent.rec[i].pz;
591 }
592 p.SetVector (vec,"car");
593 t.Set3Momentum(p);
594
595 // Use the fit number as track id
596 tid=fEvent.rec[i].tag;
597 t.SetId(abs(tid));
598
599 idf2k=fEvent.rec[i].id;
600 idpdg=0;
601 if (idf2k>1000)
602 {
603 idpdg=idf2k+10000000;
604 }
605 else if (idf2k <= 48)
606 {
607 idpdg=fPdg->ConvertGeant3ToPdg(idf2k);
608 }
609 else
610 {
611 if (idf2k==201) idpdg=12;
612 if (idf2k==202) idpdg=14;
613 if (idf2k==203) idpdg=16;
614 if (idf2k==204) idpdg=-12;
615 if (idf2k==205) idpdg=-14;
616 if (idf2k==206) idpdg=-16;
617 }
618
619 t.SetParticleCode(idpdg);
620 t.SetName(fPdg->GetParticle(idpdg)->GetName());
621 t.SetTitle("RECO track");
622 t.SetMass(fPdg->GetParticle(idpdg)->Mass());
623 t.SetCharge(fPdg->GetParticle(idpdg)->Charge()/3.);
624
625 // Retrieve the various fit parameters for this track
626 AliSignal* fitdata=fFitdefs->GetIdHit(i);
627 for (Int_t jval=0; jval<fEvent.fresult[i].nval; jval++)
628 {
629 fitdata->SetSignal(fEvent.fresult[i].val[jval],jval+1);
630 }
631
632 // Store the various fit parameters for this track
633 t.SetFitDetails(fitdata);
634
635 // Store the various reco tracks as track hypotheses.
636 // A copy of the first reco track is entered as a new track instance
637 // into the event and all reco tracks (incl. the first one) are
638 // stored as hypotheses linked to this new reco track.
639 if (i==0)
640 {
641 evt->AddTrack(t);
642 AliTrack* tx=evt->GetTrack(evt->GetNtracks());
643 Int_t nrec=evt->GetNtracks(1);
644 tx->SetId(nrec+1);
645 }
646 AliTrack* tx=evt->GetTrack(evt->GetNtracks());
647 if (tx) tx->AddTrackHypothesis(t);
648 }
649}
650///////////////////////////////////////////////////////////////////////////
5481c137 651void IceF2k::PutHits()
f67e2651 652{
653// Get the hit and waveform info from the F2000 file into the IcePack structure.
654// This memberfunction is based on the original code by Adam Bouchta.
655
5481c137 656 IceEvent* evt=(IceEvent*)GetMainObject();
f67e2651 657 if (!evt) return;
658
659 // Loop over all the hits and add them to the current event
660 IceAOM om;
661 AliSignal s;
662 s.SetSlotName("ADC",1);
663 s.SetSlotName("LE",2);
664 s.SetSlotName("TOT",3);
665 Int_t chan=0;
666 Int_t maxchan=800;
667 if (fOmdb) maxchan=fHeader.nch;
668 IceAOM* omx=0;
669 AliSignal* sx=0;
670 Int_t tid=0;
671 AliTrack* tx=0;
672 for (Int_t i=0; i<fEvent.nhits; i++)
673 {
674 chan=fEvent.h[i].ch+1;
675 if (chan>maxchan) continue; // Channels 9001, 9002 etc are trigger channels
676
677 // Get corresponding device from the current event structure
678 omx=(IceAOM*)evt->GetIdDevice(chan);
679 if (!omx)
680 {
681 if (fOmdb)
682 {
683 omx=(IceAOM*)fOmdb->GetObject(chan,1);
684 evt->AddDevice(omx);
685 }
686 else
687 {
688 om.Reset(1);
689 om.SetUniqueID(chan);
690 evt->AddDevice(om);
691 }
692 omx=(IceAOM*)evt->GetIdDevice(chan);
693 }
694
695 if (!omx) continue;
696
697 s.Reset();
698 s.SetUniqueID(fEvent.h[i].id);
699 s.SetSignal(fEvent.h[i].amp,1);
700 s.SetSignal(fEvent.h[i].t,2);
701 s.SetSignal(fEvent.h[i].tot,3);
702
703 omx->AddHit(s);
704
705 sx=omx->GetHit(omx->GetNhits());
706 if (!sx) continue;
707
708 // Bi-directional link between this hit and the track that caused the ADC value.
709 // This F2K info is probably only present for MC tracks.
710 tid=fEvent.h[i].ma;
711 if (tid > 0)
712 {
713 tx=evt->GetIdTrack(tid); // Reco tracks
714 if (!tx) tx=evt->GetIdTrack(-tid); // MC tracks
715 if (tx) sx->AddLink(tx);
716 }
717 else
718 {
719 if (tid == -2) s.SetNameTitle("N","Noise");
720 if (tid == -3) s.SetNameTitle("A","Afterpulse");
721 }
722 }
723
724 // Loop over all the waveforms and add the histo(s) to the corresponding OM's
725 TH1F histo;
726 Int_t nbins=0;
727 Float_t xlow=0;
728 Float_t xup=0;
729 TString hname;
730 for (Int_t iwf=0; iwf<fEvent.nwf; iwf++)
731 {
732 chan=fEvent.wf[iwf].om;
733 if (chan<=0 || chan>maxchan) continue; // Skip trigger channels
734
735 // Get corresponding device from the current event structure
736 omx=(IceAOM*)evt->GetIdDevice(chan);
737 if (!omx)
738 {
739 if (fOmdb)
740 {
741 omx=(IceAOM*)fOmdb->GetObject(chan,1);
742 evt->AddDevice(omx);
743 }
744 else
745 {
746 om.Reset(1);
747 om.SetUniqueID(chan);
748 evt->AddDevice(om);
749 }
750 omx=(IceAOM*)evt->GetIdDevice(chan);
751 }
752
753 if (!omx) continue;
754
755 omx->SetSlotName("BASELINE",omx->GetNnames()+1);
756 omx->SetSignal(fEvent.wf[iwf].baseline,"BASELINE");
757
758 // Fill the waveform histogram
759 hname="OM";
760 hname+=chan;
761 hname+="-WF";
762 hname+=omx->GetNwaveforms()+1;
763
764 histo.Reset();
765 histo.SetName(hname.Data());
766 nbins=fEvent.wf[iwf].ndigi;
767 xlow=fEvent.wf[iwf].t_start;
768 xup=xlow+float(nbins)*fEvent.wf[iwf].t_bin;
769 histo.SetBins(nbins,xlow,xup);
770
771 for (Int_t jbin=1; jbin<=fEvent.wf[iwf].ndigi; jbin++)
772 {
773 histo.SetBinContent(jbin,fEvent.wf[iwf].digi[jbin]);
774 }
775
776 omx->SetWaveform(&histo,omx->GetNwaveforms()+1);
777 }
778
779 // Set bi-directional links between hits and reco track hypotheses.
780 // Note : Reco tracks are recognised by their positive id.
781 Int_t hid=0;
782 TObjArray* rectracks=evt->GetTracks(1);
783 for (Int_t jtk=0; jtk<rectracks->GetEntries(); jtk++)
784 {
785 tx=(AliTrack*)rectracks->At(jtk);
786 if (!tx) continue;
787
788 for (Int_t jhyp=1; jhyp<=tx->GetNhypotheses(); jhyp++)
789 {
790 AliTrack* hypx=tx->GetTrackHypothesis(jhyp);
791 if (!hypx) continue;
792
793 // Loop over all combinations of F2K fits and used OM hits
794 for (Int_t k=0; k<fEvent.nfit_uses; k++)
795 {
796 if (fEvent.fit_uses[k].useid != hypx->GetId()) continue;
797 hid=fEvent.fit_uses[k].hitid;
798 sx=evt->GetIdHit(hid,"IceAOM");
799 if (sx) sx->AddLink(hypx);
800 }
801 }
802 }
803}
804///////////////////////////////////////////////////////////////////////////