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