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