EVE/Alieve/VSDCreator.cxx

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