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 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.
  96 // ofile->Write();
  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
 105 ClassImp(IceF2k) // Class implementation to enable ROOT I/O
 106
 107 IceF2k::IceF2k(const char* name,const char* title) : AliJob(name,title)
 108 {
 109 // Default constructor.
 110 // By default maxevent=-1, split=99, bsize=32000, printfreq=1.
 111
 112  fSplit=99;
 113  fBsize=32000;
 114  fMaxevt=-1;
 115  fPrintfreq=1;
 116  fInfile="";
 117  fOutfile=0;
 118
 119  fPdg=0;
 120  fOmdb=0;
 121  fFitdefs=0;
 122 }
 123 ///////////////////////////////////////////////////////////////////////////
 124 IceF2k::~IceF2k()
 125 {
 126 // Default destructor.
 127
 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 ///////////////////////////////////////////////////////////////////////////
 147 void 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 ///////////////////////////////////////////////////////////////////////////
 155 void 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.
 159  if (f>=0) fPrintfreq=f;
 160 }
 161 ///////////////////////////////////////////////////////////////////////////
 162 void 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 ///////////////////////////////////////////////////////////////////////////
 169 void 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 ///////////////////////////////////////////////////////////////////////////
 176 void IceF2k::SetInputFile(TString name)
 177 {
 178 // Set the name of the F2K input file.
 179  fInfile=name;
 180 }
 181 ///////////////////////////////////////////////////////////////////////////
 182 void IceF2k::SetOutputFile(TFile* ofile)
 183 {
 184 // Set the output file for the ROOT data.
 185  fOutfile=ofile;
 186 }
 187 ///////////////////////////////////////////////////////////////////////////
 188 TDatabasePDG* IceF2k::GetPDG()
 189 {
 190 // Provide pointer to the PDG database
 191  return fPdg;
 192 }
 193 ///////////////////////////////////////////////////////////////////////////
 194 AliObjMatrix* IceF2k::GetOMdbase()
 195 {
 196 // Provide pointer to the OM geometry, calib. etc... database
 197  return fOmdb;
 198 }
 199 ///////////////////////////////////////////////////////////////////////////
 200 AliDevice* IceF2k::GetFitdefs()
 201 {
 202 // Provide pointer to the fit definitions
 203  return fFitdefs;
 204 }
 205 ///////////////////////////////////////////////////////////////////////////
 206 void IceF2k::Exec(Option_t* opt)
 207 {
 208 // Job to loop over the specified number of events and convert the
 209 // F2K data into the IceEvent structure.
 210 // If maxevents<0 (default) all the entries of the input file
 211 // will be processed.
 212 // Every "printfreq" events a short event summary will be printed.
 213 // The default value is printfreq=1.
 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  }
 230
 231  // Open the input file in the default ascii format (autodetection) for reading
 232  fInput=rdmc_mcopen(fInfile.Data(),"r",RDMC_DEFAULT_ASCII_F);
 233
 234  if (!fInput)
 235  {
 236   cout << " *IceF2k Exec* No input file found with name : " << fInfile.Data() << endl;
 237   return;
 238  }
 239
 240  // Initialise the event structure
 241  rdmc_init_mevt(&fEvent);
 242
 243  // Read the file header information
 244  rdmc_rarr(fInput,&fHeader);
 245
 246  TTree* otree=0;
 247  if (fOutfile)
 248  {
 249   otree=new TTree("T","F2K Data converted to IceEvent structures");
 250   otree->SetDirectory(fOutfile);
 251  }
 252
 253  IceEvent* evt=new IceEvent();
 254  evt->SetTrackCopy(1);
 255  evt->SetDevCopy(1);
 256
 257  // Branch in the tree for the event structure
 258  if (otree) otree->Branch("IceEvent","IceEvent",&evt,fBsize,fSplit);
 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
 282  // Initialise the job working environment
 283  SetMainObject(evt);
 284  if (fOutfile)
 285  {
 286   AddObject(fOutfile);
 287   AddObject(otree);
 288  }
 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;
 296  if (fOutfile)
 297  {
 298   cout << " ROOT output file : " << fOutfile->GetName() << endl;
 299   cout << " Output characteristics : splitlevel = " << fSplit << " buffersize = " << fBsize << endl;
 300  }
 301
 302  ListEnvironment();
 303
 304  Int_t nevt=0;
 305  while (!rdmc_revt(fInput,&fHeader,&fEvent))
 306  {
 307   if (fMaxevt>-1 && nevt>=fMaxevt) break;
 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
 316   PutMcTracks();
 317
 318   PutRecoTracks();
 319
 320   PutHits();
 321
 322   // Invoke all available sub-tasks (if any)
 323   CleanTasks();
 324   ExecuteTasks(opt);
 325
 326   if (fPrintfreq)
 327   {
 328    if (!(nevt%fPrintfreq)) evt->HeaderData();
 329   }
 330
 331   // Write the complete structure to the output Tree
 332   if (otree) otree->Fill();
 333
 334   // Update event counter
 335   nevt++;
 336  }
 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  }
 345 }
 346 ///////////////////////////////////////////////////////////////////////////
 347 void 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
 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]");
 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
 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]");
 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
 378  TF1 ftotcal("ftotcal","x*[0]");
 379  TF1 ftotdecal("ftotdecal","x/[0]");
 380  ftotcal.SetParName(0,"BETA-TOT");
 381  ftotdecal.SetParName(0,"BETA-TOT");
 382
 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);
 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);
 418
 419   pos[0]=fHeader.x[i];
 420   pos[1]=fHeader.y[i];
 421   pos[2]=fHeader.z[i];
 422   dev->SetPosition(pos,"car");
 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
 485   fOmdb->EnterObject(i+1,1,dev);
 486  }
 487 }
 488 ///////////////////////////////////////////////////////////////////////////
 489 void 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 ///////////////////////////////////////////////////////////////////////////
 559 void IceF2k::PutMcTracks()
 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
 565  IceEvent* evt=(IceEvent*)GetMainObject();
 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 ///////////////////////////////////////////////////////////////////////////
 662 void IceF2k::PutRecoTracks()
 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
 668  IceEvent* evt=(IceEvent*)GetMainObject();
 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 ///////////////////////////////////////////////////////////////////////////
 769 void IceF2k::PutHits()
 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
 774  IceEvent* evt=(IceEvent*)GetMainObject();
 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;
 790  Float_t adc=0;
 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
 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   {
 842    if (fcal) fcal->SetParameter(3,1.e20);
 843    if (fdecal) fdecal->SetParameter(3,1.e20);
 844   }
 845
 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   {
 911    histo.SetBinContent(jbin,fEvent.wf[iwf].digi[jbin-1]);
 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 ///////////////////////////////////////////////////////////////////////////