010fe3133b2c35c29191124c8990f0027dc697fd
[u/mrichter/AliRoot.git] / PWG1 / AliRecInfoMaker.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16
17 ///////////////////////////////////////////////////////////////////////////////
18 //                                                                           //
19 //  Time Projection Chamber                                                  //
20 //  Comparison macro for ESD                                                 //
21 //  responsible: 
22 //  marian.ivanov@cern.ch                                                    //
23 //
24 //
25
26 /* 
27 marian.ivanov@cern.ch 
28 Usage:
29  
30
31
32 gSystem->Load("libPWG1.so");
33 //
34 AliRecInfoMaker *t2 = new AliRecInfoMaker("genTracks.root","cmpESDTracks.root","galice.root",0,0);
35 t2->Exec();
36
37
38 TFile f("cmpESDTracks.root");
39 TTree * tree = (TTree*)f.Get("ESDcmpTracks");
40
41 AliTreeDraw comp;
42 comp.SetTree(tree)
43
44
45
46 //
47 //some cuts definition
48 TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.01&&abs(MC.fVDist[2])<0.01")
49 //TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.5&&abs(MC.fVDist[2])<0.5")
50 //TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<3.9");
51 TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<5");
52 TCut csec("csec","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)>0.5");
53
54
55 TCut crec("crec","fReconstructed==1");
56 TCut cteta1("cteta1","abs(MC.fParticle.Theta()/3.1415-0.5)<0.25");
57 TCut cteta05("cteta05","abs(MC.fParticle.Theta()/3.1415-0.5)<0.1");
58
59 TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
60 TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50");
61 TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)");
62 TCut cchi2("cchi2","fESDtrack.fITSchi2MIP[0]<7.&&fESDtrack.fITSchi2MIP[1]<5.&&fESDtrack.fITSchi2MIP[2]<7.&&fESDtrack.fITSchi2MIP[3]<7.5&&fESDtrack.fITSchi2MIP[4]<6.")
63
64
65 //
66 //example
67 comp.T()->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)");
68 comp.T()->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy");
69 comp.T()->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy");
70 comp.T()->SetAlias("theta","MC.fTrackRef.Theta()");
71 comp.T()->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)");
72 comp.T()->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
73 comp.T()->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
74
75
76 TH1F his("his","his",100,0,20);
77 TH1F hpools("hpools","hpools",100,-7,7);
78 TH1F hfake("hfake","hfake",1000,0,150);
79 TProfile profp0("profp0","profp0",20,-0.4,0.9)
80
81 comp.DrawXY("fTPCinP0[3]","fTPCDelta[4]/fTPCinP1[3]","fReconstructed==1"+cprim,"1",4,0.2,1.5,-0.06,0.06)
82 comp.fRes->Draw();
83 comp.fMean->Draw();  
84
85 comp.DrawXY("fITSinP0[3]","fITSDelta[4]/fITSinP1[3]","fReconstructed==1&&fITSOn"+cprim,"1",4,0.2,1.5,-0.06,0.06)
86 comp.fRes->Draw();
87
88 comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn",20,0.2,1.5)
89 comp.fRes->Draw();
90
91 comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn&&fITSOn&&fESDtrack.fITSFakeRatio<0.1",10,0.2,1.5)
92 comp.fRes->Draw();
93 comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn&&fITSOn&&fESDtrack.fITSFakeRatio>0.1",10,0.2,1.5)
94 comp.fRes->Draw();
95
96 comp.T()->Draw("fESDtrack.fITSsignal/fESDtrack.fTPCsignal","fITSOn&&fTPCOn&&fESDtrack.fITSFakeRatio==0") 
97
98 TH1F his("his","his",100,0,20);
99 TH1F hpools("hpools","hpools",100,-7,7);
100
101 TH2F * hdedx0 = new TH2F("dEdx0","dEdx0",100, 0,2,200,0,550); hdedx0->SetMarkerColor(1);
102 TH2F * hdedx1 = new TH2F("dEdx1","dEdx1",100, 0,2,200,0,550); hdedx1->SetMarkerColor(4);
103 TH2F * hdedx2 = new TH2F("dEdx2","dEdx2",100, 0,2,200,0,550); hdedx2->SetMarkerColor(3);
104 TH2F * hdedx3 = new TH2F("dEdx3","dEdx3",100, 0,2,200,0,550); hdedx3->SetMarkerColor(2);
105
106 comp.T()->Draw("fESDtrack.fITSsignal:MC.fParticle.P()>>dEdx0","fITSOn&&abs(fPdg)==211&&fITStrack.fN==6"+cprim) 
107 comp.T()->Draw("fESDtrack.fITSsignal:MC.fParticle.P()>>dEdx1","fITSOn&&abs(fPdg)==2212&&fITStrack.fN==6"+cprim) 
108 comp.T()->Draw("fESDtrack.fITSsignal:MC.fParticle.P()>>dEdx2","fITSOn&&abs(fPdg)==321&&fITStrack.fN==6"+cprim) 
109 comp.T()->Draw("fESDtrack.fITSsignal:MC.fParticle.P()>>dEdx3","fITSOn&&abs(fPdg)==11&&fITStrack.fN==6"+cprim) 
110
111
112 comp.T()->Draw("fESDtrack.fTRDsignal:MC.fParticle.P()>>dEdx0","fTRDOn&&abs(fPdg)==211&&fTRDtrack.fN>40&&fStatus[2]>1") 
113 comp.T()->Draw("fESDtrack.fTRDsignal:MC.fParticle.P()>>dEdx1","fTRDOn&&abs(fPdg)==2212&&fTRDtrack.fN>40&&fStatus[2]>1") 
114 comp.T()->Draw("fESDtrack.fTRDsignal:MC.fParticle.P()>>dEdx2","fTRDOn&&abs(fPdg)==321&&fTRDtrack.fN>40&&fStatus[2]>1") 
115 comp.T()->Draw("fESDtrack.fTRDsignal:MC.fParticle.P()>>dEdx3","fTRDOn&&abs(fPdg)==11&&fTRDtrack.fN>40&&fStatus[2]>1") 
116
117 comp.T()->Draw("fESDtrack.fTPCsignal:fTPCinP0[4]>>dEdx0","fTPCOn&&abs(fPdg)==211&&fESDtrack.fTPCncls>180&&fESDtrack.fTPCsignal>10"+cteta1); 
118 comp.T()->Draw("fESDtrack.fTPCsignal:fTPCinP0[4]>>dEdx1","fTPCOn&&abs(fPdg)==2212&&fESDtrack.fTPCncls>180&&fESDtrack.fTPCsignal>10"+cteta1); 
119 comp.T()->Draw("fESDtrack.fTPCsignal:fTPCinP0[4]>>dEdx2","fTPCOn&&abs(fPdg)==321&&fESDtrack.fTPCncls>180&&fESDtrack.fTPCsignal>10"+cteta1); 
120 comp.T()->Draw("fESDtrack.fTPCsignal:fTPCinP0[4]>>dEdx3","fTPCOn&&abs(fPdg)==11&&fESDtrack.fTPCncls>180&&fESDtrack.fTPCsignal>10"+cteta1); 
121
122 hdedx3->SetXTitle("P(GeV/c)");
123 hdedx3->SetYTitle("dEdx(unit)");
124 hdedx3->Draw(); hdedx1->Draw("same"); hdedx2->Draw("same"); hdedx0->Draw("same");
125
126 comp.DrawXY("fITSinP0[3]","fITSPools[4]","fReconstructed==1&&fPdg==-211&&fITSOn"+cprim,"1",4,0.2,1.0,-8,8)
127
128 TProfile prof("prof","prof",10,0.5,5);
129
130
131
132
133 */
134
135
136 #include <stdio.h>
137 #include <string.h>
138 //ROOT includes
139 #include "Rtypes.h"
140 #include "TFile.h"
141 #include "TTree.h"
142 #include "TStopwatch.h"
143 #include "TVector3.h"
144 #include "TGeoManager.h"
145 //#include "Getline.h"
146 //
147 //ALIROOT includes
148 //
149 #include "AliRun.h"
150 #include "AliESDtrack.h"
151 #include "AliTPCParam.h"
152 #include "AliTPC.h"
153 #include "AliTrackReference.h"
154 #include "AliTPCParamSR.h"
155 #include "AliTracker.h"
156 #include "AliESDEvent.h"
157 #include "AliESD.h"
158 #include "AliESDfriend.h"
159 #include "AliESDtrack.h"
160 #include "AliTPCseed.h"
161 #include "AliITStrackMI.h"
162 #include "AliESDVertex.h"
163 #include "AliExternalTrackParam.h"
164 #include "AliESDkink.h"
165 #include "AliESDv0.h"
166 #include "AliV0.h"
167 //
168 #include "AliTreeDraw.h"
169 #include "AliMCInfo.h"
170 #include "AliGenKinkInfo.h"
171 #include "AliGenV0Info.h"
172
173
174 #include "AliESDRecInfo.h"
175 #include "AliESDRecV0Info.h"
176 #include "AliESDRecKinkInfo.h"
177 #include "AliRecInfoMaker.h"
178
179
180
181 ClassImp(AliRecInfoMaker)
182
183
184
185
186 AliTPCParam * AliRecInfoMaker::GetTPCParam(){
187   //
188   // create TPC param
189   //
190   AliTPCParamSR * par = new AliTPCParamSR;
191   par->Update();
192   return par;
193 }
194
195
196
197 void  AliRecInfoMaker::MakeAliases(TTree * tree)
198 {
199   //
200   // aliases definition
201   //
202   tree->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)");
203   tree->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy");
204   tree->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy");
205   tree->SetAlias("theta","MC.fTrackRef.Theta()");
206   tree->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)");
207   tree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
208   tree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
209   
210   tree->SetAlias("trddedx","(RC.fESDtrack.fTRDsignals[0]+RC.fESDtrack.fTRDsignals[1]+RC.fESDtrack.fTRDsignals[2]+RC.fESDtrack.fTRDsignals[3]+RC.fESDtrack.fTRDsignals[4]+RC.fESDtrack.fTRDsignals[5])/6.");
211   
212   tree->SetAlias("dtofmc2","fESDtrack.fTrackTime[2]-(10^12*MC.fTOFReferences[0].fTime)");
213   tree->SetAlias("dtofrc2","(fESDtrack.fTrackTime[2]-fESDtrack.fTOFsignal)");
214
215   tree->SetAlias("psum","fESDtrack.fTOFr[4]+fESDtrack.fTOFr[3]+fESDtrack.fTOFr[2]+fESDtrack.fTOFr[1]+fESDtrack.fTOFr[0]");
216   tree->SetAlias("P0","fESDtrack.fTOFr[0]/psum");
217   tree->SetAlias("P1","fESDtrack.fTOFr[1]/psum");
218   tree->SetAlias("P2","fESDtrack.fTOFr[2]/psum");
219   tree->SetAlias("P3","fESDtrack.fTOFr[3]/psum");
220   tree->SetAlias("P4","fESDtrack.fTOFr[4]/psum");
221   tree->SetAlias("MaxP","max(max(max(P0,P1),max(P2,P3)),P4)");
222 }
223
224
225 ////////////////////////////////////////////////////////////////////////
226 AliRecInfoMaker::AliRecInfoMaker(const char* fnGenTracks,
227                                  const char* fnCmp,
228                                  const char* fnGalice,
229                                  Int_t nEvents, Int_t firstEvent):
230   
231   fEventNr(0),                 //! current event number
232   fNEvents(0),                 //! number of events to process
233   fFirstEventNr(0),            //! first event to process
234   fFileCmp(0),                //! output file with cmp tracks
235   fTreeCmp(0),                //! output tree with cmp tracks
236   fTreeCmpKinks(0),                //! output tree with cmp Kinks
237   fTreeCmpV0(0),                //! output tree with cmp V0
238   //
239   fFileGenTracks(0),                //! input files with generated tracks   
240   fTreeGenTracks(0),           //! tree with generated tracks
241   fTreeGenKinks(0),            // tree with gen kinks
242   fTreeGenV0(0),            // tree with gen V0
243   //
244   fLoader(0),         //! pointer to the run loader
245   //
246   fIndexRecTracks(0),         //! index of particle label in the TreeT_ESD
247   fFakeRecTracks(0),          //! number of fake tracks
248   fMultiRecTracks(0),         //! number of multiple reconstructions
249   //
250   fIndexRecKinks(0),         //! index of particle label in treeesd
251   fMultiRecKinks(0),         //! number of multiple reconstructions
252   fSignedKinks(0),           //! indicator that kink was not fake
253   //
254   fIndexRecV0(0),         //! index of particle label in treeesd
255   fMultiRecV0(0),         //! number of multiple reconstructions
256   fSignedV0(0),                //! indicator that kink was not fake
257   //
258   fRecArray(0),           // container with rec infos
259   fEvent(0),             //!event
260   fESDfriend(0),              //!event friend
261   //
262   fParamTPC(0),         //! AliTPCParam
263   fNParticles(0),              //! number of particles in the input tree genTracks
264   fDebug(0),                   //! debug flag  
265   fNextTreeGenEntryToRead(0),    //! last entry already read from genTracks tree
266   fNextKinkToRead(0),            //! last entry already read from genKinks tree
267   fNextV0ToRead(0),            //! last entry already read from genV0 tree
268   //
269   fMCInfo(0),           //! MC information writen per particle
270   fGenKinkInfo(0),      //! MC information writen per Kink
271   fGenV0Info(0),      //! MC information writen per Kink
272   fRecInfo(0),          //! Rec. information writen per particle
273   fFriend(0),          //! friend track
274   fRecKinkInfo(0),    //! reconstructed kink info
275   fRecV0Info(0)    //! reconstructed kink info
276 {
277   // AliRecInfoMaker - connencts the MC information with reconstructed information
278   // fnGenTracks  - file with MC to be created before using AliGenInfoMaker
279   // fnCmp        - file name  to be created  
280   // fnGalice     - file with Loaders - usualy galice.root 
281   //  
282   // nEvent       - number of event s to be analyzed
283   // AliRecInfoMaker *t2 = new AliRecInfoMaker("genTracks.root","cmpESDTracks.root","galice.root",0,0);
284   //
285
286
287   Reset();
288   //  fFnGenTracks = fnGenTracks;
289   //  fFnCmp = fnCmp;
290   sprintf(fFnGenTracks,"%s",fnGenTracks);
291   sprintf(fFnCmp,"%s",fnCmp);
292
293   fFirstEventNr = firstEvent;
294   fEventNr = firstEvent;
295   fNEvents = nEvents;
296   //
297   fLoader = AliRunLoader::Open(fnGalice);
298   if (gAlice){
299     //delete gAlice->GetRunLoader();
300     delete gAlice;
301     gAlice = 0x0;
302   }
303   if (fLoader->LoadgAlice()){
304     cerr<<"Error occured while l"<<endl;
305   }
306   Int_t nall = fLoader->GetNumberOfEvents();
307   if (nEvents==0) {
308     nEvents =nall;
309     fNEvents=nall;
310     fFirstEventNr=0;
311   }    
312
313   if (nall<=0){
314     cerr<<"no events available"<<endl;
315     fEventNr = 0;
316     return;
317   }
318   if (firstEvent+nEvents>nall) {
319     fEventNr = nall-firstEvent;
320     cerr<<"restricted number of events availaible"<<endl;
321   }
322   AliMagF * magf = gAlice->Field();
323   AliTracker::SetFieldMap(magf,0);
324   TGeoManager::Import("geometry.root");
325
326
327 }
328 ////////////////////////////////////////////////////////////////////////
329 AliRecInfoMaker::AliRecInfoMaker(const AliRecInfoMaker& /*info*/):
330   
331   fEventNr(0),                 //! current event number
332   fNEvents(0),                 //! number of events to process
333   fFirstEventNr(0),            //! first event to process
334   fFileCmp(0),                //! output file with cmp tracks
335   fTreeCmp(0),                //! output tree with cmp tracks
336   fTreeCmpKinks(0),                //! output tree with cmp Kinks
337   fTreeCmpV0(0),                //! output tree with cmp V0
338   //
339   fFileGenTracks(0),                //! input files with generated tracks   
340   fTreeGenTracks(0),           //! tree with generated tracks
341   fTreeGenKinks(0),            // tree with gen kinks
342   fTreeGenV0(0),            // tree with gen V0
343   //
344   fLoader(0),         //! pointer to the run loader
345   //
346   fIndexRecTracks(0),         //! index of particle label in the TreeT_ESD
347   fFakeRecTracks(0),          //! number of fake tracks
348   fMultiRecTracks(0),         //! number of multiple reconstructions
349   //
350   fIndexRecKinks(0),         //! index of particle label in treeesd
351   fMultiRecKinks(0),         //! number of multiple reconstructions
352   fSignedKinks(0),           //! indicator that kink was not fake
353   //
354   fIndexRecV0(0),         //! index of particle label in treeesd
355   fMultiRecV0(0),         //! number of multiple reconstructions
356   fSignedV0(0),                //! indicator that kink was not fake
357   //
358   fRecArray(0),           // container with rec infos
359   fEvent(0),             //!event
360   fESDfriend(0),              //!event friend
361   //
362   fParamTPC(0),         //! AliTPCParam
363   fNParticles(0),              //! number of particles in the input tree genTracks
364   fDebug(0),                   //! debug flag  
365   fNextTreeGenEntryToRead(0),    //! last entry already read from genTracks tree
366   fNextKinkToRead(0),            //! last entry already read from genKinks tree
367   fNextV0ToRead(0),            //! last entry already read from genV0 tree
368   //
369   fMCInfo(0),           //! MC information writen per particle
370   fGenKinkInfo(0),      //! MC information writen per Kink
371   fGenV0Info(0),      //! MC information writen per Kink
372   fRecInfo(0),          //! Rec. information writen per particle
373   fFriend(0),          //! friend track
374   fRecKinkInfo(0),    //! reconstructed kink info
375   fRecV0Info(0)    //! reconstructed kink info
376 {
377   //
378   // Dummy copu constructor
379   //
380 }
381
382
383
384
385 ////////////////////////////////////////////////////////////////////////
386 AliRecInfoMaker::~AliRecInfoMaker()
387 {
388   //
389   // Destructor
390   //
391   if (fLoader) {
392     delete fLoader;
393   }
394 }
395
396 //////////////////////////////////////////////////////////////
397 Int_t AliRecInfoMaker::SetIO()
398 {
399   //
400   // SetIO  - Create the input trees
401   //
402   CreateTreeCmp();
403   if (!fTreeCmp) return 1;
404   fParamTPC = GetTPCParam();
405   //
406   if (!ConnectGenTree()) {
407     cerr<<"Cannot connect tree with generated tracks"<<endl;
408     return 1;
409   }
410   return 0;
411 }
412
413 //////////////////////////////////////////////////////////////
414
415 Int_t AliRecInfoMaker::SetIO(Int_t eventNr)
416 {
417   //
418   // 
419   // SET INPUT
420   //
421   TFile f("AliESDs.root");
422   //
423  
424   TTree* tree = (TTree*) f.Get("esdTree");
425   tree->SetBranchStatus("*",1);
426   fEvent = new AliESDEvent;
427   
428   if (tree->GetBranch("ESD")){
429     //    tree->SetBranchAddress("ESD", &fEvent);
430     // tree->SetBranchAddress("ESDfriend.",&fESDfriend);
431     // tree->GetEntry(eventNr);
432     // fEvent->SetESDfriend(fESDfriend);    
433   }else{
434     fEvent->ReadFromTree(tree);
435     fESDfriend = (AliESDfriend*)fEvent->FindListObject("AliESDfriend"); 
436     tree->GetEntry(eventNr);
437     fEvent->SetESDfriend(fESDfriend); 
438   }    
439   
440
441
442   if (!fEvent) return 1;
443
444   return 0;
445 }
446
447
448
449 ////////////////////////////////////////////////////////////////////////
450 void AliRecInfoMaker::Reset()
451 {
452   //
453   // Reset the class
454   //
455   fEventNr = 0;
456   fNEvents = 0;
457   fTreeCmp = 0;
458   fTreeCmpKinks =0;
459   fTreeCmpV0 =0;
460   //  fFnCmp = "cmpTracks.root";
461   fFileGenTracks = 0;
462   fDebug = 0;
463   //
464   fParamTPC = 0;
465   fEvent =0;
466 }
467
468 ////////////////////////////////////////////////////////////////////////
469 Int_t AliRecInfoMaker::Exec(Int_t nEvents, Int_t firstEventNr)
470 {
471   //
472   // Exec comparison for subrange of events
473   //
474   fNEvents = nEvents;
475   fFirstEventNr = firstEventNr;
476   return Exec();
477 }
478
479 ////////////////////////////////////////////////////////////////////////
480 Int_t AliRecInfoMaker::Exec()
481 {
482   //
483   // Exec comparison
484   //
485   TStopwatch timer;
486   timer.Start();
487
488   if (SetIO()==1) 
489     return 1;
490    
491   fNextTreeGenEntryToRead = 0;
492   fNextKinkToRead = 0;
493   fNextV0ToRead   =0;
494   cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
495   for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
496        eventNr++) {
497     fEventNr = eventNr;
498     SetIO(fEventNr);
499     fNParticles = gAlice->GetEvent(fEventNr);    
500
501     fIndexRecTracks = new Short_t[fNParticles*20];  //write at maximum 4 tracks corresponding to particle
502     fIndexRecKinks  = new Short_t[fNParticles*20];  //write at maximum 20 tracks corresponding to particle
503     fIndexRecV0  = new Short_t[fNParticles*20];  //write at maximum 20 tracks corresponding to particle
504
505     fFakeRecTracks = new Short_t[fNParticles];
506     fMultiRecTracks = new Short_t[fNParticles];
507     fMultiRecKinks = new Short_t[fNParticles];
508     fMultiRecV0 = new Short_t[fNParticles];
509
510     for (Int_t i = 0; i<fNParticles; i++) {
511       for (Int_t j=0;j<20;j++){
512         fIndexRecTracks[20*i+j] = -1;
513         fIndexRecKinks[20*i+j]  = -1;
514         fIndexRecV0[20*i+j]  = -1;
515       }
516       fFakeRecTracks[i] = 0;
517       fMultiRecTracks[i] = 0;
518       fMultiRecKinks[i] = 0;
519       fMultiRecV0[i] = 0;      
520     }
521   
522     cout<<"Start to process event "<<fEventNr<<endl;
523     cout<<"\tfNParticles = "<<fNParticles<<endl;
524     if (fDebug>2) cout<<"\tStart loop over TreeT"<<endl;
525     if (TreeTLoop()>0) return 1;
526
527     if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl;
528     if (TreeGenLoop(eventNr)>0) return 1;
529     //BuildKinkInfo0(eventNr);
530     BuildV0Info(eventNr); // no V0 info for a moment
531     fRecArray->Delete();
532
533     if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl;
534
535     delete [] fIndexRecTracks;
536     delete [] fIndexRecKinks;
537     delete [] fIndexRecV0;
538     delete [] fFakeRecTracks;
539     delete [] fMultiRecTracks;
540     delete [] fMultiRecKinks;
541     delete [] fMultiRecV0;
542   }
543
544   CloseOutputFile();
545
546   cerr<<"Exec finished"<<endl;
547   timer.Stop();
548   timer.Print();
549   return 0;
550
551 }
552 ////////////////////////////////////////////////////////////////////////
553 Bool_t AliRecInfoMaker::ConnectGenTree()
554 {
555 //
556 // connect all branches from the genTracksTree
557 // use the same variables as for the new cmp tree, it may work
558 //
559   fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
560   if (!fFileGenTracks) {
561     cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
562     return kFALSE;
563   }
564   fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
565   if (!fTreeGenTracks) {
566     cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
567         <<fFnGenTracks<<endl;
568     return kFALSE;
569   }
570   //
571   fMCInfo = new AliMCInfo;
572   fTreeGenTracks->SetBranchAddress("MC",&fMCInfo);
573   //
574   //
575   fTreeGenKinks = (TTree*)fFileGenTracks->Get("genKinksTree");
576   if (!fTreeGenKinks) {
577     cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
578         <<fFnGenTracks<<endl;
579     //return kFALSE;
580   }
581   else{
582     fGenKinkInfo = new AliGenKinkInfo;
583     fTreeGenKinks->SetBranchAddress("MC",&fGenKinkInfo);
584   }
585
586   fTreeGenV0 = (TTree*)fFileGenTracks->Get("genV0Tree");
587   if (!fTreeGenV0) {
588     cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
589         <<fFnGenTracks<<endl;
590     //return kFALSE;
591   }
592   else{
593     fGenV0Info = new AliGenV0Info;
594     fTreeGenV0->SetBranchAddress("MC",&fGenV0Info);
595   }
596   //
597   if (fDebug > 1) {
598     cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
599   }
600   return kTRUE;
601 }
602
603
604 ////////////////////////////////////////////////////////////////////////
605 void AliRecInfoMaker::CreateTreeCmp() 
606 {
607   //
608   // Create file and tree with comparison information 
609   //
610   fFileCmp = TFile::Open(fFnCmp,"RECREATE");
611   if (!fFileCmp) {
612     cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
613     return;
614   }
615   //
616   //
617   fTreeCmp    = new TTree("ESDcmpTracks","ESDcmpTracks");
618   fMCInfo = new AliMCInfo;
619   fRecInfo = new AliESDRecInfo;
620   AliESDtrack * esdTrack = new AliESDtrack; 
621   //  AliITStrackMI * itsTrack = new AliITStrackMI;  
622   fTreeCmp->Branch("MC","AliMCInfo",&fMCInfo,256000);
623   fTreeCmp->Branch("RC","AliESDRecInfo",&fRecInfo,256000);
624   //  fTreeCmp->Branch("ITS","AliITStrackMI",&itsTrack);
625   delete esdTrack;
626   //
627   //
628   fTreeCmpKinks    = new TTree("ESDcmpKinks","ESDcmpKinks"); 
629   fGenKinkInfo     = new AliGenKinkInfo;
630   fRecKinkInfo     = new AliESDRecKinkInfo;
631   fTreeCmpKinks->Branch("MC.","AliGenKinkInfo",&fGenKinkInfo,256000);
632   fTreeCmpKinks->Branch("RC.","AliESDRecKinkInfo",&fRecKinkInfo,256000);
633   //
634   //
635   fTreeCmpV0       = new TTree("ESDcmpV0","ESDcmpV0"); 
636   fGenV0Info     = new AliGenV0Info;
637   fRecV0Info     = new AliESDRecV0Info;
638   fTreeCmpV0->Branch("MC.","AliGenV0Info",   &fGenV0Info,256000);
639   fTreeCmpV0->Branch("RC.","AliESDRecV0Info",&fRecV0Info,256000);
640   //
641   fTreeCmp->AutoSave(); 
642   fTreeCmpKinks->AutoSave(); 
643   fTreeCmpV0->AutoSave(); 
644 }
645
646 ////////////////////////////////////////////////////////////////////////
647 void AliRecInfoMaker::CloseOutputFile()  
648 {
649   //
650   // Close output file
651   //
652
653   if (!fFileCmp) {
654     cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
655     return;
656   }
657   fFileCmp->cd();
658   fTreeCmp->Write();    
659   delete fTreeCmp;
660   
661   fFileCmp->Close();
662   delete fFileCmp;
663   return;
664 }
665 ////////////////////////////////////////////////////////////////////////
666
667 TVector3 AliRecInfoMaker::TR2Local(AliTrackReference *trackRef,
668                             AliTPCParam *paramTPC) {
669
670   //
671   // Transform position to the local coord frame
672   //
673   
674   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
675   Int_t index[4];
676   paramTPC->Transform0to1(x,index);
677   paramTPC->Transform1to2Ideal(x,index);
678   return TVector3(x);
679 }
680 ////////////////////////////////////////////////////////////////////////
681
682 Int_t AliRecInfoMaker::TreeTLoop()
683 {
684   //
685   // loop over all ESD reconstructed tracks and store info in memory
686   //
687   // + loop over all reconstructed kinks
688   TStopwatch  timer;
689   timer.Start();
690   //  
691   Int_t nEntries = (Int_t)fEvent->GetNumberOfTracks();  
692   Int_t nKinks = (Int_t) fEvent->GetNumberOfKinks();
693   Int_t nV0MIs = (Int_t) fEvent->GetNumberOfV0s();
694   fSignedKinks = new Short_t[nKinks];
695   fSignedV0    = new Short_t[nV0MIs];
696   //
697   // load kinks to the memory
698   for (Int_t i=0; i<nKinks;i++){
699     //    AliESDkink * kink =
700     fEvent->GetKink(i);
701     fSignedKinks[i]=0;
702   }
703   //
704   for (Int_t i=0; i<nV0MIs;i++){
705     //AliV0 * v0MI = 
706     (AliV0*)fEvent->GetV0(i);
707     fSignedV0[i]=0;
708   }
709   
710   //
711   //
712   AliESDtrack * track=0;
713   for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
714     track = (AliESDtrack*)fEvent->GetTrack(iEntry);
715     //
716     Int_t label = track->GetLabel();
717     Int_t absLabel = abs(label);
718     if (absLabel < fNParticles) {
719       //      fIndexRecTracks[absLabel] =  iEntry;
720       if (label < 0) fFakeRecTracks[absLabel]++;      
721       if (fMultiRecTracks[absLabel]>0){
722         if (fMultiRecTracks[absLabel]<20)
723           fIndexRecTracks[absLabel*20+fMultiRecTracks[absLabel]] =  iEntry;     
724       }
725       else      
726         fIndexRecTracks[absLabel*20] =  iEntry;
727       fMultiRecTracks[absLabel]++;
728     }
729   }
730   // sort reconstructed kinks  
731   //
732   AliESDkink * kink=0;
733   for (Int_t iEntry=0; iEntry<nKinks;iEntry++){
734     kink = (AliESDkink*)fEvent->GetKink(iEntry);
735     if (!kink) continue;
736     //
737     Int_t label0 = TMath::Abs(kink->GetLabel(0));
738     Int_t label1 = TMath::Abs(kink->GetLabel(1));
739     Int_t absLabel = TMath::Min(label0,label1);
740     if (absLabel < fNParticles) {
741       if (fMultiRecKinks[absLabel]>0){
742         if (fMultiRecKinks[absLabel]<20)
743           fIndexRecKinks[absLabel*20+fMultiRecKinks[absLabel]] =  iEntry;       
744       }
745       else      
746         fIndexRecKinks[absLabel*20] =  iEntry;
747       fMultiRecKinks[absLabel]++;
748     }
749   }  
750   // --sort reconstructed V0
751   //
752   AliV0 * v0MI=0;
753   for (Int_t iEntry=0; iEntry<nV0MIs;iEntry++){
754     v0MI = (AliV0*)fEvent->GetV0(iEntry);
755     if (!v0MI) continue;
756     //
757     //
758     //
759     //Int_t label0 = TMath::Abs(v0MI->GetLabel(0));
760     //Int_t label1 = TMath::Abs(v0MI->GetLabel(1));
761     AliESDtrack * trackn = fEvent->GetTrack((v0MI->GetNindex()));
762     AliESDtrack * trackp = fEvent->GetTrack((v0MI->GetPindex()));
763     Int_t labels[2]={-1,-1};
764     labels[0] = (trackn==0) ? -1 : trackn->GetLabel(); 
765     labels[1] = (trackp==0) ? -1 : trackp->GetLabel(); 
766     //
767     for (Int_t i=0;i<2;i++){
768       Int_t absLabel =  labels[i];
769       if (absLabel < fNParticles) {
770         if (fMultiRecV0[absLabel]>0){
771           if (fMultiRecV0[absLabel]<20)
772             fIndexRecV0[absLabel*20+fMultiRecV0[absLabel]] =  iEntry;   
773         }
774         else      
775           fIndexRecV0[absLabel*20] =  iEntry;
776         fMultiRecV0[absLabel]++;
777       }
778     }
779   }  
780
781
782   printf("Time spended in TreeTLoop\n");
783   timer.Print();
784   
785   if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;  
786   return 0;
787 }
788
789 ////////////////////////////////////////////////////////////////////////
790 Int_t AliRecInfoMaker::TreeGenLoop(Int_t eventNr)
791 {
792 //
793 // loop over all entries for a given event, find corresponding 
794 // rec. track and store in the fTreeCmp
795 //
796   TStopwatch timer;
797   timer.Start();
798   Int_t entry = fNextTreeGenEntryToRead;
799   Double_t nParticlesTR = fTreeGenTracks->GetEntriesFast();
800   cerr<<"fNParticles, nParticlesTR, fNextTreeGenEntryToRead: "<<fNParticles<<" "
801       <<nParticlesTR<<" "<<fNextTreeGenEntryToRead<<endl;
802   TBranch * branch = fTreeCmp->GetBranch("RC");
803   TBranch * branchF = fTreeCmp->GetBranch("F");
804   
805   branch->SetAddress(&fRecInfo); // set all pointers
806   fRecArray = new TObjArray(fNParticles);
807   AliESDtrack dummytrack;  //
808   AliESDfriendTrack dummytrackF;  //
809
810   while (entry < nParticlesTR) {
811     fTreeGenTracks->GetEntry(entry);
812     entry++;
813     if (eventNr < fMCInfo->fEventNr) continue;
814     if (eventNr > fMCInfo->fEventNr) continue;
815     if (fMCInfo->GetCharge()==0) continue;
816     //
817     fNextTreeGenEntryToRead = entry-1;
818     if (fDebug > 2 && fMCInfo->fLabel < 10) {
819       cerr<<"Fill track with a label "<<fMCInfo->fLabel<<endl;
820     }
821     //    if (fMCInfo->fNTPCRef<1) continue; // not TPCref
822     //
823     fRecInfo->Reset();
824     AliESDtrack * track=0;
825     fRecInfo->fReconstructed =0;
826     TVector3 local = TR2Local(&(fMCInfo->fTrackRef),fParamTPC);
827     local.GetXYZ(fRecInfo->fTRLocalCoord);      
828     //
829     if (fIndexRecTracks[fMCInfo->fLabel*20] >= 0) {
830       track= (AliESDtrack*)fEvent->GetTrack(fIndexRecTracks[fMCInfo->fLabel*20]);
831       //
832       //
833       // find nearest track if multifound
834       //Int_t sign = Int_t(track->GetSign()*fMCInfo->fCharge);
835       //
836       Int_t status = 0;
837       if  ((track->GetStatus()&AliESDtrack::kITSrefit)>0) status++;
838       if  ((track->GetStatus()&AliESDtrack::kTPCrefit)>0) status++;
839       if  ((track->GetStatus()&AliESDtrack::kTRDrefit)>0) status++;
840
841       //
842       if (fIndexRecTracks[fMCInfo->fLabel*20+1]>0){
843         //
844         Double_t p[3];
845         track->GetInnerPxPyPz(p);
846         Float_t maxp = p[0]*p[0]+p[1]*p[1]+p[2]*p[2];
847         //
848         for (Int_t i=1;i<20;i++){
849           if (fIndexRecTracks[fMCInfo->fLabel*20+i]>=0){
850             AliESDtrack * track2 = (AliESDtrack*)fEvent->GetTrack(fIndexRecTracks[fMCInfo->fLabel*20+i]);
851             if (!track2) continue;
852             //Int_t sign2 = track->GetSign()*fMCInfo->fCharge; //           
853             //if (sign2<0) continue;
854             track2->GetInnerPxPyPz(p);
855             Float_t mom = p[0]*p[0]+p[1]*p[1]+p[2]*p[2];
856             /*
857             if (sign<0){
858               sign = sign2;
859               track = track2;
860               maxp = mom;
861               continue;
862             }
863             */
864             //
865             Int_t status2 = 0;
866             if  ((track2->GetStatus()&AliESDtrack::kITSrefit)>0) status2++;
867             if  ((track2->GetStatus()&AliESDtrack::kTPCrefit)>0) status2++;
868             if  ((track2->GetStatus()&AliESDtrack::kTRDrefit)>0) status2++;
869             if (status2<status) continue;
870             //
871             if (mom<maxp) continue;
872             maxp = mom;
873             track = track2;
874             //
875           }
876         }
877       } 
878       //
879       if (track) {
880         fRecInfo->SetESDtrack(track);
881       }else{
882         fRecInfo->SetESDtrack(&dummytrack);
883       }
884       //
885
886       fRecInfo->fReconstructed = 1;
887       fRecInfo->fFake     = fFakeRecTracks[fMCInfo->fLabel];
888       fRecInfo->fMultiple = fMultiRecTracks[fMCInfo->fLabel];
889       //
890       fRecInfo->Update(fMCInfo,fParamTPC,kTRUE);          
891     }
892     else{
893       fRecInfo->SetESDtrack(&dummytrack);
894       fRecInfo->Update(fMCInfo,fParamTPC,kFALSE);
895     }
896     fRecArray->AddAt(new AliESDRecInfo(*fRecInfo),fMCInfo->fLabel);
897     fTreeCmp->Fill();
898   }
899   fTreeCmp->AutoSave();
900   printf("Time spended in TreeGenLoop\n");
901   timer.Print();
902   if (fDebug > 2) cerr<<"end of TreeGenLoop"<<endl;
903
904   return 0;
905 }
906
907
908
909 ////////////////////////////////////////////////////////////////////////
910 ////////////////////////////////////////////////////////////////////////
911 ////////////////////////////////////////////////////////////////////////
912 Int_t AliRecInfoMaker::BuildKinkInfo0(Int_t eventNr)
913 {
914 //
915 // loop over all entries for a given event, find corresponding 
916 // rec. track and store in the fTreeCmp
917 //
918   TStopwatch timer;
919   timer.Start();
920   Int_t entry = fNextKinkToRead;
921   Double_t nParticlesTR = fTreeGenKinks->GetEntriesFast();
922   cerr<<"fNParticles, nParticlesTR, fNextKinkToRead: "<<fNParticles<<" "
923       <<nParticlesTR<<" "<<fNextKinkToRead<<endl;
924   //
925   TBranch * branch = fTreeCmpKinks->GetBranch("RC.");
926   branch->SetAddress(&fRecKinkInfo); // set all pointers
927   
928   //
929   while (entry < nParticlesTR) {
930     fTreeGenKinks->GetEntry(entry);
931     entry++;
932     if (eventNr < fGenKinkInfo->GetMinus().fEventNr) continue;
933     if (eventNr > fGenKinkInfo->GetMinus().fEventNr) continue;;
934     //
935     fNextKinkToRead = entry-1;
936     //
937     //
938     AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(fGenKinkInfo->GetMinus().fLabel);
939     AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(fGenKinkInfo->GetPlus().fLabel);
940     fRecKinkInfo->fT1 = (*fRecInfo1);
941     fRecKinkInfo->fT2 = (*fRecInfo2);
942     fRecKinkInfo->fStatus =0;
943     if (fRecInfo1 && fRecInfo1->fTPCOn) fRecKinkInfo->fStatus+=1;
944     if (fRecInfo2 && fRecInfo2->fTPCOn) fRecKinkInfo->fStatus+=2;
945     if (fRecKinkInfo->fStatus==3&&fRecInfo1->fSign!=fRecInfo2->fSign) fRecKinkInfo->fStatus*=-1;
946     
947     if (fRecKinkInfo->fStatus==3){
948       fRecKinkInfo->Update();    
949     }
950     Int_t label =  TMath::Min(fGenKinkInfo->GetMinus().fLabel,fGenKinkInfo->GetPlus().fLabel);
951     Int_t label2 = TMath::Max(fGenKinkInfo->GetMinus().fLabel,fGenKinkInfo->GetPlus().fLabel);
952     
953     AliESDkink *kink=0;
954     fRecKinkInfo->fRecStatus   =0;
955     fRecKinkInfo->fMultiple    = fMultiRecKinks[label];
956     fRecKinkInfo->fKinkMultiple=0;
957     //
958     if (fMultiRecKinks[label]>0){
959
960       //      for (Int_t j=0;j<TMath::Min(fMultiRecKinks[label],100);j++){
961       for (Int_t j=TMath::Min(fMultiRecKinks[label],Short_t(20))-1;j>=0;j--){
962         Int_t index = fIndexRecKinks[label*20+j];
963         //AliESDkink *kink2  = (AliESDkink*)fKinks->At(index);
964         AliESDkink *kink2  = (AliESDkink*)fEvent->GetKink(index);
965         if (TMath::Abs(kink2->GetLabel(0))==label &&TMath::Abs(kink2->GetLabel(1))==label2) {
966           fRecKinkInfo->fKinkMultiple++;
967           fSignedKinks[index]=1;
968           Int_t c0=0;
969           if (kink){
970             //      if (kink->fTRDOn) c0++;
971             //if (kink->fITSOn) c0++;
972             if (kink->GetStatus(2)>0) c0++;
973             if (kink->GetStatus(0)>0) c0++;
974           }
975           Int_t c2=0;
976           //      if (kink2->fTRDOn) c2++;
977           //if (kink2->fITSOn) c2++;
978           if (kink2->GetStatus(2)>0) c2++;
979           if (kink2->GetStatus(0)>0) c2++;
980
981           if (c2<c0) continue;
982           kink =kink2;
983         }
984         if (TMath::Abs(kink2->GetLabel(1))==label &&TMath::Abs(kink2->GetLabel(0))==label2) {
985           fRecKinkInfo->fKinkMultiple++;
986           fSignedKinks[index]=1;
987           Int_t c0=0;
988           if (kink){
989             //if (kink->fTRDOn) c0++;
990             //if (kink->fITSOn) c0++;
991             if (kink->GetStatus(2)>0) c0++;
992             if (kink->GetStatus(0)>0) c0++;
993
994           }
995           Int_t c2=0;
996           //      if (kink2->fTRDOn) c2++;
997           //if (kink2->fITSOn) c2++;
998           if (kink2->GetStatus(2)>0) c2++;
999           if (kink2->GetStatus(0)>0) c2++;
1000
1001           if (c2<c0) continue;
1002           kink =kink2;
1003         }
1004       }
1005     }
1006     if (kink){
1007       fRecKinkInfo->fKink = *kink;
1008       fRecKinkInfo->fRecStatus=1;
1009     }
1010     fTreeCmpKinks->Fill();
1011   }
1012   //  Int_t nkinks = fKinks->GetEntriesFast();
1013   Int_t nkinks = fEvent->GetNumberOfKinks();
1014   for (Int_t i=0;i<nkinks;i++){
1015     if (fSignedKinks[i]==0){
1016       //      AliESDkink *kink  = (AliESDkink*)fKinks->At(i);
1017       AliESDkink *kink  = (AliESDkink*)fEvent->GetKink(i);
1018       if (!kink) continue;
1019       //
1020       fRecKinkInfo->fKink = *kink;
1021       fRecKinkInfo->fRecStatus =-2;
1022       //
1023       AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(kink->GetLabel(0)));
1024       AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(kink->GetLabel(1)));
1025       if (fRecInfo1 && fRecInfo2){
1026         fRecKinkInfo->fT1 = (*fRecInfo1);
1027         fRecKinkInfo->fT2 = (*fRecInfo2);
1028         fRecKinkInfo->fRecStatus =-1;
1029       }
1030       fTreeCmpKinks->Fill();
1031     }
1032   }
1033
1034
1035   fTreeCmpKinks->AutoSave();
1036   printf("Time spended in BuilKinkInfo Loop\n");
1037   timer.Print();
1038   if (fDebug > 2) cerr<<"end of BuildKinkInfo Loop"<<endl;
1039   return 0;
1040 }
1041
1042
1043
1044 ////////////////////////////////////////////////////////////////////////
1045 ////////////////////////////////////////////////////////////////////////
1046 ////////////////////////////////////////////////////////////////////////
1047
1048
1049
1050 Int_t AliRecInfoMaker::BuildV0Info(Int_t eventNr)
1051 {
1052 //
1053 // loop over all entries for a given event, find corresponding 
1054 // rec. track and store in the fTreeCmp
1055 //
1056   static TDatabasePDG pdgtable;
1057
1058   TStopwatch timer;
1059   timer.Start();
1060   Int_t entry = fNextV0ToRead;
1061   Double_t nParticlesTR = fTreeGenV0->GetEntriesFast();
1062   cerr<<"fNParticles, nParticlesTR, fNextV0ToRead: "<<fNParticles<<" "
1063       <<nParticlesTR<<" "<<fNextV0ToRead<<endl;
1064   //
1065   TBranch * branch = fTreeCmpV0->GetBranch("RC.");
1066   branch->SetAddress(&fRecV0Info); // set all pointers
1067   const AliESDVertex * esdvertex = fEvent->GetVertex();
1068   Float_t vertex[3]= {esdvertex->GetXv(), esdvertex->GetYv(),esdvertex->GetZv()};
1069   
1070   //
1071   while (entry < nParticlesTR) {
1072     fTreeGenV0->GetEntry(entry);
1073     entry++;
1074     if (eventNr < fGenV0Info->GetMinus().fEventNr) continue;
1075     if (eventNr > fGenV0Info->GetMinus().fEventNr) continue;;
1076     //
1077     fNextV0ToRead = entry-1;
1078     //
1079     //
1080     AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(fGenV0Info->GetMinus().fLabel);
1081     AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(fGenV0Info->GetPlus().fLabel);
1082     if (fGenV0Info->GetMinus().fCharge*fGenV0Info->GetPlus().fCharge>0) continue;  // interactions
1083     if (!fRecInfo1 || !fRecInfo2) continue;
1084     fRecV0Info->fT1 = (*fRecInfo1);
1085     fRecV0Info->fT2 = (*fRecInfo2);
1086     fRecV0Info->fV0Status =0;
1087     if (fRecInfo1 && fRecInfo1->fStatus[1]>0) fRecV0Info->fV0Status+=1;
1088     if (fRecInfo2 && fRecInfo2->fStatus[1]>0) fRecV0Info->fV0Status+=2;
1089
1090     if (fRecV0Info->fV0Status==3&&fRecInfo1->fSign==fRecInfo2->fSign) fRecV0Info->fV0Status*=-1;
1091
1092
1093     if (abs(fRecV0Info->fV0Status)==3){
1094       fRecV0Info->Update(vertex);
1095       {
1096         //
1097         // TPC V0 Info
1098         Double_t x,alpha, param[5],cov[15];
1099         if ( fRecV0Info->fT1.GetESDtrack()->GetInnerParam() && fRecV0Info->fT2.GetESDtrack()->GetInnerParam()){
1100           fRecV0Info->fT1.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
1101           fRecV0Info->fT1.GetESDtrack()->GetInnerExternalCovariance(cov);
1102           AliExternalTrackParam paramP(x,alpha,param,cov);
1103           //
1104           fRecV0Info->fT2.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
1105           fRecV0Info->fT2.GetESDtrack()->GetInnerExternalCovariance(cov);
1106           AliExternalTrackParam paramM(x,alpha,param,cov);
1107           //
1108           fRecV0Info->fV0tpc->SetParamN(paramM);
1109           fRecV0Info->fV0tpc->SetParamP(paramP);
1110           Double_t pid1[5],pid2[5];
1111           fRecV0Info->fT1.GetESDtrack()->GetESDpid(pid1);
1112           fRecV0Info->fT1.GetESDtrack()->GetESDpid(pid2);
1113           //
1114           //fRecV0Info->fV0tpc.UpdatePID(pid1,pid2);
1115           fRecV0Info->fV0tpc->Update(vertex);
1116         
1117           //
1118           //
1119           fRecV0Info->fT1.GetESDtrack()->GetExternalParameters(x,param);
1120           fRecV0Info->fT1.GetESDtrack()->GetExternalCovariance(cov);
1121           alpha = fRecV0Info->fT1.GetESDtrack()->GetAlpha();
1122           new (&paramP) AliExternalTrackParam(x,alpha,param,cov);
1123           //
1124           fRecV0Info->fT2.GetESDtrack()->GetExternalParameters(x,param);
1125           fRecV0Info->fT2.GetESDtrack()->GetExternalCovariance(cov);
1126           alpha = fRecV0Info->fT2.GetESDtrack()->GetAlpha();
1127           new (&paramM) AliExternalTrackParam(x,alpha,param,cov);
1128           //
1129           fRecV0Info->fV0its->SetParamN(paramM);
1130           fRecV0Info->fV0its->SetParamP(paramP);
1131           //      fRecV0Info->fV0its.UpdatePID(pid1,pid2);
1132           fRecV0Info->fV0its->Update(vertex);
1133         }
1134       }
1135       //
1136       // ????
1137       // 
1138       if (TMath::Abs(fGenV0Info->GetMinus().fPdg)==11 &&TMath::Abs(fGenV0Info->GetPlus().fPdg)==11){
1139         if (fRecV0Info->fDist2>10){
1140           fRecV0Info->Update(vertex);
1141         }
1142         if (fRecV0Info->fDist2>10){
1143           fRecV0Info->Update(vertex);
1144         }
1145       }
1146     }   
1147     //
1148     // take the V0 from reconstruction
1149  
1150     Int_t label =  TMath::Min(fGenV0Info->GetMinus().fLabel,fGenV0Info->GetPlus().fLabel);
1151     Int_t label2 = TMath::Max(fGenV0Info->GetMinus().fLabel,fGenV0Info->GetPlus().fLabel);    
1152     AliV0 *v0MI=0;
1153     AliV0 *v0MIOff=0;
1154     fRecV0Info->fRecStatus   =0;
1155     fRecV0Info->fMultiple    = fMultiRecV0[label];
1156     fRecV0Info->fV0MultipleOn=0;
1157     fRecV0Info->fV0MultipleOff=0;
1158     //
1159     if (fMultiRecV0[label]>0 || fMultiRecV0[label2]>0){
1160
1161       //      for (Int_t j=0;j<TMath::Min(fMultiRecV0s[label],100);j++){
1162       for (Int_t j=TMath::Min(fMultiRecV0[label],Short_t(20))-1;j>=0;j--){
1163         Int_t index = fIndexRecV0[label*20+j];
1164         if (index<0) continue;
1165         AliV0 *v0MI2  = (AliV0*)fEvent->GetV0(index);
1166         // get track labels
1167         AliESDtrack * trackn = fEvent->GetTrack((v0MI2->GetNindex()));
1168         AliESDtrack * trackp = fEvent->GetTrack((v0MI2->GetPindex()));
1169         Int_t vlabeln = (trackn==0) ? -1 : trackn->GetLabel(); 
1170         Int_t vlabelp = (trackp==0) ? -1 : trackp->GetLabel(); 
1171         //
1172         if (TMath::Abs(vlabeln)==label &&TMath::Abs(vlabelp)==label2) {
1173           if (v0MI2->GetOnFlyStatus()) {
1174             v0MI =v0MI2;
1175             fRecV0Info->fV0MultipleOn++;
1176           }else  {
1177             v0MIOff = v0MI2;
1178             fRecV0Info->fV0MultipleOff++;
1179           }
1180           fSignedV0[index]=1;
1181         }
1182         if (TMath::Abs(vlabelp)==label &&TMath::Abs(vlabeln)==label2) {
1183           if (v0MI2->GetOnFlyStatus()){
1184             v0MI =v0MI2;
1185             fRecV0Info->fV0MultipleOn++;
1186           }else  {
1187             v0MIOff = v0MI2;
1188             fRecV0Info->fV0MultipleOff++;
1189           }
1190           fSignedV0[index]=1;
1191         }
1192       }
1193     }
1194     if (v0MI){
1195       new (fRecV0Info->fV0rec) AliV0(*v0MI);
1196       fRecV0Info->fRecStatus=1;
1197     }
1198     if (v0MIOff){
1199       new (fRecV0Info->fV0recOff) AliV0(*v0MIOff);
1200       fRecV0Info->fRecStatus=1;
1201     }
1202     Int_t mpdg = fGenV0Info->GetMother().GetPdgCode();
1203     Float_t mass = ( pdgtable.GetParticle(mpdg)==0) ? 0 :pdgtable.GetParticle(mpdg)->Mass();
1204     fRecV0Info->UpdateKF(*esdvertex,
1205                          fGenV0Info->GetPlus().GetPdg(),
1206                          fGenV0Info->GetMinus().GetPdg(),
1207                          mass);
1208     fTreeCmpV0->Fill();
1209   }
1210   //
1211   // write fake v0s
1212   //
1213   Int_t nV0MIs = fEvent->GetNumberOfV0s();
1214   for (Int_t i=0;i<nV0MIs;i++){
1215     if (fSignedV0[i]==0){
1216       AliV0 *v0MI  = (AliV0*)fEvent->GetV0(i);
1217       if (!v0MI) continue;
1218       //
1219       new (fRecV0Info->fV0rec) AliV0(*v0MI);
1220       fRecV0Info->fV0Status  =-10;
1221       fRecV0Info->fRecStatus =-2;
1222       //
1223       AliESDtrack * trackn = fEvent->GetTrack((v0MI->GetNindex()));
1224       AliESDtrack * trackp = fEvent->GetTrack((v0MI->GetPindex()));
1225       Int_t vlabeln = (trackn==0) ? -1 : trackn->GetLabel(); 
1226       Int_t vlabelp = (trackp==0) ? -1 : trackp->GetLabel(); 
1227
1228       AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(vlabeln));
1229       AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(vlabelp));
1230       if (fRecInfo1 && fRecInfo2){
1231         fRecV0Info->fT1 = (*fRecInfo1);
1232         fRecV0Info->fT2 = (*fRecInfo2);
1233         fRecV0Info->fRecStatus =-1;
1234       }
1235       fRecV0Info->Update(vertex);
1236       fRecV0Info->UpdateKF(*esdvertex,211,211,0.49767);
1237       fTreeCmpV0->Fill();
1238     }
1239   }
1240
1241
1242
1243   fTreeCmpV0->AutoSave();
1244   printf("Time spended in BuilV0Info Loop\n");
1245   timer.Print();
1246   if (fDebug > 2) cerr<<"end of BuildV0Info Loop"<<endl;
1247   return 0;
1248 }
1249 ////////////////////////////////////////////////////////////////////////
1250 ////////////////////////////////////////////////////////////////////////
1251
1252