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