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