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