]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliESDComparisonMI.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Removing unnecessary warning lines
[u/mrichter/AliRoot.git] / STEER / AliESDComparisonMI.C
CommitLineData
ae7d73d2 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
e1eaa672 26/*
ae7d73d2 27marian.ivanov@cern.ch
28Usage:
e1eaa672 29
ae7d73d2 30
b3bf848f 31.L $ALICE_ROOT/STEER/AliGenInfo.C+
ae7d73d2 32//be sure you created genTracks file before
b3bf848f 33.L $ALICE_ROOT/STEER/AliESDComparisonMI.C+
ae7d73d2 34//
e1eaa672 35ESDCmpTr *t2 = new ESDCmpTr("genTracks.root","cmpESDTracks.root","galice.root",-1,0,0);
ae7d73d2 36t2->Exec();
37
38//
39//some cuts definition
51ad6848 40TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.01&&abs(MC.fVDist[2])<0.01")
ae7d73d2 41//TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.5&&abs(MC.fVDist[2])<0.5")
51ad6848 42//TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<3.9");
f007cb5f 43TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<5");
ae7d73d2 44TCut csec("csec","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)>0.5");
45
46
47TCut crec("crec","fReconstructed==1");
51ad6848 48TCut cteta1("cteta1","abs(MC.fParticle.Theta()/3.1415-0.5)<0.25");
f007cb5f 49TCut cteta05("cteta05","abs(MC.fParticle.Theta()/3.1415-0.5)<0.1");
50
ae7d73d2 51TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
52TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50");
53TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)");
54TCut cchi2("cchi2","fESDTrack.fITSchi2MIP[0]<7.&&fESDTrack.fITSchi2MIP[1]<5.&&fESDTrack.fITSchi2MIP[2]<7.&&fESDTrack.fITSchi2MIP[3]<7.5&&fESDTrack.fITSchi2MIP[4]<6.")
55
56AliESDComparisonDraw comp;
57comp.SetIO();
51ad6848 58TFile f("genHits.root");
59TTree * treel = (TTree*)f.Get("HitLines");
60if (treel) comp->fTree->AddFriend(treel,"L");
ae7d73d2 61
62//
63//example
64comp.fTree->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)");
51ad6848 65comp.fTree->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy");
66comp.fTree->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy");
67comp.fTree->SetAlias("theta","MC.fTrackRef.Theta()");
68comp.fTree->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)");
51ad6848 69comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
f007cb5f 70comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
71
ae7d73d2 72
73TH1F his("his","his",100,0,20);
74TH1F hpools("hpools","hpools",100,-7,7);
75TH1F hfake("hfake","hfake",1000,0,150);
76TProfile profp0("profp0","profp0",20,-0.4,0.9)
77
f007cb5f 78comp.DrawXY("fTPCinP0[3]","fTPCDelta[4]/fTPCinP1[3]","fReconstructed==1"+cprim,"1",4,0.2,1.5,-0.06,0.06)
ae7d73d2 79comp.fRes->Draw();
80comp.fMean->Draw();
51ad6848 81
f007cb5f 82comp.DrawXY("fITSinP0[3]","fITSDelta[4]/fITSinP1[3]","fReconstructed==1&&fITSOn"+cprim,"1",4,0.2,1.5,-0.06,0.06)
51ad6848 83comp.fRes->Draw();
84
85comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn",20,0.2,1.5)
86comp.fRes->Draw();
87
88comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn&&fITSOn&&fESDTrack.fITSFakeRatio<0.1",10,0.2,1.5)
89comp.fRes->Draw();
90comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn&&fITSOn&&fESDTrack.fITSFakeRatio>0.1",10,0.2,1.5)
ae7d73d2 91comp.fRes->Draw();
92
93comp.fTree->Draw("fESDTrack.fITSsignal/fESDTrack.fTPCsignal","fITSOn&&fTPCOn&&fESDTrack.fITSFakeRatio==0")
94
95TH1F his("his","his",100,0,20);
96TH1F hpools("hpools","hpools",100,-7,7);
97
f007cb5f 98TH2F * hdedx0 = new TH2F("dEdx0","dEdx0",100, 0,2,200,0,550); hdedx0->SetMarkerColor(1);
99TH2F * hdedx1 = new TH2F("dEdx1","dEdx1",100, 0,2,200,0,550); hdedx1->SetMarkerColor(4);
100TH2F * hdedx2 = new TH2F("dEdx2","dEdx2",100, 0,2,200,0,550); hdedx2->SetMarkerColor(3);
101TH2F * hdedx3 = new TH2F("dEdx3","dEdx3",100, 0,2,200,0,550); hdedx3->SetMarkerColor(2);
102
103comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx0","fITSOn&&abs(fPdg)==211&&fITStrack.fN==6"+cprim)
104comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx1","fITSOn&&abs(fPdg)==2212&&fITStrack.fN==6"+cprim)
105comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx2","fITSOn&&abs(fPdg)==321&&fITStrack.fN==6"+cprim)
106comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx3","fITSOn&&abs(fPdg)==11&&fITStrack.fN==6"+cprim)
51ad6848 107
51ad6848 108
f007cb5f 109comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx0","fTRDOn&&abs(fPdg)==211&&fTRDtrack.fN>40&&fStatus[2]>1")
110comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx1","fTRDOn&&abs(fPdg)==2212&&fTRDtrack.fN>40&&fStatus[2]>1")
111comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx2","fTRDOn&&abs(fPdg)==321&&fTRDtrack.fN>40&&fStatus[2]>1")
112comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx3","fTRDOn&&abs(fPdg)==11&&fTRDtrack.fN>40&&fStatus[2]>1")
51ad6848 113
e1eaa672 114comp.fTree->Draw("fESDTrack.fTPCsignal:fTPCinP0[4]>>dEdx0","fTPCOn&&abs(fPdg)==211&&fESDTrack.fTPCncls>180&&fESDTrack.fTPCsignal>10"+cteta1);
115comp.fTree->Draw("fESDTrack.fTPCsignal:fTPCinP0[4]>>dEdx1","fTPCOn&&abs(fPdg)==2212&&fESDTrack.fTPCncls>180&&fESDTrack.fTPCsignal>10"+cteta1);
116comp.fTree->Draw("fESDTrack.fTPCsignal:fTPCinP0[4]>>dEdx2","fTPCOn&&abs(fPdg)==321&&fESDTrack.fTPCncls>180&&fESDTrack.fTPCsignal>10"+cteta1);
117comp.fTree->Draw("fESDTrack.fTPCsignal:fTPCinP0[4]>>dEdx3","fTPCOn&&abs(fPdg)==11&&fESDTrack.fTPCncls>180&&fESDTrack.fTPCsignal>10"+cteta1);
51ad6848 118
e1eaa672 119hdedx3->SetXTitle("P(GeV/c)");
120hdedx3->SetYTitle("dEdx(unit)");
f007cb5f 121hdedx3->Draw(); hdedx1->Draw("same"); hdedx2->Draw("same"); hdedx0->Draw("same");
ae7d73d2 122
123comp.DrawXY("fITSinP0[3]","fITSPools[4]","fReconstructed==1&&fPdg==-211&&fITSOn"+cprim,"1",4,0.2,1.0,-8,8)
124
f007cb5f 125TProfile prof("prof","prof",10,0.5,5);
126
127
128
129
ae7d73d2 130*/
131
132
133#if !defined(__CINT__) || defined(__MAKECINT__)
134#include <stdio.h>
135#include <string.h>
136//ROOT includes
137#include "Rtypes.h"
138#include "TFile.h"
139#include "TTree.h"
140#include "TChain.h"
141#include "TCut.h"
142#include "TString.h"
143#include "TBenchmark.h"
144#include "TStopwatch.h"
145#include "TParticle.h"
146#include "TSystem.h"
147#include "TTimer.h"
148#include "TVector3.h"
149#include "TPad.h"
150#include "TCanvas.h"
151#include "TH1F.h"
152#include "TH2F.h"
153#include "TF1.h"
154#include "TText.h"
155#include "Getline.h"
156#include "TStyle.h"
157
158//ALIROOT includes
159#include "AliRun.h"
160#include "AliStack.h"
161#include "AliESDtrack.h"
162#include "AliSimDigits.h"
163#include "AliTPCParam.h"
164#include "AliTPC.h"
165#include "AliTPCLoader.h"
166#include "AliDetector.h"
167#include "AliTrackReference.h"
168#include "AliRun.h"
169#include "AliTPCParamSR.h"
170#include "AliTracker.h"
171#include "AliComplexCluster.h"
86ad5fcb 172#include "AliESDEvent.h"
ae7d73d2 173#include "AliESDtrack.h"
51ad6848 174#include "AliITStrackMI.h"
ae7d73d2 175#include "AliTRDtrack.h"
51ad6848 176#include "AliHelix.h"
177#include "AliESDVertex.h"
178#include "AliExternalTrackParam.h"
179#include "AliESDkink.h"
d6a49f20 180#include "AliESDv0.h"
51ad6848 181
ae7d73d2 182#endif
183#include "AliGenInfo.h"
184#include "AliESDComparisonMI.h"
185
186
187
81e97e0d 188
189
190void MakeAliases(AliESDComparisonDraw&comp)
191{
192 //
193 // aliases definition
194 //
195 comp.fTree->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)");
196 comp.fTree->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy");
197 comp.fTree->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy");
198 comp.fTree->SetAlias("theta","MC.fTrackRef.Theta()");
199 comp.fTree->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)");
200 comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
201 comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
202
203 comp.fTree->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.");
204
205 comp.fTree->SetAlias("dtofmc2","fESDTrack.fTrackTime[2]-(10^12*MC.fTOFReferences[0].fTime)");
206 comp.fTree->SetAlias("dtofrc2","(fESDTrack.fTrackTime[2]-fESDTrack.fTOFsignal)");
207
208 comp.fTree->SetAlias("psum","fESDTrack.fTOFr[4]+fESDTrack.fTOFr[3]+fESDTrack.fTOFr[2]+fESDTrack.fTOFr[1]+fESDTrack.fTOFr[0]");
209 comp.fTree->SetAlias("P0","fESDTrack.fTOFr[0]/psum");
210 comp.fTree->SetAlias("P1","fESDTrack.fTOFr[1]/psum");
211 comp.fTree->SetAlias("P2","fESDTrack.fTOFr[2]/psum");
212 comp.fTree->SetAlias("P3","fESDTrack.fTOFr[3]/psum");
213 comp.fTree->SetAlias("P4","fESDTrack.fTOFr[4]/psum");
214 comp.fTree->SetAlias("MaxP","max(max(max(P0,P1),max(P2,P3)),P4)");
215}
216
217
218
51ad6848 219
220void AliESDRecInfo::UpdatePoints(AliESDtrack*track)
221{
222 //
223 //
224 Int_t iclusters[200];
225 Float_t density[160];
226 for (Int_t i=0;i<160;i++) density[i]=-1.;
227 fTPCPoints[0]= 160;
228 fTPCPoints[1] = -1;
229 //
230 if (fTPCPoints[0]<fTPCPoints[1]) return;
231 // Int_t nclusters=track->GetTPCclusters(iclusters);
232
233 Int_t ngood=0;
234 Int_t undeff=0;
235 Int_t nall =0;
236 Int_t range=20;
237 for (Int_t i=0;i<160;i++){
238 Int_t last = i-range;
239 if (nall<range) nall++;
240 if (last>=0){
241 if (iclusters[last]>0&& (iclusters[last]&0x8000)==0) ngood--;
242 if (iclusters[last]==-1) undeff--;
243 }
244 if (iclusters[i]>0&& (iclusters[i]&0x8000)==0) ngood++;
245 if (iclusters[i]==-1) undeff++;
246 if (nall==range &&undeff<range/2) density[i-range/2] = Float_t(ngood)/Float_t(nall-undeff);
247 }
248 Float_t maxdens=0;
249 Int_t indexmax =0;
250 for (Int_t i=0;i<160;i++){
251 if (density[i]<0) continue;
252 if (density[i]>maxdens){
253 maxdens=density[i];
254 indexmax=i;
255 }
256 }
257 //
258 //max dens point
259 fTPCPoints[3] = maxdens;
260 fTPCPoints[1] = indexmax;
261 //
262 // last point
263 for (Int_t i=indexmax;i<160;i++){
264 if (density[i]<0) continue;
265 if (density[i]<maxdens/2.) {
266 break;
267 }
268 fTPCPoints[2]=i;
269 }
270 //
271 // first point
272 for (Int_t i=indexmax;i>0;i--){
273 if (density[i]<0) continue;
274 if (density[i]<maxdens/2.) {
275 break;
276 }
277 fTPCPoints[0]=i;
278 }
279 //
f007cb5f 280 // Density at the last 30 padrows
281 //
282 //
283 nall = 0;
284 ngood = 0;
285 for (Int_t i=159;i>0;i--){
286 if (iclusters[i]==-1) continue; //dead zone
287 nall++;
288 if (iclusters[i]>0) ngood++;
289 if (nall>20) break;
290 }
291 fTPCPoints[4] = Float_t(ngood)/Float_t(nall);
292 //
51ad6848 293 if ((track->GetStatus()&AliESDtrack::kITSrefit)>0) fTPCPoints[0]=-1;
f007cb5f 294 //
295 //
296 // check TRDPoints
297 /*
298 nclusters=track->GetTRDclusters(iclusters);
299 for (Int_t i=nclusters;i>0;i--){
300
301 }
302 */
303
304
51ad6848 305}
306
ae7d73d2 307//
308//
86ad5fcb 309void AliESDRecInfo::Update(AliMCInfo* info,AliTPCParam * /*par*/, Bool_t reconstructed, AliESDEvent *event)
ae7d73d2 310{
311 //
312 //
313 //calculates derived variables
314 //
315 //
51ad6848 316 UpdatePoints(&fESDTrack);
317 fBestTOFmatch=1000;
ae7d73d2 318 AliTrackReference * ref = &(info->fTrackRef);
319 fTPCinR0[0] = info->fTrackRef.X();
320 fTPCinR0[1] = info->fTrackRef.Y();
321 fTPCinR0[2] = info->fTrackRef.Z();
322 fTPCinR0[3] = TMath::Sqrt(fTPCinR0[0]*fTPCinR0[0]+fTPCinR0[1]*fTPCinR0[1]);
323 fTPCinR0[4] = TMath::ATan2(fTPCinR0[1],fTPCinR0[0]);
324 //
325 fTPCinP0[0] = ref->Px();
326 fTPCinP0[1] = ref->Py();
327 fTPCinP0[2] = ref->Pz();
328 fTPCinP0[3] = ref->Pt();
329 fTPCinP0[4] = ref->P();
330 fDeltaP = (ref->P()-info->fParticle.P())/info->fParticle.P();
331 //
332 //
333 if (fTPCinP0[3]>0.0000001){
334 //
335 fTPCAngle0[0] = TMath::ATan2(fTPCinP0[1],fTPCinP0[0]);
336 fTPCAngle0[1] = TMath::ATan(fTPCinP0[2]/fTPCinP0[3]);
337 }
338 //
339 //
340 fITSinP0[0]=info->fParticle.Px();
341 fITSinP0[1]=info->fParticle.Py();
342 fITSinP0[2]=info->fParticle.Pz();
343 fITSinP0[3]=info->fParticle.Pt();
344 //
345 fITSinR0[0]=info->fParticle.Vx();
346 fITSinR0[1]=info->fParticle.Vy();
347 fITSinR0[2]=info->fParticle.Vz();
348 fITSinR0[3] = TMath::Sqrt(fITSinR0[0]*fITSinR0[0]+fITSinR0[1]*fITSinR0[1]);
349 fITSinR0[4] = TMath::ATan2(fITSinR0[1],fITSinR0[0]);
350 //
351 //
352 if (fITSinP0[3]>0.0000001){
353 fITSAngle0[0] = TMath::ATan2(fITSinP0[1],fITSinP0[0]);
354 fITSAngle0[1] = TMath::ATan(fITSinP0[2]/fITSinP0[3]);
355 }
356 //
51ad6848 357 for (Int_t i=0;i<4;i++) fStatus[i] =0;
358 fReconstructed = kFALSE;
359 fTPCOn = kFALSE;
360 fITSOn = kFALSE;
361 fTRDOn = kFALSE;
362 if (reconstructed==kFALSE) return;
363
364 fLabels[0] = info->fLabel;
365 fLabels[1] = info->fPrimPart;
ae7d73d2 366 fReconstructed = kTRUE;
367 fTPCOn = ((fESDTrack.GetStatus()&AliESDtrack::kTPCrefit)>0) ? kTRUE : kFALSE;
368 fITSOn = ((fESDTrack.GetStatus()&AliESDtrack::kITSrefit)>0) ? kTRUE : kFALSE;
369 fTRDOn = ((fESDTrack.GetStatus()&AliESDtrack::kTRDrefit)>0) ? kTRUE : kFALSE;
51ad6848 370 //
371 //
372 if ((fESDTrack.GetStatus()&AliESDtrack::kTPCrefit)>0){
373 fStatus[1] =3;
374 }
375 else{
376 if ((fESDTrack.GetStatus()&AliESDtrack::kTPCout)>0){
377 fStatus[1] =2;
378 }
379 else{
98808dce 380 if ((fESDTrack.GetStatus()&AliESDtrack::kTPCin)>0)
381 fStatus[1]=1;
51ad6848 382 }
383 }
f007cb5f 384 //
385 if ((fESDTrack.GetStatus()&AliESDtrack::kITSout)>0){
386 fStatus[0] =2;
387 }
388 else{
389 if ((fESDTrack.GetStatus()&AliESDtrack::kITSrefit)>0){
390 fStatus[0] =1;
391 }
392 else{
393 fStatus[0]=0;
394 }
395 }
396
51ad6848 397 //
398 //
399 if ((fESDTrack.GetStatus()&AliESDtrack::kTRDrefit)>0){
400 fStatus[2] =2;
401 }
402 else{
403 if ((fESDTrack.GetStatus()&AliESDtrack::kTRDout)>0){
404 fStatus[2] =1;
405 }
406 }
407 if ((fESDTrack.GetStatus()&AliESDtrack::kTRDStop)>0){
408 fStatus[2] =10;
409 }
410
411 //
412 //TOF
413 //
414 if (((fESDTrack.GetStatus()&AliESDtrack::kTOFout)>0)){
415 //
416 // best tof match
417 Double_t times[5];
418 fESDTrack.GetIntegratedTimes(times);
419 for (Int_t i=0;i<5;i++){
420 if ( TMath::Abs(fESDTrack.GetTOFsignal()-times[i]) <TMath::Abs(fBestTOFmatch) ){
421 fBestTOFmatch = fESDTrack.GetTOFsignal()-times[i];
422 }
423 }
424 Int_t toflabel[3];
425 fESDTrack.GetTOFLabel(toflabel);
426 Bool_t toffake=kTRUE;
f007cb5f 427 Bool_t tofdaughter=kFALSE;
51ad6848 428 for (Int_t i=0;i<3;i++){
f007cb5f 429 if (toflabel[i]<0) continue;
51ad6848 430 if (toflabel[i]== TMath::Abs(fESDTrack.GetLabel())) toffake=kFALSE;
f007cb5f 431 if (toflabel[i]==info->fParticle.GetDaughter(0) || (toflabel[i]==info->fParticle.GetDaughter(1))) tofdaughter=kTRUE; // decay product of original particle
51ad6848 432 fStatus[3]=1;
433 }
f007cb5f 434 if (toffake) fStatus[3] =3; //total fake
435 if (tofdaughter) fStatus[3]=2; //fake because of decay
51ad6848 436 }else{
437 fStatus[3]=0;
438 }
439
440
81e97e0d 441 if (fStatus[1]>0 &&info->fNTPCRef>0&&TMath::Abs(fTPCinP0[3])>0.0001){
ae7d73d2 442 //TPC
443 fESDTrack.GetInnerXYZ(fTPCinR1);
444 fTPCinR1[3] = TMath::Sqrt(fTPCinR1[0]*fTPCinR1[0]+fTPCinR1[1]*fTPCinR1[1]);
445 fTPCinR1[4] = TMath::ATan2(fTPCinR1[1],fTPCinR1[0]);
446 fESDTrack.GetInnerPxPyPz(fTPCinP1);
447 fTPCinP1[3] = TMath::Sqrt(fTPCinP1[0]*fTPCinP1[0]+fTPCinP1[1]*fTPCinP1[1]);
448 fTPCinP1[4] = TMath::Sqrt(fTPCinP1[3]*fTPCinP1[3]+fTPCinP1[2]*fTPCinP1[2]);
449 //
450 //
0a1e73ec 451 if (fTPCinP1[3]>0.000000000000001){
ae7d73d2 452 fTPCAngle1[0] = TMath::ATan2(fTPCinP1[1],fTPCinP1[0]);
453 fTPCAngle1[1] = TMath::ATan(fTPCinP1[2]/fTPCinP1[3]);
454 }
e1d4c1b5 455 Double_t cov[15], param[5],x, alpha;
ae7d73d2 456 fESDTrack.GetInnerExternalCovariance(cov);
e1d4c1b5 457 fESDTrack.GetInnerExternalParameters(alpha, x,param);
f007cb5f 458 if (x<50) return ;
ae7d73d2 459 //
460 fTPCDelta[0] = (fTPCinR0[4]-fTPCinR1[4])*fTPCinR1[3]; //delta rfi
461 fTPCPools[0] = fTPCDelta[0]/TMath::Sqrt(cov[0]);
462 fTPCDelta[1] = (fTPCinR0[2]-fTPCinR1[2]); //delta z
463 fTPCPools[1] = fTPCDelta[1]/TMath::Sqrt(cov[2]);
464 fTPCDelta[2] = (fTPCAngle0[0]-fTPCAngle1[0]);
465 fTPCPools[2] = fTPCDelta[2]/TMath::Sqrt(cov[5]);
466 fTPCDelta[3] = (TMath::Tan(fTPCAngle0[1])-TMath::Tan(fTPCAngle1[1]));
467 fTPCPools[3] = fTPCDelta[3]/TMath::Sqrt(cov[9]);
468 fTPCDelta[4] = (fTPCinP0[3]-fTPCinP1[3]);
469 Double_t sign = (param[4]>0)? 1.:-1;
51ad6848 470 fSign =sign;
98808dce 471 fTPCPools[4] = sign*(1./fTPCinP0[3]-1./fTPCinP1[3])/TMath::Sqrt(TMath::Abs(cov[14]));
ae7d73d2 472 }
473 if (fITSOn){
474 // ITS
475 Double_t param[5],x;
f007cb5f 476 fESDTrack.GetExternalParameters(x,param);
477 // fESDTrack.GetConstrainedExternalParameters(x,param);
ae7d73d2 478 Double_t cov[15];
479 fESDTrack.GetExternalCovariance(cov);
f007cb5f 480 //fESDTrack.GetConstrainedExternalCovariance(cov);
ae7d73d2 481 if (TMath::Abs(param[4])<0.0000000001) return;
482
483 fESDTrack.GetXYZ(fITSinR1);
484 fESDTrack.GetPxPyPz(fITSinP1);
485 fITSinP1[3] = TMath::Sqrt(fITSinP1[0]*fITSinP1[0]+fITSinP1[1]*fITSinP1[1]);
486 //
487 fITSinR1[3] = TMath::Sqrt(fITSinR1[0]*fITSinR1[0]+fITSinR1[1]*fITSinR1[1]);
488 fITSinR1[4] = TMath::ATan2(fITSinR1[1],fITSinR1[0]);
489 //
490 //
491 if (fITSinP1[3]>0.0000001){
492 fITSAngle1[0] = TMath::ATan2(fITSinP1[1],fITSinP1[0]);
493 fITSAngle1[1] = TMath::ATan(fITSinP1[2]/fITSinP1[3]);
494 }
495 //
496 //
497 fITSDelta[0] = (fITSinR0[4]-fITSinR1[4])*fITSinR1[3]; //delta rfi
498 fITSPools[0] = fITSDelta[0]/TMath::Sqrt(cov[0]);
499 fITSDelta[1] = (fITSinR0[2]-fITSinR1[2]); //delta z
500 fITSPools[1] = fITSDelta[1]/TMath::Sqrt(cov[2]);
501 fITSDelta[2] = (fITSAngle0[0]-fITSAngle1[0]);
502 fITSPools[2] = fITSDelta[2]/TMath::Sqrt(cov[5]);
503 fITSDelta[3] = (TMath::Tan(fITSAngle0[1])-TMath::Tan(fITSAngle1[1]));
504 fITSPools[3] = fITSDelta[3]/TMath::Sqrt(cov[9]);
505 fITSDelta[4] = (fITSinP0[3]-fITSinP1[3]);
506 Double_t sign = (param[4]>0) ? 1:-1;
51ad6848 507 fSign = sign;
ae7d73d2 508 fITSPools[4] = sign*(1./fITSinP0[3]-1./fITSinP1[3])/TMath::Sqrt(cov[14]);
509 }
510
511}
512
513
514void AliESDRecV0Info::Update(Float_t vertex[3])
51ad6848 515{
516
517 if ( (fT1.fStatus[1]>0)&& (fT2.fStatus[1]>0)){
518 Float_t distance1,distance2;
519 Double_t xx[3],pp[3];
520 //
521 Double_t xd[3],pd[3],signd;
522 Double_t xm[3],pm[3],signm;
f007cb5f 523 //
524 //
525 if (fT1.fITSOn&&fT2.fITSOn){
526 for (Int_t i=0;i<3;i++){
527 xd[i] = fT2.fITSinR1[i];
528 pd[i] = fT2.fITSinP1[i];
529 xm[i] = fT1.fITSinR1[i];
530 pm[i] = fT1.fITSinP1[i];
531 }
532 }
533 else{
534
535 for (Int_t i=0;i<3;i++){
536 xd[i] = fT2.fTPCinR1[i];
537 pd[i] = fT2.fTPCinP1[i];
538 xm[i] = fT1.fTPCinR1[i];
539 pm[i] = fT1.fTPCinP1[i];
540 }
51ad6848 541 }
f007cb5f 542 //
543 //
51ad6848 544 signd = fT2.fSign<0 ? -1:1;
545 signm = fT1.fSign<0 ? -1:1;
546
547 AliHelix dhelix1(xd,pd,signd);
548 dhelix1.GetMomentum(0,pp,0);
549 dhelix1.Evaluate(0,xx);
550 //
551 // Double_t x2[3],p2[3];
552 //
553 AliHelix mhelix(xm,pm,signm);
554 //
555 //find intersection linear
556 //
557 Double_t phase[2][2],radius[2];
558 Int_t points = dhelix1.GetRPHIintersections(mhelix, phase, radius,200);
559 Double_t delta1=10000,delta2=10000;
560
561 if (points==1){
f007cb5f 562 fRs[0] = TMath::Sqrt(radius[0]);
563 fRs[1] = TMath::Sqrt(radius[0]);
51ad6848 564 }
565 if (points==2){
f007cb5f 566 fRs[0] =TMath::Min(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1]));
567 fRs[1] =TMath::Max(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1]));
51ad6848 568 }
569
570 if (points>0){
571 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
572 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
573 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
574 }
575 if (points==2){
576 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
577 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
578 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
579 }
580 if (points==1){
f007cb5f 581 fRs[0] = TMath::Sqrt(radius[0]);
582 fRs[1] = TMath::Sqrt(radius[0]);
51ad6848 583 fDistMinR = delta1;
584 }
585 if (points==2){
586 if (radius[0]<radius[1]){
f007cb5f 587 fRs[0] = TMath::Sqrt(radius[0]);
588 fRs[1] = TMath::Sqrt(radius[1]);
51ad6848 589 fDistMinR = delta1;
590 }
591 else{
f007cb5f 592 fRs[0] = TMath::Sqrt(radius[1]);
593 fRs[1] = TMath::Sqrt(radius[0]);
51ad6848 594 fDistMinR = delta2;
595 }
596 }
597 //
598 //
599 distance1 = TMath::Min(delta1,delta2);
600 //
601 //find intersection parabolic
602 //
603 points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
604 delta1=10000,delta2=10000;
605
606 if (points>0){
607 dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
608 }
609 if (points==2){
610 dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
611 }
612
613 distance2 = TMath::Min(delta1,delta2);
f007cb5f 614 if (distance2>100) fDist2 =100;
615 return;
51ad6848 616 if (delta1<delta2){
617 //get V0 info
618 dhelix1.Evaluate(phase[0][0],fXr);
619 dhelix1.GetMomentum(phase[0][0],fPdr);
620 mhelix.GetMomentum(phase[0][1],fPm);
621 dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle);
622 fRr = TMath::Sqrt(radius[0]);
623 }
624 else{
625 dhelix1.Evaluate(phase[1][0],fXr);
626 dhelix1.GetMomentum(phase[1][0], fPdr);
627 mhelix.GetMomentum(phase[1][1], fPm);
628 dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle);
629 fRr = TMath::Sqrt(radius[1]);
630 }
631 fDist1 = TMath::Sqrt(distance1);
632 fDist2 = TMath::Sqrt(distance2);
633
634 if (fDist2<10.5){
635 Double_t x,alpha,param[5],cov[15];
636 //
c0b978f0 637 fT1.fESDTrack.GetInnerExternalParameters(alpha,x,param);
51ad6848 638 fT1.fESDTrack.GetInnerExternalCovariance(cov);
51ad6848 639 AliExternalTrackParam paramm(x,alpha,param,cov);
640 //
c0b978f0 641 fT2.fESDTrack.GetInnerExternalParameters(alpha,x,param);
51ad6848 642 fT2.fESDTrack.GetInnerExternalCovariance(cov);
51ad6848 643 AliExternalTrackParam paramd(x,alpha,param,cov);
644 }
645 //
646 //
647
648 Float_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]};
649 Float_t p[3] = {fPdr[0]+fPm[0], fPdr[1]+fPm[1],fPdr[2]+fPm[2]};
650
651 Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
652 Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
653 vnorm2 = TMath::Sqrt(vnorm2);
654 Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
655 Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
656 pnorm2 = TMath::Sqrt(pnorm2);
657
658 fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
659 fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);
660 fPointAngle = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
661 }
662}
663
664////
665void AliESDRecKinkInfo::Update()
ae7d73d2 666{
51ad6848 667
668 if ( (fT1.fTPCOn)&& (fT2.fTPCOn)){
669 //
670 // IF BOTH RECONSTRUCTED
671 Float_t distance1,distance2;
672 Double_t xx[3],pp[3];
673 //
674 Double_t xd[3],pd[3],signd;
675 Double_t xm[3],pm[3],signm;
676 for (Int_t i=0;i<3;i++){
677 xd[i] = fT2.fTPCinR1[i];
678 pd[i] = fT2.fTPCinP1[i];
679 xm[i] = fT1.fTPCinR1[i];
680 pm[i] = fT1.fTPCinP1[i];
681 }
682 signd = fT2.fSign<0 ? -1:1;
683 signm = fT1.fSign<0 ? -1:1;
684
685 AliHelix dhelix1(xd,pd,signd);
686 dhelix1.GetMomentum(0,pp,0);
687 dhelix1.Evaluate(0,xx);
688 //
689 // Double_t x2[3],p2[3];
690 //
691 AliHelix mhelix(xm,pm,signm);
692 //
693 //find intersection linear
694 //
695 Double_t phase[2][2],radius[2];
696 Int_t points = dhelix1.GetRPHIintersections(mhelix, phase, radius,200);
697 Double_t delta1=10000,delta2=10000;
698
699 if (points==1){
700 fMinR = TMath::Sqrt(radius[0]);
701 }
702 if (points==2){
703 fMinR =TMath::Min(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1]));
704 }
705
706 if (points>0){
707 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
708 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
709 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
710 }
711 if (points==2){
712 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
713 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
714 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
715 }
716 if (points==1){
717 fMinR = TMath::Sqrt(radius[0]);
718 fDistMinR = delta1;
719 }
720 if (points==2){
721 if (radius[0]<radius[1]){
722 fMinR = TMath::Sqrt(radius[0]);
723 fDistMinR = delta1;
724 }
725 else{
726 fMinR = TMath::Sqrt(radius[1]);
727 fDistMinR = delta2;
728 }
729 }
730 //
731 //
732 distance1 = TMath::Min(delta1,delta2);
733 //
734 //find intersection parabolic
735 //
736 points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
737 delta1=10000,delta2=10000;
738
739 if (points>0){
740 dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
741 }
742 if (points==2){
743 dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
744 }
745
746 distance2 = TMath::Min(delta1,delta2);
747 if (delta1<delta2){
748 //get V0 info
749 dhelix1.Evaluate(phase[0][0],fXr);
750 dhelix1.GetMomentum(phase[0][0],fPdr);
751 mhelix.GetMomentum(phase[0][1],fPm);
752 dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle);
753 fRr = TMath::Sqrt(radius[0]);
754 }
755 else{
756 dhelix1.Evaluate(phase[1][0],fXr);
757 dhelix1.GetMomentum(phase[1][0], fPdr);
758 mhelix.GetMomentum(phase[1][1], fPm);
759 dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle);
760 fRr = TMath::Sqrt(radius[1]);
761 }
762 fDist1 = TMath::Sqrt(distance1);
763 fDist2 = TMath::Sqrt(distance2);
764
765 if (fDist2<10.5){
766 Double_t x,alpha,param[5],cov[15];
767 //
c0b978f0 768 fT1.fESDTrack.GetInnerExternalParameters(alpha,x,param);
51ad6848 769 fT1.fESDTrack.GetInnerExternalCovariance(cov);
51ad6848 770 AliExternalTrackParam paramm(x,alpha,param,cov);
771 //
c0b978f0 772 fT2.fESDTrack.GetInnerExternalParameters(alpha,x,param);
51ad6848 773 fT2.fESDTrack.GetInnerExternalCovariance(cov);
51ad6848 774 AliExternalTrackParam paramd(x,alpha,param,cov);
775 /*
776 AliESDkink kink;
777 kink.Update(&paramm,&paramd);
778 // kink.Dump();
779 Double_t diff = kink.fRr-fRr;
780 Double_t diff2 = kink.fDist2-fDist2;
781 printf("Diff\t%f\t%f\n",diff,diff2);
782 */
783 }
784
785 //
786 //
787 }
ae7d73d2 788
789}
790
791
792////////////////////////////////////////////////////////////////////////
793ESDCmpTr::ESDCmpTr()
794{
795 Reset();
796}
797
798////////////////////////////////////////////////////////////////////////
799ESDCmpTr::ESDCmpTr(const char* fnGenTracks,
800 const char* fnCmp,
801 const char* fnGalice, Int_t direction,
802 Int_t nEvents, Int_t firstEvent)
803{
804 Reset();
805 // fFnGenTracks = fnGenTracks;
806 // fFnCmp = fnCmp;
807 sprintf(fFnGenTracks,"%s",fnGenTracks);
808 sprintf(fFnCmp,"%s",fnCmp);
809
810 fFirstEventNr = firstEvent;
811 fEventNr = firstEvent;
812 fNEvents = nEvents;
813 fDirection = direction;
814 //
815 fLoader = AliRunLoader::Open(fnGalice);
816 if (gAlice){
33c3c91a 817 //delete AliRunLoader::Instance();
ae7d73d2 818 delete gAlice;
819 gAlice = 0x0;
820 }
821 if (fLoader->LoadgAlice()){
822 cerr<<"Error occured while l"<<endl;
823 }
824 Int_t nall = fLoader->GetNumberOfEvents();
825 if (nEvents==0) {
826 nEvents =nall;
827 fNEvents=nall;
828 fFirstEventNr=0;
829 }
830
831 if (nall<=0){
832 cerr<<"no events available"<<endl;
833 fEventNr = 0;
834 return;
835 }
836 if (firstEvent+nEvents>nall) {
837 fEventNr = nall-firstEvent;
838 cerr<<"restricted number of events availaible"<<endl;
839 }
ae7d73d2 840
841}
842
843
844////////////////////////////////////////////////////////////////////////
845ESDCmpTr::~ESDCmpTr()
846{
847 if (fLoader) {
848 delete fLoader;
849 }
850}
851
852//////////////////////////////////////////////////////////////
853Int_t ESDCmpTr::SetIO()
854{
855 //
856 //
857 CreateTreeCmp();
858 if (!fTreeCmp) return 1;
859 fParamTPC = GetTPCParam();
860 //
861 if (!ConnectGenTree()) {
862 cerr<<"Cannot connect tree with generated tracks"<<endl;
863 return 1;
864 }
865 return 0;
866}
867
868//////////////////////////////////////////////////////////////
869
870Int_t ESDCmpTr::SetIO(Int_t eventNr)
871{
872 //
873 //
874 // SET INPUT
875 //
876 TFile f("AliESDs.root");
877 //
878
879 TTree* tree = (TTree*) f.Get("esdTree");
86ad5fcb 880 if(fEvent)delete fEvent;
881 fEvent = new AliESDEvent();
882 fEvent->ReadFromTree(tree);
883 tree->GetEntry(eventNr);
884
ae7d73d2 885
886
887 /*
888 Char_t ename[100];
889 sprintf(ename,"%d",eventNr);
890 fEvent = (AliESD*)f.Get(ename);
891 if (!fEvent){
892 sprintf(ename,"ESD%d",eventNr);
893 fEvent = (AliESD*)f.Get(ename);
894 }
895
896 TTree* tree = (TTree*) f.Get("esdTree");
897 if (!tree) {
898 Error("CheckESD", "no ESD tree found");
899 return kFALSE;
900 }
901 tree->SetBranchAddress("ESD", &fEvent);
902 tree->GetEntry(eventNr);
903 */
904
905 if (!fEvent) return 1;
906
907 return 0;
908}
909
910
911
912////////////////////////////////////////////////////////////////////////
913void ESDCmpTr::Reset()
914{
915 fEventNr = 0;
916 fNEvents = 0;
917 fTreeCmp = 0;
51ad6848 918 fTreeCmpKinks =0;
919 fTreeCmpV0 =0;
ae7d73d2 920 // fFnCmp = "cmpTracks.root";
921 fFileGenTracks = 0;
922 fDebug = 0;
923 //
924 fParamTPC = 0;
925 fEvent =0;
926}
927
928////////////////////////////////////////////////////////////////////////
929Int_t ESDCmpTr::Exec(Int_t nEvents, Int_t firstEventNr)
930{
931 fNEvents = nEvents;
932 fFirstEventNr = firstEventNr;
933 return Exec();
934}
935
936////////////////////////////////////////////////////////////////////////
937Int_t ESDCmpTr::Exec()
938{
939 TStopwatch timer;
940 timer.Start();
941
942 if (SetIO()==1)
943 return 1;
944
945 fNextTreeGenEntryToRead = 0;
51ad6848 946 fNextKinkToRead = 0;
81e97e0d 947 fNextV0ToRead =0;
ae7d73d2 948 cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
949 for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
950 eventNr++) {
951 fEventNr = eventNr;
952 SetIO(fEventNr);
953 fNParticles = gAlice->GetEvent(fEventNr);
954
f007cb5f 955 fIndexRecTracks = new Short_t[fNParticles*20]; //write at maximum 4 tracks corresponding to particle
956 fIndexRecKinks = new Short_t[fNParticles*20]; //write at maximum 20 tracks corresponding to particle
957 fIndexRecV0 = new Short_t[fNParticles*20]; //write at maximum 20 tracks corresponding to particle
51ad6848 958
f007cb5f 959 fFakeRecTracks = new Short_t[fNParticles];
960 fMultiRecTracks = new Short_t[fNParticles];
961 fMultiRecKinks = new Short_t[fNParticles];
962 fMultiRecV0 = new Short_t[fNParticles];
ae7d73d2 963
51ad6848 964 for (Int_t i = 0; i<fNParticles; i++) {
965 for (Int_t j=0;j<20;j++){
966 fIndexRecTracks[20*i+j] = -1;
967 fIndexRecKinks[20*i+j] = -1;
968 fIndexRecV0[20*i+j] = -1;
969 }
ae7d73d2 970 fFakeRecTracks[i] = 0;
971 fMultiRecTracks[i] = 0;
51ad6848 972 fMultiRecKinks[i] = 0;
973 fMultiRecV0[i] = 0;
ae7d73d2 974 }
975
976 cout<<"Start to process event "<<fEventNr<<endl;
977 cout<<"\tfNParticles = "<<fNParticles<<endl;
978 if (fDebug>2) cout<<"\tStart loop over TreeT"<<endl;
979 if (TreeTLoop()>0) return 1;
980
981 if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl;
982 if (TreeGenLoop(eventNr)>0) return 1;
51ad6848 983 BuildKinkInfo0(eventNr);
984 BuildV0Info(eventNr);
985 fRecArray->Delete();
986
ae7d73d2 987 if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl;
988
989 delete [] fIndexRecTracks;
51ad6848 990 delete [] fIndexRecKinks;
991 delete [] fIndexRecV0;
ae7d73d2 992 delete [] fFakeRecTracks;
993 delete [] fMultiRecTracks;
51ad6848 994 delete [] fMultiRecKinks;
995 delete [] fMultiRecV0;
ae7d73d2 996 }
997
998 CloseOutputFile();
999
1000 cerr<<"Exec finished"<<endl;
1001 timer.Stop();
1002 timer.Print();
1003 return 0;
1004
1005}
1006////////////////////////////////////////////////////////////////////////
1007Bool_t ESDCmpTr::ConnectGenTree()
1008{
1009//
1010// connect all branches from the genTracksTree
1011// use the same variables as for the new cmp tree, it may work
1012//
1013 fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
1014 if (!fFileGenTracks) {
1015 cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
1016 return kFALSE;
1017 }
1018 fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
1019 if (!fTreeGenTracks) {
1020 cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
1021 <<fFnGenTracks<<endl;
1022 return kFALSE;
1023 }
1024 //
1025 fMCInfo = new AliMCInfo;
1026 fTreeGenTracks->SetBranchAddress("MC",&fMCInfo);
1027 //
51ad6848 1028 //
1029 fTreeGenKinks = (TTree*)fFileGenTracks->Get("genKinksTree");
1030 if (!fTreeGenKinks) {
1031 cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
1032 <<fFnGenTracks<<endl;
1033 //return kFALSE;
1034 }
1035 else{
1036 fGenKinkInfo = new AliGenKinkInfo;
1037 fTreeGenKinks->SetBranchAddress("MC",&fGenKinkInfo);
1038 }
1039
1040 fTreeGenV0 = (TTree*)fFileGenTracks->Get("genV0Tree");
1041 if (!fTreeGenV0) {
1042 cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
1043 <<fFnGenTracks<<endl;
1044 //return kFALSE;
1045 }
1046 else{
1047 fGenV0Info = new AliGenV0Info;
1048 fTreeGenV0->SetBranchAddress("MC",&fGenV0Info);
1049 }
1050 //
ae7d73d2 1051 if (fDebug > 1) {
1052 cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
1053 }
1054 return kTRUE;
1055}
1056
1057
1058////////////////////////////////////////////////////////////////////////
1059void ESDCmpTr::CreateTreeCmp()
1060{
1061 fFileCmp = TFile::Open(fFnCmp,"RECREATE");
1062 if (!fFileCmp) {
1063 cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
1064 return;
1065 }
ae7d73d2 1066 //
51ad6848 1067 //
1068 fTreeCmp = new TTree("ESDcmpTracks","ESDcmpTracks");
ae7d73d2 1069 fMCInfo = new AliMCInfo;
1070 fRecInfo = new AliESDRecInfo;
ae7d73d2 1071 AliESDtrack * esdTrack = new AliESDtrack;
51ad6848 1072 // AliITStrackMI * itsTrack = new AliITStrackMI;
e1eaa672 1073 fTreeCmp->Branch("MC","AliMCInfo",&fMCInfo,256000);
1074 fTreeCmp->Branch("RC","AliESDRecInfo",&fRecInfo,256000);
f007cb5f 1075 // fTreeCmp->Branch("fESDTrack","AliESDtrack",&esdTrack);
51ad6848 1076 // fTreeCmp->Branch("ITS","AliITStrackMI",&itsTrack);
ae7d73d2 1077 delete esdTrack;
1078 //
51ad6848 1079 //
1080 fTreeCmpKinks = new TTree("ESDcmpKinks","ESDcmpKinks");
1081 fGenKinkInfo = new AliGenKinkInfo;
1082 fRecKinkInfo = new AliESDRecKinkInfo;
e1eaa672 1083 fTreeCmpKinks->Branch("MC.","AliGenKinkInfo",&fGenKinkInfo,256000);
1084 fTreeCmpKinks->Branch("RC.","AliESDRecKinkInfo",&fRecKinkInfo,256000);
51ad6848 1085 //
1086 //
1087 fTreeCmpV0 = new TTree("ESDcmpV0","ESDcmpV0");
1088 fGenV0Info = new AliGenV0Info;
1089 fRecV0Info = new AliESDRecV0Info;
e1eaa672 1090 fTreeCmpV0->Branch("MC.","AliGenV0Info", &fGenV0Info,256000);
1091 fTreeCmpV0->Branch("RC.","AliESDRecV0Info",&fRecV0Info,256000);
51ad6848 1092 //
ae7d73d2 1093 fTreeCmp->AutoSave();
51ad6848 1094 fTreeCmpKinks->AutoSave();
1095 fTreeCmpV0->AutoSave();
ae7d73d2 1096}
1097////////////////////////////////////////////////////////////////////////
1098void ESDCmpTr::CloseOutputFile()
1099{
1100 if (!fFileCmp) {
1101 cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
1102 return;
1103 }
1104 fFileCmp->cd();
1105 fTreeCmp->Write();
1106 delete fTreeCmp;
1107
1108 fFileCmp->Close();
1109 delete fFileCmp;
1110 return;
1111}
1112////////////////////////////////////////////////////////////////////////
1113
1114TVector3 ESDCmpTr::TR2Local(AliTrackReference *trackRef,
1115 AliTPCParam *paramTPC) {
1116
1117 Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
1118 Int_t index[4];
1119 paramTPC->Transform0to1(x,index);
9502e15e 1120 paramTPC->Transform1to2Ideal(x,index);
ae7d73d2 1121 return TVector3(x);
1122}
1123////////////////////////////////////////////////////////////////////////
1124
1125Int_t ESDCmpTr::TreeTLoop()
1126{
1127 //
1128 // loop over all ESD reconstructed tracks and store info in memory
1129 //
51ad6848 1130 // + loop over all reconstructed kinks
ae7d73d2 1131 TStopwatch timer;
1132 timer.Start();
1133 //
1134 Int_t nEntries = (Int_t)fEvent->GetNumberOfTracks();
51ad6848 1135 Int_t nKinks = (Int_t) fEvent->GetNumberOfKinks();
d6a49f20 1136 Int_t nV0MIs = (Int_t) fEvent->GetNumberOfV0s();
f007cb5f 1137 fSignedKinks = new Short_t[nKinks];
1138 fSignedV0 = new Short_t[nV0MIs];
ae7d73d2 1139 //
51ad6848 1140 // load kinks to the memory
1141 for (Int_t i=0; i<nKinks;i++){
1142 AliESDkink * kink =fEvent->GetKink(i);
1143 fSignedKinks[i]=0;
1144 if (kink->fStatus<0) continue;
1145 }
ae7d73d2 1146 //
51ad6848 1147 for (Int_t i=0; i<nV0MIs;i++){
d6a49f20 1148 AliESDv0 * v0MI =fEvent->GetV0(i);
51ad6848 1149 fSignedV0[i]=0;
1150 if (v0MI->fStatus<0) continue;
ae7d73d2 1151 }
51ad6848 1152
ae7d73d2 1153 //
1154 //
1155 AliESDtrack * track=0;
1156 for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
51ad6848 1157 //track = (AliESDtrack*)fTracks->UncheckedAt(iEntry);
1158 track = (AliESDtrack*)fEvent->GetTrack(iEntry);
ae7d73d2 1159 //
1160 Int_t label = track->GetLabel();
1161 Int_t absLabel = abs(label);
1162 if (absLabel < fNParticles) {
1163 // fIndexRecTracks[absLabel] = iEntry;
1164 if (label < 0) fFakeRecTracks[absLabel]++;
1165 if (fMultiRecTracks[absLabel]>0){
51ad6848 1166 if (fMultiRecTracks[absLabel]<20)
1167 fIndexRecTracks[absLabel*20+fMultiRecTracks[absLabel]] = iEntry;
ae7d73d2 1168 }
1169 else
51ad6848 1170 fIndexRecTracks[absLabel*20] = iEntry;
ae7d73d2 1171 fMultiRecTracks[absLabel]++;
1172 }
51ad6848 1173 }
1174 // sort reconstructed kinks
1175 //
1176 AliESDkink * kink=0;
1177 for (Int_t iEntry=0; iEntry<nKinks;iEntry++){
1178 kink = (AliESDkink*)fEvent->GetKink(iEntry);
1179 if (!kink) continue;
1180 //
1181 Int_t label0 = TMath::Abs(kink->fLab[0]);
1182 Int_t label1 = TMath::Abs(kink->fLab[1]);
51ad6848 1183 Int_t absLabel = TMath::Min(label0,label1);
1184 if (absLabel < fNParticles) {
1185 if (fMultiRecKinks[absLabel]>0){
1186 if (fMultiRecKinks[absLabel]<20)
1187 fIndexRecKinks[absLabel*20+fMultiRecKinks[absLabel]] = iEntry;
1188 }
1189 else
1190 fIndexRecKinks[absLabel*20] = iEntry;
1191 fMultiRecKinks[absLabel]++;
1192 }
ae7d73d2 1193 }
51ad6848 1194 // --sort reconstructed V0
1195 //
d6a49f20 1196 AliESDv0 * v0MI=0;
51ad6848 1197 for (Int_t iEntry=0; iEntry<nV0MIs;iEntry++){
d6a49f20 1198 v0MI = fEvent->GetV0(iEntry);
51ad6848 1199 if (!v0MI) continue;
1200 //
f007cb5f 1201 // Int_t label0 = TMath::Abs(v0MI->fLab[0]);
1202 //Int_t label1 = TMath::Abs(v0MI->fLab[1]);
1203 //
1204 for (Int_t i=0;i<2;i++){
1205 // Int_t absLabel = TMath::Min(label0,label1);
1206 Int_t absLabel = TMath::Abs(v0MI->fLab[i]);
1207 if (absLabel < fNParticles) {
1208 if (fMultiRecV0[absLabel]>0){
1209 if (fMultiRecV0[absLabel]<20)
1210 fIndexRecV0[absLabel*20+fMultiRecV0[absLabel]] = iEntry;
1211 }
1212 else
1213 fIndexRecV0[absLabel*20] = iEntry;
1214 fMultiRecV0[absLabel]++;
51ad6848 1215 }
51ad6848 1216 }
1217 }
1218
1219
ae7d73d2 1220 printf("Time spended in TreeTLoop\n");
1221 timer.Print();
1222
1223 if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;
1224 return 0;
1225}
1226
1227////////////////////////////////////////////////////////////////////////
1228Int_t ESDCmpTr::TreeGenLoop(Int_t eventNr)
1229{
1230//
1231// loop over all entries for a given event, find corresponding
1232// rec. track and store in the fTreeCmp
1233//
1234 TStopwatch timer;
1235 timer.Start();
1236 Int_t entry = fNextTreeGenEntryToRead;
1237 Double_t nParticlesTR = fTreeGenTracks->GetEntriesFast();
1238 cerr<<"fNParticles, nParticlesTR, fNextTreeGenEntryToRead: "<<fNParticles<<" "
1239 <<nParticlesTR<<" "<<fNextTreeGenEntryToRead<<endl;
1240 TBranch * branch = fTreeCmp->GetBranch("RC");
1241 branch->SetAddress(&fRecInfo); // set all pointers
51ad6848 1242 fRecArray = new TObjArray(fNParticles);
98808dce 1243 AliESDtrack dummytrack; //
ae7d73d2 1244
1245 while (entry < nParticlesTR) {
1246 fTreeGenTracks->GetEntry(entry);
1247 entry++;
1248 if (eventNr < fMCInfo->fEventNr) continue;
1249 if (eventNr > fMCInfo->fEventNr) continue;;
1250 //
1251 fNextTreeGenEntryToRead = entry-1;
1252 if (fDebug > 2 && fMCInfo->fLabel < 10) {
1253 cerr<<"Fill track with a label "<<fMCInfo->fLabel<<endl;
1254 }
51ad6848 1255 // if (fMCInfo->fNTPCRef<1) continue; // not TPCref
ae7d73d2 1256 //
1257 fRecInfo->Reset();
1258 AliESDtrack * track=0;
1259 fRecInfo->fReconstructed =0;
1260 TVector3 local = TR2Local(&(fMCInfo->fTrackRef),fParamTPC);
1261 local.GetXYZ(fRecInfo->fTRLocalCoord);
1262 //
51ad6848 1263 if (fIndexRecTracks[fMCInfo->fLabel*20] >= 0) {
1264 //track= (AliESDtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4]);
1265 track= (AliESDtrack*)fEvent->GetTrack(fIndexRecTracks[fMCInfo->fLabel*20]);
ae7d73d2 1266 //
1267 //
1268 // find nearest track if multifound
51ad6848 1269 //Int_t sign = Int_t(track->GetSign()*fMCInfo->fCharge);
1270 //
1271 Int_t status = 0;
1272 if ((track->GetStatus()&AliESDtrack::kITSrefit)>0) status++;
1273 if ((track->GetStatus()&AliESDtrack::kTPCrefit)>0) status++;
1274 if ((track->GetStatus()&AliESDtrack::kTRDrefit)>0) status++;
1275
1276 //
1277 if (fIndexRecTracks[fMCInfo->fLabel*20+1]>0){
1278 //
1279 Double_t p[3];
1280 track->GetInnerPxPyPz(p);
1281 Float_t maxp = p[0]*p[0]+p[1]*p[1]+p[2]*p[2];
1282 //
1283 for (Int_t i=1;i<20;i++){
1284 if (fIndexRecTracks[fMCInfo->fLabel*20+i]>=0){
1285 AliESDtrack * track2 = (AliESDtrack*)fEvent->GetTrack(fIndexRecTracks[fMCInfo->fLabel*20+i]);
1286 if (!track2) continue;
1287 //Int_t sign2 = track->GetSign()*fMCInfo->fCharge; //
1288 //if (sign2<0) continue;
1289 track2->GetInnerPxPyPz(p);
1290 Float_t mom = p[0]*p[0]+p[1]*p[1]+p[2]*p[2];
1291 /*
1292 if (sign<0){
1293 sign = sign2;
1294 track = track2;
1295 maxp = mom;
1296 continue;
1297 }
1298 */
1299 //
1300 Int_t status2 = 0;
1301 if ((track2->GetStatus()&AliESDtrack::kITSrefit)>0) status2++;
1302 if ((track2->GetStatus()&AliESDtrack::kTPCrefit)>0) status2++;
1303 if ((track2->GetStatus()&AliESDtrack::kTRDrefit)>0) status2++;
1304 if (status2<status) continue;
1305 //
1306 if (mom<maxp) continue;
1307 maxp = mom;
1308 track = track2;
1309 //
ae7d73d2 1310 }
1311 }
1312 }
1313 //
98808dce 1314 if (track) {
0a1e73ec 1315 new (&(fRecInfo->fESDTrack)) AliESDtrack(*track);
98808dce 1316 }else{
1317 fRecInfo->fESDTrack = dummytrack;
1318 }
1319
ae7d73d2 1320 if (track->GetITStrack())
51ad6848 1321 fRecInfo->fITStrack = *((AliITStrackMI*)track->GetITStrack());
81e97e0d 1322 else{
1323 fRecInfo->fITStrack = *track;
1324 }
ae7d73d2 1325 if (track->GetTRDtrack()){
1326 fRecInfo->fTRDtrack = *((AliTRDtrack*)track->GetTRDtrack());
1327 }
1328 else{
1329 fRecInfo->fTRDtrack.SetdEdx(-1);
1330 }
1331 fRecInfo->fReconstructed = 1;
1332 fRecInfo->fFake = fFakeRecTracks[fMCInfo->fLabel];
1333 fRecInfo->fMultiple = fMultiRecTracks[fMCInfo->fLabel];
1334 //
98808dce 1335 fRecInfo->Update(fMCInfo,fParamTPC,kTRUE, fEvent);
ae7d73d2 1336 }
1337 else{
98808dce 1338 fRecInfo->fESDTrack = dummytrack;
1339 fRecInfo->Update(fMCInfo,fParamTPC,kFALSE, fEvent);
ae7d73d2 1340 }
51ad6848 1341 fRecArray->AddAt(new AliESDRecInfo(*fRecInfo),fMCInfo->fLabel);
ae7d73d2 1342 fTreeCmp->Fill();
1343 }
1344 fTreeCmp->AutoSave();
51ad6848 1345 //fTracks->Delete();
ae7d73d2 1346 printf("Time spended in TreeGenLoop\n");
1347 timer.Print();
1348 if (fDebug > 2) cerr<<"end of TreeGenLoop"<<endl;
1349
1350 return 0;
1351}
51ad6848 1352
1353
1354
1355////////////////////////////////////////////////////////////////////////
1356////////////////////////////////////////////////////////////////////////
1357////////////////////////////////////////////////////////////////////////
1358Int_t ESDCmpTr::BuildKinkInfo0(Int_t eventNr)
1359{
1360//
1361// loop over all entries for a given event, find corresponding
1362// rec. track and store in the fTreeCmp
1363//
1364 TStopwatch timer;
1365 timer.Start();
1366 Int_t entry = fNextKinkToRead;
1367 Double_t nParticlesTR = fTreeGenKinks->GetEntriesFast();
1368 cerr<<"fNParticles, nParticlesTR, fNextKinkToRead: "<<fNParticles<<" "
1369 <<nParticlesTR<<" "<<fNextKinkToRead<<endl;
1370 //
1371 TBranch * branch = fTreeCmpKinks->GetBranch("RC.");
1372 branch->SetAddress(&fRecKinkInfo); // set all pointers
1373
1374 //
1375 while (entry < nParticlesTR) {
1376 fTreeGenKinks->GetEntry(entry);
1377 entry++;
1378 if (eventNr < fGenKinkInfo->fMCm.fEventNr) continue;
1379 if (eventNr > fGenKinkInfo->fMCm.fEventNr) continue;;
1380 //
1381 fNextKinkToRead = entry-1;
1382 //
1383 //
1384 AliESDRecInfo* fRecInfo1 = (AliESDRecInfo*)fRecArray->At(fGenKinkInfo->fMCm.fLabel);
1385 AliESDRecInfo* fRecInfo2 = (AliESDRecInfo*)fRecArray->At(fGenKinkInfo->fMCd.fLabel);
1386 fRecKinkInfo->fT1 = (*fRecInfo1);
1387 fRecKinkInfo->fT2 = (*fRecInfo2);
1388 fRecKinkInfo->fStatus =0;
1389 if (fRecInfo1 && fRecInfo1->fTPCOn) fRecKinkInfo->fStatus+=1;
1390 if (fRecInfo2 && fRecInfo2->fTPCOn) fRecKinkInfo->fStatus+=2;
1391 if (fRecKinkInfo->fStatus==3&&fRecInfo1->fSign!=fRecInfo2->fSign) fRecKinkInfo->fStatus*=-1;
1392
1393 if (fRecKinkInfo->fStatus==3){
1394 fRecKinkInfo->Update();
1395 }
1396 Int_t label = TMath::Min(fGenKinkInfo->fMCm.fLabel,fGenKinkInfo->fMCd.fLabel);
1397 Int_t label2 = TMath::Max(fGenKinkInfo->fMCm.fLabel,fGenKinkInfo->fMCd.fLabel);
1398
1399 AliESDkink *kink=0;
1400 fRecKinkInfo->fRecStatus =0;
1401 fRecKinkInfo->fMultiple = fMultiRecKinks[label];
1402 fRecKinkInfo->fKinkMultiple=0;
1403 //
1404 if (fMultiRecKinks[label]>0){
1405
1406 // for (Int_t j=0;j<TMath::Min(fMultiRecKinks[label],100);j++){
f007cb5f 1407 for (Int_t j=TMath::Min(fMultiRecKinks[label],Short_t(20))-1;j>=0;j--){
51ad6848 1408 Int_t index = fIndexRecKinks[label*20+j];
1409 //AliESDkink *kink2 = (AliESDkink*)fKinks->At(index);
1410 AliESDkink *kink2 = (AliESDkink*)fEvent->GetKink(index);
1411 if (TMath::Abs(kink2->fLab[0])==label &&TMath::Abs(kink2->fLab[1])==label2) {
51ad6848 1412 fRecKinkInfo->fKinkMultiple++;
1413 fSignedKinks[index]=1;
f007cb5f 1414 Int_t c0=0;
1415 if (kink){
1416 // if (kink->fTRDOn) c0++;
1417 //if (kink->fITSOn) c0++;
1418 if (kink->GetStatus(2)>0) c0++;
1419 if (kink->GetStatus(0)>0) c0++;
1420 }
1421 Int_t c2=0;
1422 // if (kink2->fTRDOn) c2++;
1423 //if (kink2->fITSOn) c2++;
1424 if (kink2->GetStatus(2)>0) c2++;
1425 if (kink2->GetStatus(0)>0) c2++;
1426
1427 if (c2<c0) continue;
1428 kink =kink2;
51ad6848 1429 }
1430 if (TMath::Abs(kink2->fLab[1])==label &&TMath::Abs(kink2->fLab[0])==label2) {
51ad6848 1431 fRecKinkInfo->fKinkMultiple++;
1432 fSignedKinks[index]=1;
f007cb5f 1433 Int_t c0=0;
1434 if (kink){
1435 //if (kink->fTRDOn) c0++;
1436 //if (kink->fITSOn) c0++;
1437 if (kink->GetStatus(2)>0) c0++;
1438 if (kink->GetStatus(0)>0) c0++;
1439
1440 }
1441 Int_t c2=0;
1442 // if (kink2->fTRDOn) c2++;
1443 //if (kink2->fITSOn) c2++;
1444 if (kink2->GetStatus(2)>0) c2++;
1445 if (kink2->GetStatus(0)>0) c2++;
1446
1447 if (c2<c0) continue;
1448 kink =kink2;
51ad6848 1449 }
1450 }
1451 }
1452 if (kink){
1453 fRecKinkInfo->fKink = *kink;
1454 fRecKinkInfo->fRecStatus=1;
1455 }
1456 fTreeCmpKinks->Fill();
1457 }
1458 // Int_t nkinks = fKinks->GetEntriesFast();
1459 Int_t nkinks = fEvent->GetNumberOfKinks();
1460 for (Int_t i=0;i<nkinks;i++){
1461 if (fSignedKinks[i]==0){
1462 // AliESDkink *kink = (AliESDkink*)fKinks->At(i);
1463 AliESDkink *kink = (AliESDkink*)fEvent->GetKink(i);
1464 if (!kink) continue;
1465 //
1466 fRecKinkInfo->fKink = *kink;
1467 fRecKinkInfo->fRecStatus =-2;
1468 //
1469 AliESDRecInfo* fRecInfo1 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(kink->fLab[0]));
1470 AliESDRecInfo* fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(kink->fLab[1]));
1471 if (fRecInfo1 && fRecInfo2){
1472 fRecKinkInfo->fT1 = (*fRecInfo1);
1473 fRecKinkInfo->fT2 = (*fRecInfo2);
1474 fRecKinkInfo->fRecStatus =-1;
1475 }
1476 fTreeCmpKinks->Fill();
1477 }
1478 }
1479
1480
1481 fTreeCmpKinks->AutoSave();
1482 printf("Time spended in BuilKinkInfo Loop\n");
1483 timer.Print();
1484 if (fDebug > 2) cerr<<"end of BuildKinkInfo Loop"<<endl;
1485 return 0;
1486}
1487
f007cb5f 1488
1489
1490void ESDCmpTr::MakePoints(AliESDtrack * track, AliPointsMI &points)
1491{
1492 //
1493 // make points in global coordinate frame
1494 //
1495 return;
1496 /*
1497 UInt_t itscl[10];
1498 Int_t nits = track->GetITSclusters(itscl);
1499 Int_t tpccl[1000];
1500 Int_t ntpc = track->GetTPcclusters(tpccl);
1501 UInt_t trdcl[1000];
1502 Int_t ntrd = track->GetTRDclusters(trdcl);
1503 //
1504 AliLoader *itsloader = fLoader->GetLoader("ITSLoader");
1505 AliLoader *tpcloader = fLoader->GetLoader("TPCLoader");
1506 AliLoader *trdloader = fLoader->GetLoader("TRDLoader");
1507 //
1508 AliITStrackerMI itstracker();
1509 */
1510}
1511
1512
51ad6848 1513////////////////////////////////////////////////////////////////////////
1514////////////////////////////////////////////////////////////////////////
1515////////////////////////////////////////////////////////////////////////
1516
f007cb5f 1517
1518
51ad6848 1519Int_t ESDCmpTr::BuildV0Info(Int_t eventNr)
1520{
1521//
1522// loop over all entries for a given event, find corresponding
1523// rec. track and store in the fTreeCmp
1524//
1525 TStopwatch timer;
1526 timer.Start();
1527 Int_t entry = fNextV0ToRead;
1528 Double_t nParticlesTR = fTreeGenV0->GetEntriesFast();
1529 cerr<<"fNParticles, nParticlesTR, fNextV0ToRead: "<<fNParticles<<" "
1530 <<nParticlesTR<<" "<<fNextV0ToRead<<endl;
1531 //
1532 TBranch * branch = fTreeCmpV0->GetBranch("RC.");
1533 branch->SetAddress(&fRecV0Info); // set all pointers
1534 const AliESDVertex * esdvertex = fEvent->GetVertex();
1535 Float_t vertex[3]= {esdvertex->GetXv(), esdvertex->GetYv(),esdvertex->GetZv()};
1536
1537 //
1538 while (entry < nParticlesTR) {
1539 fTreeGenV0->GetEntry(entry);
1540 entry++;
1541 if (eventNr < fGenV0Info->fMCm.fEventNr) continue;
1542 if (eventNr > fGenV0Info->fMCm.fEventNr) continue;;
1543 //
1544 fNextV0ToRead = entry-1;
1545 //
1546 //
1547 AliESDRecInfo* fRecInfo1 = (AliESDRecInfo*)fRecArray->At(fGenV0Info->fMCm.fLabel);
1548 AliESDRecInfo* fRecInfo2 = (AliESDRecInfo*)fRecArray->At(fGenV0Info->fMCd.fLabel);
1549 if (fGenV0Info->fMCm.fCharge*fGenV0Info->fMCd.fCharge>0) continue; // interactions
1550 if (!fRecInfo1 || !fRecInfo2) continue;
1551 fRecV0Info->fT1 = (*fRecInfo1);
1552 fRecV0Info->fT2 = (*fRecInfo2);
1553 fRecV0Info->fV0Status =0;
1554 if (fRecInfo1 && fRecInfo1->fStatus[1]>0) fRecV0Info->fV0Status+=1;
1555 if (fRecInfo2 && fRecInfo2->fStatus[1]>0) fRecV0Info->fV0Status+=2;
1556
1557 if (fRecV0Info->fV0Status==3&&fRecInfo1->fSign==fRecInfo2->fSign) fRecV0Info->fV0Status*=-1;
1558
1559
1560 if (abs(fRecV0Info->fV0Status)==3){
1561 fRecV0Info->Update(vertex);
1562 {
1563 //
1564 // TPC V0 Info
1565 Double_t x,alpha, param[5],cov[15];
0a1e73ec 1566 if ( fRecV0Info->fT1.fESDTrack.GetInnerParam() && fRecV0Info->fT2.fESDTrack.GetInnerParam()){
1567 fRecV0Info->fT1.fESDTrack.GetInnerExternalParameters(alpha,x,param);
1568 fRecV0Info->fT1.fESDTrack.GetInnerExternalCovariance(cov);
1569 AliExternalTrackParam paramP(x,alpha,param,cov);
1570 //
1571 fRecV0Info->fT2.fESDTrack.GetInnerExternalParameters(alpha,x,param);
1572 fRecV0Info->fT2.fESDTrack.GetInnerExternalCovariance(cov);
1573 AliExternalTrackParam paramM(x,alpha,param,cov);
1574 //
1575 fRecV0Info->fV0tpc.SetM(paramM);
1576 fRecV0Info->fV0tpc.SetP(paramP);
1577 Double_t pid1[5],pid2[5];
1578 fRecV0Info->fT1.fESDTrack.GetESDpid(pid1);
1579 fRecV0Info->fT1.fESDTrack.GetESDpid(pid2);
1580 //
1581 fRecV0Info->fV0tpc.UpdatePID(pid1,pid2);
1582 fRecV0Info->fV0tpc.Update(vertex);
1583
1584 //
1585 //
1586 fRecV0Info->fT1.fESDTrack.GetExternalParameters(x,param);
1587 fRecV0Info->fT1.fESDTrack.GetExternalCovariance(cov);
1588 alpha = fRecV0Info->fT1.fESDTrack.GetAlpha();
1589 new (&paramP) AliExternalTrackParam(x,alpha,param,cov);
1590 //
1591 fRecV0Info->fT2.fESDTrack.GetExternalParameters(x,param);
1592 fRecV0Info->fT2.fESDTrack.GetExternalCovariance(cov);
1593 alpha = fRecV0Info->fT2.fESDTrack.GetAlpha();
1594 new (&paramM) AliExternalTrackParam(x,alpha,param,cov);
1595 //
1596 fRecV0Info->fV0its.SetM(paramM);
1597 fRecV0Info->fV0its.SetP(paramP);
1598 fRecV0Info->fV0its.UpdatePID(pid1,pid2);
1599 fRecV0Info->fV0its.Update(vertex);
1600 }
51ad6848 1601 }
1602 if (TMath::Abs(fGenV0Info->fMCm.fPdg)==11 &&TMath::Abs(fGenV0Info->fMCd.fPdg)==11){
1603 if (fRecV0Info->fDist2>10){
1604 fRecV0Info->Update(vertex);
1605 }
1606 if (fRecV0Info->fDist2>10){
1607 fRecV0Info->Update(vertex);
1608 }
1609 }
1610 }
1611 //
1612 // take the V0 from reconstruction
1613
1614 Int_t label = TMath::Min(fGenV0Info->fMCm.fLabel,fGenV0Info->fMCd.fLabel);
1615 Int_t label2 = TMath::Max(fGenV0Info->fMCm.fLabel,fGenV0Info->fMCd.fLabel);
d6a49f20 1616 AliESDv0 *v0MI=0;
51ad6848 1617 fRecV0Info->fRecStatus =0;
1618 fRecV0Info->fMultiple = fMultiRecV0[label];
1619 fRecV0Info->fV0Multiple=0;
1620 //
f007cb5f 1621 if (fMultiRecV0[label]>0 || fMultiRecV0[label2]>0){
51ad6848 1622
1623 // for (Int_t j=0;j<TMath::Min(fMultiRecV0s[label],100);j++){
f007cb5f 1624 for (Int_t j=TMath::Min(fMultiRecV0[label],Short_t(20))-1;j>=0;j--){
51ad6848 1625 Int_t index = fIndexRecV0[label*20+j];
e1eaa672 1626 if (index<0) continue;
d6a49f20 1627 AliESDv0 *v0MI2 = fEvent->GetV0(index);
51ad6848 1628 if (TMath::Abs(v0MI2->fLab[0])==label &&TMath::Abs(v0MI2->fLab[1])==label2) {
1629 v0MI =v0MI2;
1630 fRecV0Info->fV0Multiple++;
1631 fSignedV0[index]=1;
1632 }
1633 if (TMath::Abs(v0MI2->fLab[1])==label &&TMath::Abs(v0MI2->fLab[0])==label2) {
1634 v0MI =v0MI2;
1635 fRecV0Info->fV0Multiple++;
1636 fSignedV0[index]=1;
1637 }
1638 }
1639 }
1640 if (v0MI){
1641 fRecV0Info->fV0rec = *v0MI;
1642 fRecV0Info->fRecStatus=1;
1643 }
1644
1645 fTreeCmpV0->Fill();
1646 }
1647 //
1648 // write fake v0s
1649
d6a49f20 1650 Int_t nV0MIs = fEvent->GetNumberOfV0s();
51ad6848 1651 for (Int_t i=0;i<nV0MIs;i++){
1652 if (fSignedV0[i]==0){
d6a49f20 1653 AliESDv0 *v0MI = fEvent->GetV0(i);
51ad6848 1654 if (!v0MI) continue;
1655 //
1656 fRecV0Info->fV0rec = *v0MI;
1657 fRecV0Info->fV0Status =-10;
1658 fRecV0Info->fRecStatus =-2;
1659 //
1660 AliESDRecInfo* fRecInfo1 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(v0MI->fLab[0]));
1661 AliESDRecInfo* fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(v0MI->fLab[1]));
1662 if (fRecInfo1 && fRecInfo2){
1663 fRecV0Info->fT1 = (*fRecInfo1);
1664 fRecV0Info->fT2 = (*fRecInfo2);
1665 fRecV0Info->fRecStatus =-1;
1666 }
f007cb5f 1667 fRecV0Info->Update(vertex);
51ad6848 1668 fTreeCmpV0->Fill();
1669 }
1670 }
1671
1672
1673
1674 fTreeCmpV0->AutoSave();
51ad6848 1675 printf("Time spended in BuilV0Info Loop\n");
1676 timer.Print();
1677 if (fDebug > 2) cerr<<"end of BuildV0Info Loop"<<endl;
1678 return 0;
1679}
ae7d73d2 1680////////////////////////////////////////////////////////////////////////
1681////////////////////////////////////////////////////////////////////////
1682
1683void AliESDComparisonDraw::SetIO(const char *fname)
1684{
1685 //
1686 TFile* file = TFile::Open(fname);
1687 //
1688 fTree = (TTree*) file->Get("ESDcmpTracks");
1689 if (!fTree) {
1690 printf("no track comparison tree found\n");
1691 file->Close();
1692 delete file;
1693 }
1694}
1695
1696