]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGPP/TPC/AliRecInfoMaker.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fixed cmake files (Constantin)
[u/mrichter/AliRoot.git] / PWGPP / TPC / AliRecInfoMaker.cxx
CommitLineData
7cc34f08 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/*
27marian.ivanov@cern.ch
28Usage:
29
30
31
2bfe5463 32gSystem->Load("libPWGPP.so");
7cc34f08 33//
34AliRecInfoMaker *t2 = new AliRecInfoMaker("genTracks.root","cmpESDTracks.root","galice.root",0,0);
35t2->Exec();
36
37
38TFile f("cmpESDTracks.root");
39TTree * tree = (TTree*)f.Get("ESDcmpTracks");
40
41AliTreeDraw comp;
42comp.SetTree(tree)
43
44
45
46//
47//some cuts definition
48TCut 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");
51TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<5");
52TCut csec("csec","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)>0.5");
53
54
55TCut crec("crec","fReconstructed==1");
56TCut cteta1("cteta1","abs(MC.fParticle.Theta()/3.1415-0.5)<0.25");
57TCut cteta05("cteta05","abs(MC.fParticle.Theta()/3.1415-0.5)<0.1");
58
59TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
60TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50");
61TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)");
62TCut 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
67comp.T()->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)");
68comp.T()->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy");
69comp.T()->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy");
70comp.T()->SetAlias("theta","MC.fTrackRef.Theta()");
71comp.T()->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)");
72comp.T()->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
73comp.T()->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
74
75
76TH1F his("his","his",100,0,20);
77TH1F hpools("hpools","hpools",100,-7,7);
78TH1F hfake("hfake","hfake",1000,0,150);
79TProfile profp0("profp0","profp0",20,-0.4,0.9)
80
81comp.DrawXY("fTPCinP0[3]","fTPCDelta[4]/fTPCinP1[3]","fReconstructed==1"+cprim,"1",4,0.2,1.5,-0.06,0.06)
82comp.fRes->Draw();
83comp.fMean->Draw();
84
85comp.DrawXY("fITSinP0[3]","fITSDelta[4]/fITSinP1[3]","fReconstructed==1&&fITSOn"+cprim,"1",4,0.2,1.5,-0.06,0.06)
86comp.fRes->Draw();
87
88comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn",20,0.2,1.5)
89comp.fRes->Draw();
90
91comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn&&fITSOn&&fESDtrack.fITSFakeRatio<0.1",10,0.2,1.5)
92comp.fRes->Draw();
93comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn&&fITSOn&&fESDtrack.fITSFakeRatio>0.1",10,0.2,1.5)
94comp.fRes->Draw();
95
96comp.T()->Draw("fESDtrack.fITSsignal/fESDtrack.fTPCsignal","fITSOn&&fTPCOn&&fESDtrack.fITSFakeRatio==0")
97
98TH1F his("his","his",100,0,20);
99TH1F hpools("hpools","hpools",100,-7,7);
100
101TH2F * hdedx0 = new TH2F("dEdx0","dEdx0",100, 0,2,200,0,550); hdedx0->SetMarkerColor(1);
102TH2F * hdedx1 = new TH2F("dEdx1","dEdx1",100, 0,2,200,0,550); hdedx1->SetMarkerColor(4);
103TH2F * hdedx2 = new TH2F("dEdx2","dEdx2",100, 0,2,200,0,550); hdedx2->SetMarkerColor(3);
104TH2F * hdedx3 = new TH2F("dEdx3","dEdx3",100, 0,2,200,0,550); hdedx3->SetMarkerColor(2);
105
106comp.T()->Draw("fESDtrack.fITSsignal:MC.fParticle.P()>>dEdx0","fITSOn&&abs(fPdg)==211&&fITStrack.fN==6"+cprim)
107comp.T()->Draw("fESDtrack.fITSsignal:MC.fParticle.P()>>dEdx1","fITSOn&&abs(fPdg)==2212&&fITStrack.fN==6"+cprim)
108comp.T()->Draw("fESDtrack.fITSsignal:MC.fParticle.P()>>dEdx2","fITSOn&&abs(fPdg)==321&&fITStrack.fN==6"+cprim)
109comp.T()->Draw("fESDtrack.fITSsignal:MC.fParticle.P()>>dEdx3","fITSOn&&abs(fPdg)==11&&fITStrack.fN==6"+cprim)
110
111
112comp.T()->Draw("fESDtrack.fTRDsignal:MC.fParticle.P()>>dEdx0","fTRDOn&&abs(fPdg)==211&&fTRDtrack.fN>40&&fStatus[2]>1")
113comp.T()->Draw("fESDtrack.fTRDsignal:MC.fParticle.P()>>dEdx1","fTRDOn&&abs(fPdg)==2212&&fTRDtrack.fN>40&&fStatus[2]>1")
114comp.T()->Draw("fESDtrack.fTRDsignal:MC.fParticle.P()>>dEdx2","fTRDOn&&abs(fPdg)==321&&fTRDtrack.fN>40&&fStatus[2]>1")
115comp.T()->Draw("fESDtrack.fTRDsignal:MC.fParticle.P()>>dEdx3","fTRDOn&&abs(fPdg)==11&&fTRDtrack.fN>40&&fStatus[2]>1")
116
117comp.T()->Draw("fESDtrack.fTPCsignal:fTPCinP0[4]>>dEdx0","fTPCOn&&abs(fPdg)==211&&fESDtrack.fTPCncls>180&&fESDtrack.fTPCsignal>10"+cteta1);
118comp.T()->Draw("fESDtrack.fTPCsignal:fTPCinP0[4]>>dEdx1","fTPCOn&&abs(fPdg)==2212&&fESDtrack.fTPCncls>180&&fESDtrack.fTPCsignal>10"+cteta1);
119comp.T()->Draw("fESDtrack.fTPCsignal:fTPCinP0[4]>>dEdx2","fTPCOn&&abs(fPdg)==321&&fESDtrack.fTPCncls>180&&fESDtrack.fTPCsignal>10"+cteta1);
120comp.T()->Draw("fESDtrack.fTPCsignal:fTPCinP0[4]>>dEdx3","fTPCOn&&abs(fPdg)==11&&fESDtrack.fTPCncls>180&&fESDtrack.fTPCsignal>10"+cteta1);
121
122hdedx3->SetXTitle("P(GeV/c)");
123hdedx3->SetYTitle("dEdx(unit)");
124hdedx3->Draw(); hdedx1->Draw("same"); hdedx2->Draw("same"); hdedx0->Draw("same");
125
126comp.DrawXY("fITSinP0[3]","fITSPools[4]","fReconstructed==1&&fPdg==-211&&fITSOn"+cprim,"1",4,0.2,1.0,-8,8)
127
128TProfile 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
c64cb1f6 181using std::cout;
182using std::cerr;
183using std::endl;
184
7cc34f08 185ClassImp(AliRecInfoMaker)
186
187
188
189
190AliTPCParam * AliRecInfoMaker::GetTPCParam(){
191 //
192 // create TPC param
193 //
194 AliTPCParamSR * par = new AliTPCParamSR;
195 par->Update();
196 return par;
197}
198
199
200
201void 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////////////////////////////////////////////////////////////////////////
230AliRecInfoMaker::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;
18e588e9 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);
7cc34f08 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////////////////////////////////////////////////////////////////////////
335AliRecInfoMaker::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 //
18e588e9 386 memset(fFnGenTracks,0,sizeof(fFnGenTracks));
387 memset(fFnCmp,0,sizeof(fFnCmp));
7cc34f08 388}
389
390
391
392
393////////////////////////////////////////////////////////////////////////
394AliRecInfoMaker::~AliRecInfoMaker()
395{
396 //
397 // Destructor
398 //
399 if (fLoader) {
400 delete fLoader;
401 }
402}
403
404//////////////////////////////////////////////////////////////
405Int_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
423Int_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////////////////////////////////////////////////////////////////////////
458void 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////////////////////////////////////////////////////////////////////////
477Int_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////////////////////////////////////////////////////////////////////////
488Int_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////////////////////////////////////////////////////////////////////////
561Bool_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////////////////////////////////////////////////////////////////////////
613void 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////////////////////////////////////////////////////////////////////////
655void 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
675TVector3 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
690Int_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////////////////////////////////////////////////////////////////////////
798Int_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]);
18e588e9 839 if(!track) continue;
7cc34f08 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 //
820afb27 888 fRecInfo->SetESDtrack(track);
889
890 /*
7cc34f08 891 if (track) {
892 fRecInfo->SetESDtrack(track);
893 }else{
894 fRecInfo->SetESDtrack(&dummytrack);
895 }
820afb27 896 */
7cc34f08 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////////////////////////////////////////////////////////////////////////
925Int_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);
18e588e9 952 if(!fRecInfo1) continue;
7cc34f08 953 AliESDRecInfo* fRecInfo2 = (AliESDRecInfo*)fRecArray->At(fGenKinkInfo->GetPlus().fLabel);
18e588e9 954 if(!fRecInfo2) continue;
7cc34f08 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
1065Int_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]= {esdvertex->GetXv(), esdvertex->GetYv(),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