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