RALICE/icepack/iceconvert/IceF2k.cxx

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