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