RICH becomes HMPID
[u/mrichter/AliRoot.git] / 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     // { "HMPID", "AliHMPIDhit",        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 }