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