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