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