EVE/EveBase/AliEveVSDCreator.cxx

   1 // $Id$
   2 // Main authors: Matevz Tadel & Alja Mrak-Tadel: 2006, 2007
   3
   4 /**************************************************************************
   5  * Copyright(c) 1998-2008, ALICE Experiment at CERN, all rights reserved. *
   6  * See http://aliceinfo.cern.ch/Offline/AliRoot/License.html for          *
   7  * full copyright notice.                                                 *
   8  **************************************************************************/
   9
  10 #include <memory>
  11 #include "AliEveVSDCreator.h"
  12
  13 #include "AliEveEventManager.h"
  14
  15 #include <AliStack.h>
  16 #include <AliITSLoader.h>
  17 #include <AliTPCTrackHitsV2.h>
  18 #include <AliPDG.h>
  19 #include <AliHit.h>
  20 #include <AliESDEvent.h>
  21 #include <AliESDv0.h>
  22 #include <AliTPCclusterMI.h>
  23 #include <AliTPCClustersRow.h>
  24 #include <AliITSclusterV2.h>
  25 #include <AliESDkink.h>
  26 #include <AliESDtrack.h>
  27
  28 #include <AliRun.h>
  29 #include <AliRunLoader.h>
  30 #include <AliTPCParam.h>
  31
  32 #include <TSystem.h>
  33 #include <TFile.h>
  34
  35 //______________________________________________________________________________
  36 //
  37 // Create VSD file from ALICE data.
  38
  39 ClassImp(AliEveVSDCreator)
  40
  41 AliEveVSDCreator::AliEveVSDCreator(const Text_t* name, const Text_t* title) :
  42   TEveVSD(name, title),
  43
  44   fTPCHitRes  (2),
  45   fTRDHitRes  (2),
  46
  47   fDebugLevel (0),
  48   fRunLoader  (0),
  49   fGenInfoMap ()
  50 {
  51   // Constructor.
  52
  53   // Particles not in ROOT's PDG database occuring in ALICE
  54   AliPDG::AddParticlesToPdgDataBase();
  55   {
  56     TDatabasePDG *pdgDB = TDatabasePDG::Instance();
  57     // const Int_t kspe=50000000;
  58     const Int_t kion=10000000;
  59
  60     const Double_t kAu2Gev=0.9314943228;
  61     const Double_t khSlash = 1.0545726663e-27;
  62     const Double_t kErg2Gev = 1/1.6021773349e-3;
  63     const Double_t khShGev = khSlash*kErg2Gev;
  64     const Double_t kYear2Sec = 3600*24*365.25;
  65
  66     pdgDB->AddParticle("Deuteron","Deuteron",2*kAu2Gev+8.071e-3,kTRUE,
  67                        0,1,"Ion",kion+10020);
  68     pdgDB->AddParticle("Triton","Triton",3*kAu2Gev+14.931e-3,kFALSE,
  69                        khShGev/(12.33*kYear2Sec),1,"Ion",kion+10030);
  70     pdgDB->AddParticle("Alpha","Alpha",4*kAu2Gev+2.424e-3,kTRUE,
  71                        khShGev/(12.33*kYear2Sec),2,"Ion",kion+20040);
  72     pdgDB->AddParticle("HE3","HE3",3*kAu2Gev+14.931e-3,kFALSE,
  73                        0,2,"Ion",kion+20030);
  74   }
  75
  76   // AliKalmanTrack::SetConvConst(1);
  77 }
  78
  79 /******************************************************************************/
  80
  81 void AliEveVSDCreator::CreateVSD(const Text_t* vsdFile)
  82 {
  83   // Create the VSD for current event in AliEveEventManager.
  84   // Result is stored in vsdFile.
  85   //
  86   // Needs to be extended to support conversion of multiple events.
  87
  88   static const TEveException kEH("AliEveVSDCreator::CreateVSD ");
  89
  90   fRunLoader = AliEveEventManager::AssertRunLoader();
  91
  92   if(fDebugLevel > 0)
  93     printf("%s open seems ok. Now loading sim data.\n", kEH.Data());
  94
  95   fRunLoader->LoadHeader();
  96   fRunLoader->LoadKinematics();
  97   fRunLoader->LoadTrackRefs();
  98   fRunLoader->LoadHits();
  99
 100   // GledNS::PushFD();
 101
 102   if(fDebugLevel > 0)
 103     printf("%s opening output TEveVSD.\n", kEH.Data());
 104
 105   TFile* file = TFile::Open(vsdFile, "RECREATE", "ALICE Visualization Summary Data");
 106   fDirectory = new TDirectoryFile("Event0", "");
 107
 108   if(fDebugLevel > 0)
 109     printf("%s creating trees now ...\n", kEH.Data());
 110
 111   CreateTrees();
 112
 113   if(fDebugLevel > 0)
 114     printf("%s trees created, closing files.\n", kEH.Data());
 115
 116   file->Write();
 117   file->Close();
 118   delete file;
 119   fDirectory =0;
 120
 121   //GledNS::PopFD();
 122
 123   // clean after the TEveVSD data was sucessfuly written
 124   fTreeK      = 0;
 125   fTreeH      = 0;
 126   //fTreeTR     = 0;
 127   fTreeC      = 0;
 128   fTreeV0     = 0;
 129   fTreeKK     = 0;
 130   fTreeR      = 0;
 131   fTreeGI     = 0;
 132
 133   fRunLoader = 0;
 134
 135   if(fDebugLevel > 0)
 136     printf("%s all done.\n", kEH.Data());
 137 }
 138
 139 void AliEveVSDCreator::CreateTrees()
 140 {
 141   // Create and fill the output trees by calling all the
 142   // ConvertXyzz() functions.
 143   // Exceptions from individual functions are displayed as warnings.
 144
 145   static const TEveException kEH("AliEveVSDCreator::CreateTrees ");
 146
 147   if (fDirectory == 0)
 148     throw(kEH + "output directory not set.");
 149
 150   try {
 151     if (fDebugLevel > 1)
 152       printf("%sConvertKinematics.\n", kEH.Data());
 153     ConvertKinematics();
 154   } catch(TEveException& exc) { Warning(kEH, exc); }
 155
 156   Warning(kEH, "Explicitly abandoning further conversion.");
 157   return;
 158
 159   try {
 160     if (fDebugLevel > 1)
 161       printf("%sConvertHits.\n", kEH.Data());
 162     ConvertHits();
 163   } catch(TEveException& exc) { Warning(kEH, exc); }
 164
 165   try {
 166     if (fDebugLevel > 1)
 167       printf("%sConvertClusters.\n", kEH.Data());
 168     ConvertClusters();
 169   } catch(TEveException& exc) { Warning(kEH, exc); }
 170
 171   try {
 172     if (fDebugLevel > 1)
 173       printf("%sConvertRecTracks.\n", kEH.Data());
 174     ConvertRecTracks();
 175   } catch(TEveException& exc) {
 176     Warning(exc, "skipping AliEveV0 extraction.");
 177     goto end_esd_processing;
 178   }
 179
 180   try {
 181     if (fDebugLevel > 1)
 182       printf("%sConvertV0.\n", kEH.Data());
 183     ConvertV0();
 184   } catch(TEveException& exc) { Warning(kEH, exc); }
 185
 186   try {
 187     if (fDebugLevel > 1)
 188       printf("%sConvertKinks.\n", kEH.Data());
 189     ConvertKinks();
 190   } catch(TEveException& exc) { Warning(kEH, exc); }
 191
 192 end_esd_processing:
 193
 194   try {
 195     if (fDebugLevel > 1)
 196       printf("%sConvertGenInfo.\n", kEH.Data());
 197     ConvertGenInfo();
 198   } catch(TEveException& exc) { Warning(kEH, exc); }
 199
 200   return;
 201 }
 202
 203 /******************************************************************************/
 204 // Kinematics
 205 /******************************************************************************/
 206
 207 void AliEveVSDCreator::ConvertKinematics()
 208 {
 209   // Convert kinematics.
 210   // Track references are not stored, they should be.
 211
 212   static const TEveException kEH("AliEveVSDCreator::ConvertKinematics ");
 213
 214   if(fTreeK != 0)
 215     throw (kEH + "kinematics already converted");
 216
 217   AliStack* stack = fRunLoader->Stack();
 218   if(stack == 0)
 219     throw(kEH + "stack is null.");
 220
 221   fDirectory->cd();
 222   fTreeK = new TTree("Kinematics", "TParticles sorted by Label");
 223
 224   Int_t nentries = stack->GetNtrack();
 225   std::vector<TEveMCTrack>  vmc(nentries);
 226   for (Int_t idx=0; idx<nentries; idx++) {
 227     TParticle*   tp = stack->Particle(idx);
 228     vmc[idx]        = *tp;
 229     vmc[idx].fLabel = idx;
 230   }
 231
 232   // read track refrences
 233   // functionality now in AliEveKineTools.
 234   /*
 235   TTree* fTreeTR =  fRunLoader->TreeTR();
 236
 237   if(fTreeTR == 0) {
 238     Warning(kEH, "no TrackRefs; some data will not be available.");
 239   } else {
 240     TClonesArray* RunArrayTR = 0;
 241     fTreeTR->SetBranchAddress("AliRun", &RunArrayTR);
 242
 243     Int_t nPrimaries = (Int_t) fTreeTR->GetEntries();
 244     for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
 245       // printf("T0 fTreeTR->GetEntry(%d) \n",iPrimPart);
 246       fTreeTR->GetEntry(iPrimPart);
 247       // printf("END fTreeTR->GetEntry(%d) \n",iPrimPart);
 248
 249       for (Int_t iTrackRef = 0; iTrackRef < RunArrayTR->GetEntriesFast(); iTrackRef++) {
 250         AliTrackReference *trackRef = (AliTrackReference*)RunArrayTR->At(iTrackRef);
 251         Int_t track = trackRef->GetTrack();
 252         if(track < nentries && track > 0){
 253           TEveMCTrack& mct = vmc[track];
 254           if(trackRef->TestBit(kNotDeleted)) {
 255             mct.decayed   = true;
 256             mct.t_decay   = trackRef->GetTime();
 257             mct.V_decay.x = trackRef->X();
 258             mct.V_decay.y = trackRef->Y();
 259             mct.V_decay.z = trackRef->Z();
 260             mct.P_decay.x = trackRef->Px();
 261             mct.P_decay.y = trackRef->Py();
 262             mct.P_decay.z = trackRef->Pz();
 263             if(TMath::Abs(mct.GetPdgCode()) == 11)
 264               mct.decayed = false; // a bug in TreeTR
 265           }
 266         }
 267       }
 268     }
 269   }
 270   */
 271
 272   fTreeK->Branch("K", "TEveMCTrack",  &fpK, fBuffSize);
 273
 274   printf("sizeofvmc = %d\n", (Int_t) vmc.size());
 275   for(std::vector<TEveMCTrack>::iterator k=vmc.begin(); k!=vmc.end(); ++k) {
 276     TEveMCTrack& mct = *k;
 277     fK = mct;
 278
 279     TParticle* m  = &mct;
 280     Int_t      mi = mct.fLabel;
 281     int cnt = 0;
 282     while(m->GetMother(0) != -1) {
 283       if(cnt > 100) {
 284         printf("cnt %d mi=%d, mo=%d\n", cnt, mi, m->GetMother(0));
 285       }
 286       mi = m->GetMother(0);
 287       m = &vmc[mi];
 288       ++cnt;
 289     }
 290     fK.fEvaLabel = mi;
 291
 292     fTreeK->Fill();
 293   }
 294
 295   fTreeK->BuildIndex("fLabel");
 296 }
 297
 298 /******************************************************************************/
 299 // Hits
 300 /******************************************************************************/
 301
 302 namespace {
 303
 304   struct Detector_t
 305   {
 306     const char*   fName;      // Detector name.
 307     const char*   fHitbranch; // Name of branch containing hits.
 308     unsigned char fDetidx;    // Index identifying the detector internally.
 309   };
 310
 311   Detector_t fgDetectors[] = {
 312     { "ITS",  "AliITShit",         0 },
 313     { "TPC",  "AliTPCTrackHitsV2", 1 },
 314     { "TRD",  "AliTRDhit",         2 },
 315     { "TOF",  "AliTOFhit",         3 }
 316     // { "HMPID", "AliHMPIDhit",        4 },
 317   };
 318
 319 }
 320
 321 /******************************************************************************/
 322
 323 void AliEveVSDCreator::ConvertHits()
 324 {
 325   // Convert MC hits.
 326   // TPC hits are handled specially as they are compressed - only mayor
 327   // hits are stored
 328
 329   static const TEveException kEH("AliEveVSDCreator::ConvertHits ");
 330
 331   if (fTreeH != 0)
 332     throw(kEH + "hits already converted.");
 333
 334   fDirectory->cd();
 335   fTreeH =  new TTree("Hits", "Combined detector hits.");
 336   fTreeH->Branch("H", "TEveHit", &fpH, fBuffSize);
 337
 338   std::map<Int_t, Int_t> hmap;
 339   // parameters for ITS, TPC hits filtering
 340   Float_t x,y,z, x1,y1,z1;
 341   Float_t tpcSqrRes = fTPCHitRes*fTPCHitRes;
 342   Float_t trdSqrRes = fTRDHitRes*fTRDHitRes;
 343
 344   int l=0;
 345   // load hits from the rest of detectors
 346   while (fgDetectors[l].fName != 0)
 347   {
 348     Detector_t& det = fgDetectors[l++];
 349
 350     switch(det.fDetidx)
 351     {
 352       case 1:
 353       {
 354         Int_t count = 0;
 355         TTree* treeh = fRunLoader->GetTreeH(det.fName, false);
 356         if(treeh == 0) {
 357           Warning(kEH, Form("no hits for %s.", det.fName));
 358           continue;
 359         }
 360         AliTPCTrackHitsV2 hv2, *hv2p = &hv2;
 361         treeh->SetBranchAddress("TPC2", &hv2p);
 362         Int_t np = treeh->GetEntries();
 363         for (Int_t i = 0; i < np; ++i)
 364         {
 365           treeh->GetEntry(i);
 366           Int_t evaIdx = np -i -1;
 367           if (hv2.First() == 0) continue;
 368           x = y = z = 0;
 369           do {
 370             AliHit* ah = hv2.GetHit();
 371             x1 = ah->X(); y1 = ah->Y(); z1 = ah->Z();
 372             if ((x-x1)*(x-x1) + (y-y1)*(y-y1) + (z-z1)*(z-z1) > tpcSqrRes)
 373             {
 374               fH.fDetId    = det.fDetidx;
 375               fH.fSubdetId = 0;
 376               fH.fLabel    = ah->Track();
 377               fH.fEvaLabel = evaIdx;
 378               fH.fV.fX = x1; fH.fV.fY = y1; fH.fV.fZ = z1;
 379               fTreeH->Fill();
 380               hmap[fH.fLabel]++;
 381               x = x1; y = y1; z = z1;
 382               count++;
 383             }
 384           } while (hv2.Next());
 385         }
 386         // printf("%d entries in TPChits \n",count);
 387         break;
 388       }
 389       default:
 390       {
 391         TTree* treeh = fRunLoader->GetTreeH(det.fName, false);
 392         if (treeh == 0) {
 393           Warning(kEH, Form("no hits for %s.", det.fName));
 394           continue;
 395         }
 396         TClonesArray *arr = new TClonesArray(det.fHitbranch);
 397         treeh->SetBranchAddress(det.fName, &arr);
 398         Int_t np = treeh->GetEntries();
 399         // in TreeH files hits are grouped in clones arrays
 400         // each eva particle has its own clone array
 401         for (Int_t i = 0; i < np; ++i)
 402         {
 403           treeh->GetEntry(i);
 404           Int_t evaIdx = np -i -1;
 405           Int_t nh=arr->GetEntriesFast();
 406           x = y = z = 0;
 407           // printf("%d entry %d hits for primary %d \n", i, nh, evaIdx);
 408           for (Int_t j = 0; j < nh; ++j)
 409           {
 410             AliHit* aliHit = (AliHit*)arr->UncheckedAt(j);
 411             fH.fDetId    = det.fDetidx;
 412             fH.fSubdetId = 0;
 413             fH.fLabel     = aliHit->GetTrack();
 414             fH.fEvaLabel = evaIdx;
 415             fH.fV.Set(aliHit->X(), aliHit->Y(), aliHit->Z());
 416             if (det.fDetidx == 2)
 417             {
 418               x1=aliHit->X();y1=aliHit->Y();z1=aliHit->Z();
 419               if ((x-x1)*(x-x1)+(y-y1)*(y-y1)+(z-z1)*(z-z1) < trdSqrRes) continue;
 420               x=x1; y=y1; z=z1;
 421             }
 422             hmap[fH.fLabel]++;
 423             fTreeH->Fill();
 424           }
 425         }
 426         delete arr;
 427         break;
 428       } // end default
 429     } // end switch
 430   } // end while
 431
 432
 433   //set geninfo
 434   for(std::map<Int_t, Int_t>::iterator j=hmap.begin(); j!=hmap.end(); ++j)
 435   {
 436     GetGeninfo(j->first)->fNHits += j->second;
 437   }
 438 }
 439
 440 /******************************************************************************/
 441 // Clusters
 442 /******************************************************************************/
 443
 444 void AliEveVSDCreator::ConvertClusters()
 445 {
 446   // Convert clusters.
 447   //
 448   // Only supported for ITS and TPC at the moment, see dedicated
 449   // functions ConvertITSClusters() and ConvertTPCClusters().
 450   //
 451   // It should be possible now to do this in a general manner (with
 452   // the alignment framework).
 453
 454   static const TEveException kEH("AliEveVSDCreator::ConvertClusters ");
 455
 456   if(fTreeC != 0)
 457     throw(kEH + "clusters already converted.");
 458
 459   fDirectory->cd();
 460   fTreeC =  new TTree("Clusters", "rec clusters");
 461   fTreeC->Branch("C", "TEveCluster", &fpC, fBuffSize);
 462
 463   try {
 464     ConvertITSClusters();
 465   } catch(TEveException& exc) { Warning(kEH, exc); }
 466
 467   try {
 468     ConvertTPCClusters();
 469   } catch(TEveException& exc) { Warning(kEH, exc); }
 470 }
 471
 472 /******************************************************************************/
 473
 474 void AliEveVSDCreator::ConvertTPCClusters()
 475 {
 476   // Convert TPC clusters and transform them to global coordinates.
 477
 478   static const TEveException kEH("AliEveVSDCreator::ConvertTPCClusters ");
 479
 480   fRunLoader->LoadRecPoints("TPC");
 481   TTree* tree = fRunLoader->GetTreeR("TPC", false);
 482   if (!tree)
 483     throw(kEH + "'TreeR' not found.");
 484
 485   AliTPCClustersRow  clrow, *clrowp = &clrow;
 486   AliTPCclusterMI   *cl;
 487   clrow.SetClass("AliTPCclusterMI");
 488   tree->SetBranchAddress("Segment", &clrowp);
 489
 490   // count clusters
 491   Int_t nClusters = 0;
 492   Int_t nEnt = tree->GetEntries();
 493   for (Int_t n = 0; n < nEnt; ++n)
 494   {
 495     tree->GetEntry(n);
 496     nClusters += clrow.GetArray()->GetEntriesFast();
 497   }
 498
 499   std::map<Int_t, Int_t> cmap;
 500
 501   for (Int_t n = 0; n < tree->GetEntries(); ++n)
 502   {
 503     tree->GetEntry(n);
 504     Int_t ncl = clrow.GetArray()->GetEntriesFast();
 505     if (ncl > 0)
 506     {
 507       while (ncl--)
 508       {
 509         if (clrow.GetArray())
 510         {
 511           // cl = new AliTPCclusterMI(*(AliTPCclusterMI*)clrow.GetArray()->UncheckedAt(ncl));
 512           cl = (AliTPCclusterMI*)clrow.GetArray()->UncheckedAt(ncl);
 513           if (cl->GetLabel(0) >= 0)
 514           {
 515             Float_t g[3]; //global coordinates
 516             cl->GetGlobalXYZ(g);
 517
 518             fC.fDetId    = 1;
 519             fC.fSubdetId = clrow.GetID();
 520             fC.fLabel[0] = cl->GetLabel(0);
 521             fC.fLabel[1] = cl->GetLabel(1);
 522             fC.fLabel[2] = cl->GetLabel(2);
 523             fC.fV.Set(g);
 524
 525             fTreeC->Fill();
 526             {
 527               int i = 0;
 528               while (i < 3 && fC.fLabel[i])
 529                 cmap[fC.fLabel[i++]]++;
 530             }
 531           }
 532         }
 533       }
 534     }
 535   }
 536   //set geninfo
 537   for (std::map<Int_t, Int_t>::iterator j=cmap.begin(); j!=cmap.end(); ++j)
 538   {
 539     GetGeninfo(j->first)->fNClus += j->second;
 540   }
 541 }
 542
 543 /******************************************************************************/
 544
 545 void AliEveVSDCreator::ConvertITSClusters()
 546 {
 547   // Convert ITS clusters and transform them to global coordinates.
 548
 549   static const TEveException kEH("AliEveVSDCreator::ConvertITSClusters ");
 550
 551   fRunLoader->LoadRecPoints("ITS");
 552   TTree* tree = fRunLoader->GetTreeR("ITS", false);
 553   if (!tree)
 554     throw(kEH + "'TreeR' not found.");
 555
 556   //
 557   AliITSLoader *itsLd = (AliITSLoader*) fRunLoader->GetLoader("ITSLoader");
 558   AliITSgeom   *geom  = itsLd->GetITSgeom();
 559
 560   //printf("alice ITS geom %p \n",geom );
 561
 562   if (!geom)
 563     throw(kEH + "can not find ITS geometry");
 564
 565   TClonesArray *arr = new TClonesArray("AliITSclusterV2");
 566   tree->SetBranchAddress("Clusters", &arr);
 567   Int_t nmods = tree->GetEntries();
 568   Float_t gc[3];
 569   std::map<Int_t, Int_t> cmap;
 570
 571   for (Int_t mod = 0; mod < nmods; ++mod)
 572   {
 573     tree->GetEntry(mod);
 574     Int_t nc=arr->GetEntriesFast();
 575     for (Int_t j = 0; j < nc; ++j)
 576     {
 577       AliITSclusterV2* recp = (AliITSclusterV2*) arr->UncheckedAt(j);
 578
 579       Double_t rot[9];
 580       geom->GetRotMatrix(mod,rot);
 581       Int_t lay,lad,det;
 582       geom->GetModuleId(mod,lay,lad,det);
 583       Float_t tx,ty,tz;
 584       geom->GetTrans(lay,lad,det,tx,ty,tz);
 585
 586       Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
 587       Double_t phi1=TMath::Pi()/2+alpha;
 588       if (lay == 1) phi1+=TMath::Pi();
 589
 590       Float_t cp=TMath::Cos(phi1), sp=TMath::Sin(phi1);
 591       Float_t  r=tx*cp+ty*sp;
 592       gc[0] = r*cp - recp->GetY()*sp;
 593       gc[1] = r*sp + recp->GetY()*cp;
 594       gc[2] = recp->GetZ();
 595
 596       fC.fDetId    = 0;
 597       fC.fSubdetId = 0;
 598       fC.fLabel[0] = recp->GetLabel(0);
 599       fC.fLabel[1] = recp->GetLabel(1);
 600       fC.fLabel[2] = recp->GetLabel(2);
 601       fC.fV.fX     = r*cp - recp->GetY()*sp;
 602       fC.fV.fY     = r*sp + recp->GetY()*cp;
 603       fC.fV.fZ     = recp->GetZ();
 604       fTreeC->Fill();
 605       { int i = 0;
 606         while(i < 3 && fC.fLabel[i])
 607           cmap[fC.fLabel[i++]]++;
 608       }
 609     }
 610
 611     for (std::map<Int_t, Int_t>::iterator j=cmap.begin(); j!=cmap.end(); ++j)
 612     {
 613       GetGeninfo(j->first)->fNClus += j->second;
 614     }
 615   }
 616   delete arr;
 617 }
 618
 619 /******************************************************************************/
 620 // ESD
 621 /******************************************************************************/
 622
 623 void AliEveVSDCreator::ConvertRecTracks()
 624 {
 625   // Convert reconstructed tracks.
 626
 627   static const TEveException kEH("AliEveVSDCreator::ConvertRecTracks ");
 628
 629   if (fTreeR != 0)
 630     throw(kEH + "tracks already converted.");
 631
 632   AliESDEvent* esdEvent = AliEveEventManager::AssertESD();
 633
 634   fDirectory->cd();
 635   fTreeR =  new TTree("RecTracks", "rec tracks");
 636
 637   fTreeR->Branch("R", "TEveRecTrack", &fpR, 512*1024,1);
 638
 639   // reconstructed tracks
 640   AliESDtrack* esdTrack;
 641   Double_t     dbuf[3];
 642   for (Int_t n = 0; n < esdEvent->GetNumberOfTracks(); ++n)
 643   {
 644     esdTrack = esdEvent->GetTrack(n);
 645
 646     fR.fLabel  = esdTrack->GetLabel();
 647     fR.fStatus = (Int_t) esdTrack->GetStatus();
 648     fR.fSign   = (Int_t) esdTrack->GetSign();
 649     esdTrack->GetXYZ(dbuf);    fR.fV.Set(dbuf);
 650     esdTrack->GetPxPyPz(dbuf); fR.fP.Set(dbuf);
 651     Double_t ep = esdTrack->GetP();
 652     fR.fBeta = ep/TMath::Sqrt(ep*ep + TMath::C()*TMath::C()*esdTrack->GetMass()*esdTrack->GetMass());
 653     fTreeR->Fill();
 654   }
 655   fTreeR->BuildIndex("label");
 656   delete esdEvent;
 657 }
 658
 659 /******************************************************************************/
 660
 661 void AliEveVSDCreator::ConvertV0()
 662 {
 663   // Convert reconstructed V0s.
 664
 665   static const TEveException kEH("AliEveVSDCreator::ConvertV0 ");
 666
 667   if(fTreeV0 != 0)
 668     throw(kEH + "AliEveV0 already converted.");
 669
 670   AliESDEvent* esdEvent = AliEveEventManager::AssertESD();
 671
 672   fDirectory->cd();
 673   fTreeV0 =  new TTree("AliEveV0", "AliEveV0 points");
 674
 675   fTreeV0->Branch("AliEveV0", "TEveRecV0", &fpV0, 512*1024,1);
 676
 677   for (Int_t n = 0; n < esdEvent->GetNumberOfV0s(); ++n)
 678   {
 679     AliESDv0    *av     = esdEvent->GetV0(n);
 680     AliESDtrack *trackN = esdEvent->GetTrack(av->GetNindex()); // negative daughter
 681     AliESDtrack *trackP = esdEvent->GetTrack(av->GetPindex()); // positive daughter
 682
 683     Double_t pos[3];
 684
 685     fV0.fStatus = av->GetStatus();
 686     // Point of closest approach
 687     av->GetXYZ(pos[0],pos[1],pos[2]);
 688     fV0.fVCa.fX = pos[0];
 689     fV0.fVCa.fY = pos[1];
 690     fV0.fVCa.fZ = pos[2];
 691     // set birth vertex of neutral particle
 692     av->GetXYZ(pos[0], pos[1], pos[2]);
 693     fV0.fV0Birth.Set(pos);
 694
 695     // momentum and position of negative particle
 696     av->GetParamN()->GetPxPyPz(pos);
 697     fV0.fPNeg.Set(pos);
 698     av->GetParamN()->GetXYZ(pos);
 699     fV0.fVNeg.Set(pos);
 700
 701     // momentum and position of positive particle
 702     av->GetParamP()->GetPxPyPz(pos);
 703     fV0.fPPos.Set(pos);
 704     av->GetParamP()->GetXYZ(pos);
 705     fV0.fVPos.Set(pos);
 706
 707     fV0.fLabel = 0; // !!!! mother label unknown
 708     fV0.fPdg   = av->GetPdgCode();
 709
 710     // daughter indices
 711     fV0.fDLabel[0] = TMath::Abs(trackN->GetLabel());
 712     fV0.fDLabel[1] = TMath::Abs(trackP->GetLabel());
 713
 714     // printf("AliEveV0 convert labels(%d,%d) index(%d,%d)\n",
 715     //     fV0.d_label[0],  fV0.d_label[1],
 716     //     av->GetNIndex(), av->GetPIndex());
 717
 718     fTreeV0->Fill();
 719   }
 720   // if (esdEvent->GetNumberOfV0s()) fTreeV0->BuildIndex("label");
 721   delete esdEvent;
 722 }
 723
 724 /******************************************************************************/
 725
 726 void AliEveVSDCreator::ConvertKinks()
 727 {
 728   // Convert reconstructed kinks.
 729
 730   static const TEveException kEH("AliEveVSDCreator::ConvertKinks ");
 731
 732   if (fTreeKK != 0)
 733     throw(kEH + "Kinks already converted.");
 734
 735   throw kEH + "Currently non-supported - TEveRecKink being updated.";
 736
 737   /*
 738   AliESDEvent* esdEvent = AliEveEventManager::AssertESD();
 739
 740   fDirectory->cd();
 741   fTreeKK =  new TTree("Kinks", "ESD Kinks");
 742
 743   fTreeKK->Branch("KK", "TEveRecKink", &fpKK, fBuffSize);
 744
 745   //  printf("CONVERT KINK Read %d entries in tree kinks \n",  esdEvent->GetNumberOfKinks());
 746   for (Int_t n = 0; n < esdEvent->GetNumberOfKinks(); ++n)
 747   {
 748     AliESDkink* kk = esdEvent->GetKink(n);
 749
 750     Double_t pos[3];
 751
 752     fKK.fLabel  = kk->GetLabel(0);
 753     fKK.fStatus = 0; // status is Char_t[12] ... have no idea how/what to extract.
 754
 755     // reconstructed kink position
 756     fKK.fLabelSec = kk->GetLabel(1);
 757     fKK.fVKink.Set(kk->GetPosition());
 758
 759     const AliExternalTrackParam& tpMother = kk->RefParamMother();
 760     // momentum and position of mother
 761     tpMother.GetPxPyPz(pos);
 762     fKK.fP.Set(pos);
 763     tpMother.GetXYZ(pos);
 764     fKK.fV.Set(pos);
 765     const Double_t* par =  tpMother.GetParameter();
 766     // printf("KINK Pt %f, %f \n",1/tpMother.Pt(),par[4] );
 767     fKK.fSign = (par[4] < 0) ? -1 : 1;
 768
 769     const AliExternalTrackParam& tpDaughter = kk->RefParamDaughter();
 770     // momentum and position of daughter
 771     tpDaughter.GetPxPyPz(pos);
 772     fKK.fPSec.Set(pos);
 773     tpDaughter.GetXYZ(pos);
 774     fKK.fVEnd.Set(pos);
 775
 776     fTreeKK->Fill();
 777   }
 778   if (esdEvent->GetNumberOfKinks()) fTreeKK->BuildIndex("label");
 779   delete esdEvent;
 780   */
 781 }
 782 /******************************************************************************/
 783 // TEveMCRecCrossRef
 784 /******************************************************************************/
 785
 786 void AliEveVSDCreator::ConvertGenInfo()
 787 {
 788   // Build simulation-reconstruction cross-reference table.
 789   // In a rather poor state at the moment.
 790
 791   static const TEveException kEH("AliEveVSDCreator::ConvertGenInfo ");
 792
 793   if(fTreeGI != 0)
 794     throw(kEH + "GI already converted.");
 795
 796   fDirectory->cd();
 797   fTreeGI = new TTree("TEveMCRecCrossRef", "Objects prepared for cross querry");
 798
 799   TEveMCRecCrossRef::Class()->IgnoreTObjectStreamer(true);
 800   fTreeGI->Branch("GI", "TEveMCRecCrossRef",  &fpGI, fBuffSize);
 801   fTreeGI->Branch("K.", "TEveMCTrack",  &fpK);
 802   fTreeGI->Branch("R.", "TEveRecTrack", &fpR);
 803
 804   for (std::map<Int_t, TEveMCRecCrossRef*>::iterator j=fGenInfoMap.begin(); j!=fGenInfoMap.end(); ++j) {
 805     fGI        = *(j->second);
 806     fGI.fLabel = j->first;
 807     fTreeK->GetEntry(j->first);
 808
 809     if (fTreeR) {
 810       Int_t re = fTreeR->GetEntryNumberWithIndex(j->first);
 811       if(re != -1)
 812         fGI.fIsRec = true;
 813     }
 814     //    Int_t hasV0 =  fTreeV0->GetEntryNumberWithIndex(j->first);
 815     //if (hasV0 != -1)
 816     //  fGI.has_AliEveV0 = true;
 817     if (fTreeKK) {
 818       Int_t hasKk =  fTreeKK->GetEntryNumberWithIndex(j->first);
 819       if (hasKk != -1)
 820         fGI.fHasKink = true;
 821     }
 822     fTreeGI->Fill();
 823   }
 824   fGenInfoMap.clear();
 825 }
 826
 827 /******************************************************************************/
 828 /******************************************************************************/
 829 // Protected methods
 830 /******************************************************************************/
 831 /******************************************************************************/
 832
 833 TEveMCRecCrossRef* AliEveVSDCreator::GetGeninfo(Int_t label)
 834 {
 835   // Return the cross-reference structure for given label.
 836   // If the object does not exist it is created.
 837
 838   // printf("get_geninfo %d\n", label);
 839   TEveMCRecCrossRef* gi;
 840   std::map<Int_t, TEveMCRecCrossRef*>::iterator i = fGenInfoMap.find(label);
 841   if (i == fGenInfoMap.end()) {
 842     gi =  new TEveMCRecCrossRef();
 843     fGenInfoMap[label] = gi;
 844   } else {
 845     gi = i->second;
 846   }
 847   return gi;
 848 }