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