]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliGenInfo.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Changes for the kink finder. Coding conventions (M.Ivanov)
[u/mrichter/AliRoot.git] / STEER / AliGenInfo.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
20Origin: marian.ivanov@cern.ch
21
22Macro to generate comples MC information - used for Comparison later on
23How to use it?
24
25.L $ALICE_ROOT/STEER/AliGenInfo.C+
26AliGenInfoMaker *t = new AliGenInfoMaker("galice.root","genTracks.root",1,0)
27t->Exec();
28
29*/
30
31#if !defined(__CINT__) || defined(__MAKECINT__)
32#include <stdio.h>
33#include <string.h>
34//ROOT includes
35#include "Rtypes.h"
36#include "TFile.h"
37#include "TTree.h"
38#include "TChain.h"
39#include "TCut.h"
40#include "TString.h"
41#include "TBenchmark.h"
42#include "TStopwatch.h"
43#include "TParticle.h"
44#include "TSystem.h"
45#include "TTimer.h"
46#include "TVector3.h"
47#include "TH1F.h"
48#include "TH2F.h"
49#include "TCanvas.h"
50#include "TPad.h"
51#include "TF1.h"
51ad6848 52#include "TView.h"
53#include "TGeometry.h"
54#include "TPolyLine3D.h"
ae7d73d2 55
56//ALIROOT includes
57#include "AliRun.h"
58#include "AliStack.h"
59#include "AliSimDigits.h"
60#include "AliTPCParam.h"
61#include "AliTPC.h"
62#include "AliTPCLoader.h"
63#include "AliDetector.h"
64#include "AliTrackReference.h"
65#include "AliTPCParamSR.h"
66#include "AliTracker.h"
67#include "AliMagF.h"
51ad6848 68#include "AliHelix.h"
69#include "AliPoints.h"
70
ae7d73d2 71#endif
72#include "AliGenInfo.h"
73//
51ad6848 74//
ae7d73d2 75
76AliTPCParam * GetTPCParam(){
77 AliTPCParamSR * par = new AliTPCParamSR;
78 par->Update();
79 return par;
80}
81
82
83//_____________________________________________________________________________
84Float_t TPCBetheBloch(Float_t bg)
85{
86 //
87 // Bethe-Bloch energy loss formula
88 //
89 const Double_t kp1=0.76176e-1;
90 const Double_t kp2=10.632;
91 const Double_t kp3=0.13279e-4;
92 const Double_t kp4=1.8631;
93 const Double_t kp5=1.9479;
94
95 Double_t dbg = (Double_t) bg;
96
97 Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
98
99 Double_t aa = TMath::Power(beta,kp4);
100 Double_t bb = TMath::Power(1./dbg,kp5);
101
102 bb=TMath::Log(kp3+bb);
103
104 return ((Float_t)((kp2-aa-bb)*kp1/aa));
105}
106
51ad6848 107AliPointsMI::AliPointsMI(){
108 fN=0;
109 fX=0;
110 fY=0;
111 fZ=0;
112 fCapacity = 0;
113 fLabel0=0;
114 fLabel1=0;
115}
116
117
118AliPointsMI::AliPointsMI(Int_t n, Float_t *x,Float_t *y, Float_t *z){
119 //
120 //
121 fN=n;
122 fCapacity = 1000+2*n;
123 fX= new Float_t[n];
124 fY= new Float_t[n];
125 fZ= new Float_t[n];
126 memcpy(fX,x,n*sizeof(Float_t));
127 memcpy(fY,y,n*sizeof(Float_t));
128 memcpy(fZ,z,n*sizeof(Float_t));
129 fLabel0=0;
130 fLabel1=0;
131}
132
133void AliPointsMI::Reset()
134{
135 fN=0;
136}
137
138void AliPointsMI::Reset(AliDetector * det, Int_t particle){
139 //
140 // write points from detector points
141 //
142 Reset();
143 if (!det) return;
144 TObjArray *array = det->Points();
145 if (!array) return;
146 for (Int_t i=0;i<array->GetEntriesFast();i++){
147 AliPoints * points = (AliPoints*) array->At(i);
148 if (!points) continue;
149 if (points->GetIndex()!=particle) continue;
150 Int_t npoints = points->GetN();
151 if (npoints<2) continue;
152 Int_t delta = npoints/100;
153 if (delta<1) delta=1;
154 if (delta>10) delta=10;
155 Int_t mypoints = npoints/delta;
156 //
157 fN = mypoints;
158 if (fN>fCapacity){
159 fCapacity = 1000+2*fN;
160 delete []fX;
161 delete []fY;
162 delete []fZ;
163 fX = new Float_t[fCapacity];
164 fY = new Float_t[fCapacity];
165 fZ = new Float_t[fCapacity];
166 }
167 Float_t *xyz = points->GetP();
168 for (Int_t ipoint=0;ipoint<mypoints;ipoint++){
169 Int_t index = 3*ipoint*delta;
170 fX[ipoint]=0;
171 fY[ipoint]=0;
172 fZ[ipoint]=0;
173 if (index+2<npoints*3){
174 fX[ipoint] = xyz[index];
175 fY[ipoint] = xyz[index+1];
176 fZ[ipoint] = xyz[index+2];
177 }
178 }
179 }
180 fLabel0 = particle;
181}
182
183
184AliPointsMI::~AliPointsMI(){
185 fN=0;
186 fCapacity =0;
187 delete[] fX;
188 delete[]fY;
189 delete []fZ;
190 fX=0;fY=0;fZ=0;
191}
192
ae7d73d2 193
194
195
196
197////////////////////////////////////////////////////////////////////////
198AliMCInfo::AliMCInfo()
199{
200 fTPCReferences = new TClonesArray("AliTrackReference",10);
201 fITSReferences = new TClonesArray("AliTrackReference",10);
202 fTRDReferences = new TClonesArray("AliTrackReference",10);
51ad6848 203 fTOFReferences = new TClonesArray("AliTrackReference",10);
ae7d73d2 204 fTRdecay.SetTrack(-1);
205}
206
207AliMCInfo::~AliMCInfo()
208{
209 if (fTPCReferences) {
210 delete fTPCReferences;
211 }
212 if (fITSReferences){
213 delete fITSReferences;
214 }
215 if (fTRDReferences){
216 delete fTRDReferences;
217 }
51ad6848 218 if (fTOFReferences){
219 delete fTOFReferences;
220 }
221
ae7d73d2 222}
223
224
225
226void AliMCInfo::Update()
227{
228 //
229 //
230 fMCtracks =1;
231 if (!fTPCReferences) {
232 fNTPCRef =0;
233 return;
234 }
235 Float_t direction=1;
236 //Float_t rlast=0;
237 fNTPCRef = fTPCReferences->GetEntriesFast();
238 fNITSRef = fITSReferences->GetEntriesFast();
51ad6848 239 fNTRDRef = fTRDReferences->GetEntriesFast();
240 fNTOFRef = fTOFReferences->GetEntriesFast();
241
ae7d73d2 242 for (Int_t iref =0;iref<fTPCReferences->GetEntriesFast();iref++){
243 AliTrackReference * ref = (AliTrackReference *) fTPCReferences->At(iref);
244 //Float_t r = (ref->X()*ref->X()+ref->Y()*ref->Y());
245 Float_t newdirection = ref->X()*ref->Px()+ref->Y()*ref->Py(); //inside or outside
246 if (iref==0) direction = newdirection;
247 if ( newdirection*direction<0){
248 //changed direction
249 direction = newdirection;
250 fMCtracks+=1;
251 }
252 //rlast=r;
253 }
254 //
255 // decay info
256 fTPCdecay=kFALSE;
257 if (fTRdecay.GetTrack()>0){
258 fDecayCoord[0] = fTRdecay.X();
259 fDecayCoord[1] = fTRdecay.Y();
260 fDecayCoord[2] = fTRdecay.Z();
261 if ( (fTRdecay.R()<250)&&(fTRdecay.R()>85) && (TMath::Abs(fTRdecay.Z())<250) ){
262 fTPCdecay=kTRUE;
263 }
264 else{
265 fDecayCoord[0] = 0;
266 fDecayCoord[1] = 0;
267 fDecayCoord[2] = 0;
268 }
269 }
270 //
271 //
272 //digits information update
273 fRowsWithDigits = fTPCRow.RowsOn();
274 fRowsWithDigitsInn = fTPCRow.RowsOn(63); // 63 = number of inner rows
275 fRowsTrackLength = fTPCRow.Last() - fTPCRow.First();
276 //
277 //
278 // calculate primary ionization per cm
279 if (fParticle.GetPDG()){
280 fMass = fParticle.GetMass();
281 fCharge = fParticle.GetPDG()->Charge();
282 if (fTPCReferences->GetEntriesFast()>0){
283 fTrackRef = *((AliTrackReference*)fTPCReferences->At(0));
284 }
285 if (fMass>0){
286 Float_t p = TMath::Sqrt(fTrackRef.Px()*fTrackRef.Px()+
287 fTrackRef.Py()*fTrackRef.Py()+
288 fTrackRef.Pz()*fTrackRef.Pz());
289 if (p>0.001){
290 Float_t betagama = p /fMass;
291 fPrim = TPCBetheBloch(betagama);
292 }else fPrim=0;
293 }
294 }else{
295 fMass =0;
296 fPrim =0;
297 }
298}
299
300/////////////////////////////////////////////////////////////////////////////////
51ad6848 301/*
ae7d73d2 302void AliGenV0Info::Update()
303{
304 fMCPd[0] = fMCd.fParticle.Px();
305 fMCPd[1] = fMCd.fParticle.Py();
306 fMCPd[2] = fMCd.fParticle.Pz();
307 fMCPd[3] = fMCd.fParticle.P();
51ad6848 308 //
ae7d73d2 309 fMCX[0] = fMCd.fParticle.Vx();
310 fMCX[1] = fMCd.fParticle.Vy();
311 fMCX[2] = fMCd.fParticle.Vz();
312 fMCR = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
51ad6848 313 //
ae7d73d2 314 fPdg[0] = fMCd.fParticle.GetPdgCode();
315 fPdg[1] = fMCm.fParticle.GetPdgCode();
316 //
317 fLab[0] = fMCd.fParticle.GetUniqueID();
318 fLab[1] = fMCm.fParticle.GetUniqueID();
51ad6848 319 //
320}
321*/
ae7d73d2 322
51ad6848 323void AliGenV0Info::Update(Float_t vertex[3])
324{
325 fMCPd[0] = fMCd.fParticle.Px();
326 fMCPd[1] = fMCd.fParticle.Py();
327 fMCPd[2] = fMCd.fParticle.Pz();
328 fMCPd[3] = fMCd.fParticle.P();
329 //
330 fMCX[0] = fMCd.fParticle.Vx();
331 fMCX[1] = fMCd.fParticle.Vy();
332 fMCX[2] = fMCd.fParticle.Vz();
333 fMCR = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
334 //
335 fPdg[0] = fMCd.fParticle.GetPdgCode();
336 fPdg[1] = fMCm.fParticle.GetPdgCode();
337 //
338 fLab[0] = fMCd.fParticle.GetUniqueID();
339 fLab[1] = fMCm.fParticle.GetUniqueID();
340 //
341 //
342 //
343 Double_t x1[3],p1[3];
344 TParticle & pdaughter = fMCd.fParticle;
345 x1[0] = pdaughter.Vx();
346 x1[1] = pdaughter.Vy();
347 x1[2] = pdaughter.Vz();
348 p1[0] = pdaughter.Px();
349 p1[1] = pdaughter.Py();
350 p1[2] = pdaughter.Pz();
351 Double_t sign = (pdaughter.GetPDG()->Charge()>0)? -1:1;
352 AliHelix dhelix1(x1,p1,sign);
353 //
354 //
355 Double_t x2[3],p2[3];
356 //
357 TParticle & pmother = fMCm.fParticle;
358 x2[0] = pmother.Vx();
359 x2[1] = pmother.Vy();
360 x2[2] = pmother.Vz();
361 p2[0] = pmother.Px();
362 p2[1] = pmother.Py();
363 p2[2] = pmother.Pz();
364 //
365 //
366 sign = (pmother.GetPDG()->Charge()>0) ? -1:1;
367 AliHelix mhelix(x2,p2,sign);
368
369 //
370 //
371 //
372 //find intersection linear
373 //
374 Double_t distance1, distance2;
375 Double_t phase[2][2],radius[2];
376 Int_t points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
377 Double_t delta1=10000,delta2=10000;
378 if (points>0){
379 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
380 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
381 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
382 }
383 else{
384 fInvMass=-1;
385 return;
386 }
387 if (points==2){
388 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
389 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
390 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
391 }
392 distance1 = TMath::Min(delta1,delta2);
393 //
394 //find intersection parabolic
395 //
396 points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
397 delta1=10000,delta2=10000;
398
399 if (points>0){
400 dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
401 }
402 if (points==2){
403 dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
404 }
405
406 distance2 = TMath::Min(delta1,delta2);
407 //
408 if (delta1<delta2){
409 //get V0 info
410 dhelix1.Evaluate(phase[0][0],fMCXr);
411 dhelix1.GetMomentum(phase[0][0],fMCPdr);
412 mhelix.GetMomentum(phase[0][1],fMCPm);
413 dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fMCAngle);
414 fMCRr = TMath::Sqrt(radius[0]);
415 }
416 else{
417 dhelix1.Evaluate(phase[1][0],fMCXr);
418 dhelix1.GetMomentum(phase[1][0],fMCPdr);
419 mhelix.GetMomentum(phase[1][1],fMCPm);
420 dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fMCAngle);
421 fMCRr = TMath::Sqrt(radius[1]);
422 }
423 //
424 //
425 fMCDist1 = TMath::Sqrt(distance1);
426 fMCDist2 = TMath::Sqrt(distance2);
427 //
428 //
429 //
430 Float_t v[3] = {fMCXr[0]-vertex[0],fMCXr[1]-vertex[1],fMCXr[2]-vertex[2]};
431 //Float_t v[3] = {fMCXr[0],fMCXr[1],fMCXr[2]};
432 Float_t p[3] = {fMCPdr[0]+fMCPm[0], fMCPdr[1]+fMCPm[1],fMCPdr[2]+fMCPm[2]};
433 Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
434 Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
435 vnorm2 = TMath::Sqrt(vnorm2);
436 Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
437 Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
438 pnorm2 = TMath::Sqrt(pnorm2);
439 //
440 fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
441 fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);
442 fPointAngle = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
443 Float_t mass1 = fMCd.fMass;
444 Float_t mass2 = fMCm.fMass;
445
446
447 Float_t e1 = TMath::Sqrt(mass1*mass1+
448 fMCPdr[0]*fMCPdr[0]+
449 fMCPdr[1]*fMCPdr[1]+
450 fMCPdr[2]*fMCPdr[2]);
451 Float_t e2 = TMath::Sqrt(mass2*mass2+
452 fMCPm[0]*fMCPm[0]+
453 fMCPm[1]*fMCPm[1]+
454 fMCPm[2]*fMCPm[2]);
455
456 fInvMass =
457 (fMCPm[0]+fMCPdr[0])*(fMCPm[0]+fMCPdr[0])+
458 (fMCPm[1]+fMCPdr[1])*(fMCPm[1]+fMCPdr[1])+
459 (fMCPm[2]+fMCPdr[2])*(fMCPm[2]+fMCPdr[2]);
460
461 fInvMass = TMath::Sqrt((e1+e2)*(e1+e2)-fInvMass);
462
463
464}
465
466
467
468
469
470/////////////////////////////////////////////////////////////////////////////////
471void AliGenKinkInfo::Update()
472{
473 fMCPd[0] = fMCd.fParticle.Px();
474 fMCPd[1] = fMCd.fParticle.Py();
475 fMCPd[2] = fMCd.fParticle.Pz();
476 fMCPd[3] = fMCd.fParticle.P();
477 //
478 fMCX[0] = fMCd.fParticle.Vx();
479 fMCX[1] = fMCd.fParticle.Vy();
480 fMCX[2] = fMCd.fParticle.Vz();
481 fMCR = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
482 //
483 fPdg[0] = fMCd.fParticle.GetPdgCode();
484 fPdg[1] = fMCm.fParticle.GetPdgCode();
485 //
486 fLab[0] = fMCd.fParticle.GetUniqueID();
487 fLab[1] = fMCm.fParticle.GetUniqueID();
488 //
489 //
490 //
491 Double_t x1[3],p1[3];
492 TParticle & pdaughter = fMCd.fParticle;
493 x1[0] = pdaughter.Vx();
494 x1[1] = pdaughter.Vy();
495 x1[2] = pdaughter.Vz();
496 p1[0] = pdaughter.Px();
497 p1[1] = pdaughter.Py();
498 p1[2] = pdaughter.Pz();
499 Double_t sign = (pdaughter.GetPDG()->Charge()>0)? -1:1;
500 AliHelix dhelix1(x1,p1,sign);
501 //
502 //
503 Double_t x2[3],p2[3];
504 //
505 TParticle & pmother = fMCm.fParticle;
506 x2[0] = pmother.Vx();
507 x2[1] = pmother.Vy();
508 x2[2] = pmother.Vz();
509 p2[0] = pmother.Px();
510 p2[1] = pmother.Py();
511 p2[2] = pmother.Pz();
512 //
513 AliTrackReference & pdecay = fMCm.fTRdecay;
514 x2[0] = pdecay.X();
515 x2[1] = pdecay.Y();
516 x2[2] = pdecay.Z();
517 p2[0] = pdecay.Px();
518 p2[1] = pdecay.Py();
519 p2[2] = pdecay.Pz();
520 //
521 sign = (pmother.GetPDG()->Charge()>0) ? -1:1;
522 AliHelix mhelix(x2,p2,sign);
523
524 //
525 //
526 //
527 //find intersection linear
528 //
529 Double_t distance1, distance2;
530 Double_t phase[2][2],radius[2];
531 Int_t points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
532 Double_t delta1=10000,delta2=10000;
533 if (points>0){
534 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
535 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
536 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
537 }
538 if (points==2){
539 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
540 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
541 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
542 }
543 distance1 = TMath::Min(delta1,delta2);
544 //
545 //find intersection parabolic
546 //
547 points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
548 delta1=10000,delta2=10000;
549
550 if (points>0){
551 dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
552 }
553 if (points==2){
554 dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
555 }
556
557 distance2 = TMath::Min(delta1,delta2);
558 //
559 if (delta1<delta2){
560 //get V0 info
561 dhelix1.Evaluate(phase[0][0],fMCXr);
562 dhelix1.GetMomentum(phase[0][0],fMCPdr);
563 mhelix.GetMomentum(phase[0][1],fMCPm);
564 dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fMCAngle);
565 fMCRr = TMath::Sqrt(radius[0]);
566 }
567 else{
568 dhelix1.Evaluate(phase[1][0],fMCXr);
569 dhelix1.GetMomentum(phase[1][0],fMCPdr);
570 mhelix.GetMomentum(phase[1][1],fMCPm);
571 dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fMCAngle);
572 fMCRr = TMath::Sqrt(radius[1]);
573 }
574 //
575 //
576 fMCDist1 = TMath::Sqrt(distance1);
577 fMCDist2 = TMath::Sqrt(distance2);
578
ae7d73d2 579}
580
581
51ad6848 582Float_t AliGenKinkInfo::GetQt(){
583 //
584 //
585 Float_t momentumd = TMath::Sqrt(fMCPd[0]*fMCPd[0]+fMCPd[1]*fMCPd[1]+fMCPd[2]*fMCPd[2]);
586 return TMath::Sin(fMCAngle[2])*momentumd;
587}
588
ae7d73d2 589
590
591
592////////////////////////////////////////////////////////////////////////
593
594////////////////////////////////////////////////////////////////////////
595//
596// End of implementation of the class AliMCInfo
597//
598////////////////////////////////////////////////////////////////////////
599
600
601
602////////////////////////////////////////////////////////////////////////
603digitRow::digitRow()
604{
605 Reset();
606}
607////////////////////////////////////////////////////////////////////////
608digitRow & digitRow::operator=(const digitRow &digOld)
609{
610 for (Int_t i = 0; i<kgRowBytes; i++) fDig[i] = digOld.fDig[i];
611 return (*this);
612}
613////////////////////////////////////////////////////////////////////////
614void digitRow::SetRow(Int_t row)
615{
616 if (row >= 8*kgRowBytes) {
617 cerr<<"digitRow::SetRow: index "<<row<<" out of bounds."<<endl;
618 return;
619 }
620 Int_t iC = row/8;
621 Int_t iB = row%8;
622 SETBIT(fDig[iC],iB);
623}
624
625////////////////////////////////////////////////////////////////////////
626Bool_t digitRow::TestRow(Int_t row)
627{
628//
629// return kTRUE if row is on
630//
631 Int_t iC = row/8;
632 Int_t iB = row%8;
633 return TESTBIT(fDig[iC],iB);
634}
635////////////////////////////////////////////////////////////////////////
636Int_t digitRow::RowsOn(Int_t upto)
637{
638//
639// returns number of rows with a digit
640// count only rows less equal row number upto
641//
642 Int_t total = 0;
643 for (Int_t i = 0; i<kgRowBytes; i++) {
644 for (Int_t j = 0; j < 8; j++) {
645 if (i*8+j > upto) return total;
646 if (TESTBIT(fDig[i],j)) total++;
647 }
648 }
649 return total;
650}
651////////////////////////////////////////////////////////////////////////
652void digitRow::Reset()
653{
654//
655// resets all rows to zero
656//
657 for (Int_t i = 0; i<kgRowBytes; i++) {
658 fDig[i] <<= 8;
659 }
660}
661////////////////////////////////////////////////////////////////////////
662Int_t digitRow::Last()
663{
664//
665// returns the last row number with a digit
666// returns -1 if now digits
667//
668 for (Int_t i = kgRowBytes-1; i>=0; i--) {
669 for (Int_t j = 7; j >= 0; j--) {
670 if TESTBIT(fDig[i],j) return i*8+j;
671 }
672 }
673 return -1;
674}
675////////////////////////////////////////////////////////////////////////
676Int_t digitRow::First()
677{
678//
679// returns the first row number with a digit
680// returns -1 if now digits
681//
682 for (Int_t i = 0; i<kgRowBytes; i++) {
683 for (Int_t j = 0; j < 8; j++) {
684 if (TESTBIT(fDig[i],j)) return i*8+j;
685 }
686 }
687 return -1;
688}
689
690////////////////////////////////////////////////////////////////////////
691//
692// end of implementation of a class digitRow
693//
694////////////////////////////////////////////////////////////////////////
695
696////////////////////////////////////////////////////////////////////////
697AliGenInfoMaker::AliGenInfoMaker()
698{
699 Reset();
700}
701
702////////////////////////////////////////////////////////////////////////
703AliGenInfoMaker::AliGenInfoMaker(const char * fnGalice, const char* fnRes,
704 Int_t nEvents, Int_t firstEvent)
705{
706 Reset();
707 fFirstEventNr = firstEvent;
708 fEventNr = firstEvent;
709 fNEvents = nEvents;
710 // fFnRes = fnRes;
711 sprintf(fFnRes,"%s",fnRes);
712 //
713 fLoader = AliRunLoader::Open(fnGalice);
714 if (gAlice){
715 delete gAlice->GetRunLoader();
716 delete gAlice;
717 gAlice = 0x0;
718 }
719 if (fLoader->LoadgAlice()){
720 cerr<<"Error occured while l"<<endl;
721 }
722 Int_t nall = fLoader->GetNumberOfEvents();
723 if (nEvents==0) {
724 nEvents =nall;
725 fNEvents=nall;
726 fFirstEventNr=0;
727 }
728
729 if (nall<=0){
730 cerr<<"no events available"<<endl;
731 fEventNr = 0;
732 return;
733 }
734 if (firstEvent+nEvents>nall) {
735 fEventNr = nall-firstEvent;
736 cerr<<"restricted number of events availaible"<<endl;
737 }
738 AliMagF * magf = gAlice->Field();
739 AliTracker::SetFieldMap(magf);
740}
741
742
743AliMCInfo * AliGenInfoMaker::MakeInfo(UInt_t i)
744{
745 //
746 if (i<fNParticles) {
747 if (fGenInfo[i]) return fGenInfo[i];
748 fGenInfo[i] = new AliMCInfo;
749 fNInfos++;
750 return fGenInfo[i];
751 }
752 else
753 return 0;
754}
755
756////////////////////////////////////////////////////////////////////////
757void AliGenInfoMaker::Reset()
758{
759 fEventNr = 0;
760 fNEvents = 0;
761 fTreeGenTracks = 0;
762 fFileGenTracks = 0;
763 fGenInfo = 0;
764 fNInfos = 0;
765 //
766 //
767 fDebug = 0;
768 fVPrim[0] = -1000.;
769 fVPrim[1] = -1000.;
770 fVPrim[2] = -1000.;
771 fParamTPC = 0;
772}
773////////////////////////////////////////////////////////////////////////
774AliGenInfoMaker::~AliGenInfoMaker()
775{
776
777 if (fLoader){
778 fLoader->UnloadgAlice();
779 gAlice = 0;
780 delete fLoader;
781 }
782}
783
784Int_t AliGenInfoMaker::SetIO()
785{
786 //
787 //
788 CreateTreeGenTracks();
789 if (!fTreeGenTracks) return 1;
790 // AliTracker::SetFieldFactor();
791
792 fParamTPC = GetTPCParam();
793 //
794 return 0;
795}
796
797////////////////////////////////////////////////////////////////////////
798Int_t AliGenInfoMaker::SetIO(Int_t eventNr)
799{
800 //
801 //
802 // SET INPUT
803 fLoader->SetEventNumber(eventNr);
804 //
805 fLoader->LoadHeader();
806 fLoader->LoadKinematics();
807 fStack = fLoader->Stack();
808 //
809 fLoader->LoadTrackRefs();
51ad6848 810 fLoader->LoadHits();
ae7d73d2 811 fTreeTR = fLoader->TreeTR();
812 //
813 AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
814 tpcl->LoadDigits();
815 fTreeD = tpcl->TreeD();
816 return 0;
817}
818
819Int_t AliGenInfoMaker::CloseIOEvent()
820{
821 fLoader->UnloadHeader();
822 fLoader->UnloadKinematics();
823 fLoader->UnloadTrackRefs();
824 AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
825 tpcl->UnloadDigits();
826 return 0;
827}
828
829Int_t AliGenInfoMaker::CloseIO()
830{
831 fLoader->UnloadgAlice();
832 return 0;
833}
834
835
836
837////////////////////////////////////////////////////////////////////////
838Int_t AliGenInfoMaker::Exec(Int_t nEvents, Int_t firstEventNr)
839{
840 fNEvents = nEvents;
841 fFirstEventNr = firstEventNr;
842 return Exec();
843}
844
845////////////////////////////////////////////////////////////////////////
846Int_t AliGenInfoMaker::Exec()
847{
848 TStopwatch timer;
849 timer.Start();
850 Int_t status =SetIO();
851 if (status>0) return status;
852 //
853
854 for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
855 fEventNr++) {
856 SetIO(fEventNr);
857 fNParticles = fStack->GetNtrack();
858 //
859 fGenInfo = new AliMCInfo*[fNParticles];
860 for (UInt_t i = 0; i<fNParticles; i++) {
861 fGenInfo[i]=0;
862 }
863 //
864 cout<<"Start to process event "<<fEventNr<<endl;
865 cout<<"\tfNParticles = "<<fNParticles<<endl;
866 if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeTR"<<endl;
867 if (TreeTRLoop()>0) return 1;
868 //
869 if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeD"<<endl;
870 if (TreeDLoop()>0) return 1;
871 //
872 if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeK"<<endl;
873 if (TreeKLoop()>0) return 1;
874 if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
51ad6848 875 //
876 if (BuildKinkInfo()>0) return 1;
877 if (BuildV0Info()>0) return 1;
878 //if (BuildHitLines()>0) return 1;
879 if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
880 //
ae7d73d2 881 for (UInt_t i = 0; i<fNParticles; i++) {
882 if (fGenInfo[i]) delete fGenInfo[i];
883 }
884 delete []fGenInfo;
885 CloseIOEvent();
886 }
887 //
888 CloseIO();
889 CloseOutputFile();
890
891 cerr<<"Exec finished"<<endl;
892
893 timer.Stop();
894 timer.Print();
895 return 0;
896}
897////////////////////////////////////////////////////////////////////////
898void AliGenInfoMaker::CreateTreeGenTracks()
899{
900 fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
901 if (!fFileGenTracks) {
902 cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
903 return;
904 }
905 fTreeGenTracks = new TTree("genTracksTree","genTracksTree");
906 AliMCInfo * info = new AliMCInfo;
ae7d73d2 907 fTreeGenTracks->Branch("MC","AliMCInfo",&info);
908 delete info;
51ad6848 909 //
910 AliGenKinkInfo *kinkinfo = new AliGenKinkInfo;
911 fTreeKinks = new TTree("genKinksTree","genKinksTree");
912 fTreeKinks->Branch("MC","AliGenKinkInfo",&kinkinfo);
913 delete kinkinfo;
914 //
915 AliGenV0Info *v0info = new AliGenV0Info;
916 fTreeV0 = new TTree("genV0Tree","genV0Tree");
917 fTreeV0->Branch("MC","AliGenV0Info",&v0info);
918 delete v0info;
919 //
920 //
921 AliPointsMI * points0 = new AliPointsMI;
922 AliPointsMI * points1 = new AliPointsMI;
923 AliPointsMI * points2 = new AliPointsMI;
924 fTreeHitLines = new TTree("HitLines","HitLines");
925 fTreeHitLines->Branch("TPC.","AliPointsMI",&points0,32000,10);
926 fTreeHitLines->Branch("TRD.","AliPointsMI",&points1,32000,10);
927 fTreeHitLines->Branch("ITS.","AliPointsMI",&points2,32000,10);
928 Int_t label=0;
929 fTreeHitLines->Branch("Label",&label,"label/I");
930 //
ae7d73d2 931 fTreeGenTracks->AutoSave();
51ad6848 932 fTreeKinks->AutoSave();
933 fTreeV0->AutoSave();
934 fTreeHitLines->AutoSave();
ae7d73d2 935}
936////////////////////////////////////////////////////////////////////////
937void AliGenInfoMaker::CloseOutputFile()
938{
939 if (!fFileGenTracks) {
940 cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
941 return;
942 }
943 fFileGenTracks->cd();
944 fTreeGenTracks->Write();
945 delete fTreeGenTracks;
51ad6848 946 fTreeKinks->Write();
947 delete fTreeKinks;
948 fTreeV0->Write();
949 delete fTreeV0;
950
ae7d73d2 951 fFileGenTracks->Close();
952 delete fFileGenTracks;
953 return;
954}
955
956////////////////////////////////////////////////////////////////////////
957Int_t AliGenInfoMaker::TreeKLoop()
958{
959//
960// open the file with treeK
961// loop over all entries there and save information about some tracks
962//
963
964 AliStack * stack = fStack;
965 if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
966
967 if (fDebug > 0) {
968 cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
969 <<fEventNr<<endl;
970 }
971 Int_t ipdg = 0;
972 TParticlePDG *ppdg = 0;
973 // not all generators give primary vertex position. Take the vertex
974 // of the particle 0 as primary vertex.
975 TDatabasePDG pdg; //get pdg table
976 //thank you very much root for this
977 TBranch * br = fTreeGenTracks->GetBranch("MC");
978 TParticle *particle = stack->ParticleFromTreeK(0);
979 fVPrim[0] = particle->Vx();
980 fVPrim[1] = particle->Vy();
981 fVPrim[2] = particle->Vz();
982 for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
983 // load only particles with TR
984 AliMCInfo * info = GetInfo(iParticle);
985 if (!info) continue;
986 //////////////////////////////////////////////////////////////////////
987 info->fLabel = iParticle;
988 //
989 info->fParticle = *(stack->Particle(iParticle));
990 info->fVDist[0] = info->fParticle.Vx()-fVPrim[0];
991 info->fVDist[1] = info->fParticle.Vy()-fVPrim[1];
992 info->fVDist[2] = info->fParticle.Vz()-fVPrim[2];
993 info->fVDist[3] = TMath::Sqrt(info->fVDist[0]*info->fVDist[0]+
994 info->fVDist[1]*info->fVDist[1]+info->fVDist[2]*info->fVDist[2]);
995 //
996 //
997 ipdg = info->fParticle.GetPdgCode();
998 info->fPdg = ipdg;
999 ppdg = pdg.GetParticle(ipdg);
1000 info->fEventNr = fEventNr;
1001 info->Update();
1002 //////////////////////////////////////////////////////////////////////
1003 br->SetAddress(&info);
51ad6848 1004 fTreeGenTracks->Fill();
ae7d73d2 1005 }
1006 fTreeGenTracks->AutoSave();
1007 if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
1008 return 0;
1009}
1010
1011
51ad6848 1012////////////////////////////////////////////////////////////////////////////////////
1013Int_t AliGenInfoMaker::BuildKinkInfo()
1014{
1015 //
1016 // Build tree with Kink Information
1017 //
1018 AliStack * stack = fStack;
1019 if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
1020
1021 if (fDebug > 0) {
1022 cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
1023 <<fEventNr<<endl;
1024 }
1025 // Int_t ipdg = 0;
1026 //TParticlePDG *ppdg = 0;
1027 // not all generators give primary vertex position. Take the vertex
1028 // of the particle 0 as primary vertex.
1029 TDatabasePDG pdg; //get pdg table
1030 //thank you very much root for this
1031 TBranch * br = fTreeKinks->GetBranch("MC");
1032 //
1033 AliGenKinkInfo * kinkinfo = new AliGenKinkInfo;
1034 //
1035 for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
1036 // load only particles with TR
1037 AliMCInfo * info = GetInfo(iParticle);
1038 if (!info) continue;
1039 if (info->fCharge==0) continue;
1040 if (info->fTRdecay.P()<0.13) continue; //momenta cut
1041 if (info->fTRdecay.R()>500) continue; //R cut - decay outside barrel
1042 TParticle & particle = info->fParticle;
1043 Int_t first = particle.GetDaughter(0);
1044 if (first ==0) continue;
1045
1046 Int_t last = particle.GetDaughter(1);
1047 if (last ==0) last=first;
1048 AliMCInfo * info2 = 0;
1049 AliMCInfo * dinfo = 0;
1050
1051 for (Int_t id2=first;id2<=last;id2++){
1052 info2 = GetInfo(id2);
1053 if (!info2) continue;
1054 TParticle &p2 = info2->fParticle;
1055 Double_t r = TMath::Sqrt(p2.Vx()*p2.Vx()+p2.Vy()*p2.Vy());
1056 if ( TMath::Abs(info->fTRdecay.R()-r)>0.01) continue;
1057 if (!(p2.GetPDG())) continue;
1058 dinfo =info2;
1059
1060 kinkinfo->fMCm = (*info);
1061 kinkinfo->fMCd = (*dinfo);
1062 if (kinkinfo->fMCm.fParticle.GetPDG()==0 || kinkinfo->fMCd.fParticle.GetPDG()==0) continue;
1063 kinkinfo->Update();
1064 br->SetAddress(&kinkinfo);
1065 fTreeKinks->Fill();
1066 }
1067 /*
1068 if (dinfo){
1069 kinkinfo->fMCm = (*info);
1070 kinkinfo->fMCd = (*dinfo);
1071 kinkinfo->Update();
1072 br->SetAddress(&kinkinfo);
1073 fTreeKinks->Fill();
1074 }
1075 */
1076 }
1077 fTreeGenTracks->AutoSave();
1078 if (fDebug > 2) cerr<<"end of Kink Loop"<<endl;
1079 return 0;
1080}
1081
1082
1083////////////////////////////////////////////////////////////////////////////////////
1084Int_t AliGenInfoMaker::BuildV0Info()
1085{
1086 //
1087 // Build tree with V0 Information
1088 //
1089 AliStack * stack = fStack;
1090 if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
1091
1092 if (fDebug > 0) {
1093 cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
1094 <<fEventNr<<endl;
1095 }
1096 // Int_t ipdg = 0;
1097 //TParticlePDG *ppdg = 0;
1098 // not all generators give primary vertex position. Take the vertex
1099 // of the particle 0 as primary vertex.
1100 TDatabasePDG pdg; //get pdg table
1101 //thank you very much root for this
1102 TBranch * br = fTreeV0->GetBranch("MC");
1103 //
1104 AliGenV0Info * v0info = new AliGenV0Info;
1105 //
1106 //
1107 for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
1108 // load only particles with TR
1109 AliMCInfo * info = GetInfo(iParticle);
1110 if (!info) continue;
1111 if (info->fCharge==0) continue;
1112 //
1113 //
1114 TParticle & particle = info->fParticle;
1115 Int_t mother = particle.GetMother(0);
1116 if (mother <=0) continue;
1117 //
1118 TParticle * motherparticle = stack->Particle(mother);
1119 if (!motherparticle) continue;
1120 //
1121 Int_t last = motherparticle->GetDaughter(1);
1122 if (last==(int)iParticle) continue;
1123 AliMCInfo * info2 = info;
1124 AliMCInfo * dinfo = GetInfo(last);
1125 if (!dinfo) continue;
1126 if (!dinfo->fParticle.GetPDG()) continue;
1127 if (!info2->fParticle.GetPDG()) continue;
1128 if (dinfo){
1129 v0info->fMCm = (*info);
1130 v0info->fMCd = (*dinfo);
1131 v0info->Update(fVPrim);
1132 br->SetAddress(&v0info);
1133 fTreeV0->Fill();
1134 }
1135 }
1136 fTreeV0->AutoSave();
1137 if (fDebug > 2) cerr<<"end of V0 Loop"<<endl;
1138 return 0;
1139}
1140
1141
1142
1143
1144////////////////////////////////////////////////////////////////////////
1145Int_t AliGenInfoMaker::BuildHitLines()
1146{
1147
1148//
1149// open the file with treeK
1150// loop over all entries there and save information about some tracks
1151//
1152
1153 AliStack * stack = fStack;
1154 if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
1155
1156 if (fDebug > 0) {
1157 cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
1158 <<fEventNr<<endl;
1159 }
1160 Int_t ipdg = 0;
1161 // TParticlePDG *ppdg = 0;
1162 // not all generators give primary vertex position. Take the vertex
1163 // of the particle 0 as primary vertex.
1164 TDatabasePDG pdg; //get pdg table
1165 //thank you very much root for this
1166 AliPointsMI *tpcp = new AliPointsMI;
1167 AliPointsMI *trdp = new AliPointsMI;
1168 AliPointsMI *itsp = new AliPointsMI;
1169 Int_t label =0;
1170 //
1171 TBranch * brtpc = fTreeHitLines->GetBranch("TPC.");
1172 TBranch * brtrd = fTreeHitLines->GetBranch("TRD.");
1173 TBranch * brits = fTreeHitLines->GetBranch("ITS.");
1174 TBranch * brlabel = fTreeHitLines->GetBranch("Label");
1175 brlabel->SetAddress(&label);
1176 brtpc->SetAddress(&tpcp);
1177 brtrd->SetAddress(&trdp);
1178 brits->SetAddress(&itsp);
1179 //
1180 AliDetector *dtpc = gAlice->GetDetector("TPC");
1181 AliDetector *dtrd = gAlice->GetDetector("TRD");
1182 AliDetector *dits = gAlice->GetDetector("ITS");
1183
1184 for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
1185 // load only particles with TR
1186 AliMCInfo * info = GetInfo(iParticle);
1187 if (!info) continue;
1188 Int_t primpart = info->fPrimPart;
1189 ipdg = info->fParticle.GetPdgCode();
1190 label = iParticle;
1191 //
1192 gAlice->ResetHits();
1193 tpcp->Reset();
1194 itsp->Reset();
1195 trdp->Reset();
1196 tpcp->fLabel1 = ipdg;
1197 trdp->fLabel1 = ipdg;
1198 itsp->fLabel1 = ipdg;
1199 if (dtpc->TreeH()->GetEvent(primpart)){
1200 dtpc->LoadPoints(primpart);
1201 tpcp->Reset(dtpc,iParticle);
1202 }
1203 if (dtrd->TreeH()->GetEvent(primpart)){
1204 dtrd->LoadPoints(primpart);
1205 trdp->Reset(dtrd,iParticle);
1206 }
1207 if (dits->TreeH()->GetEvent(primpart)){
1208 dits->LoadPoints(primpart);
1209 itsp->Reset(dits,iParticle);
1210 }
1211 //
1212 fTreeHitLines->Fill();
1213 dtpc->ResetPoints();
1214 dtrd->ResetPoints();
1215 dits->ResetPoints();
1216 }
1217 delete tpcp;
1218 delete trdp;
1219 delete itsp;
1220 fTreeHitLines->AutoSave();
1221 if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
1222 return 0;
1223}
ae7d73d2 1224
1225
1226////////////////////////////////////////////////////////////////////////
1227Int_t AliGenInfoMaker::TreeDLoop()
1228{
1229 //
1230 // open the file with treeD
1231 // loop over all entries there and save information about some tracks
1232 //
1233
1234 Int_t nInnerSector = fParamTPC->GetNInnerSector();
1235 Int_t rowShift = 0;
1236 Int_t zero=fParamTPC->GetZeroSup()+6;
1237 //
1238 //
1239 AliSimDigits digitsAddress, *digits=&digitsAddress;
1240 fTreeD->GetBranch("Segment")->SetAddress(&digits);
1241
1242 Int_t sectorsByRows=(Int_t)fTreeD->GetEntries();
1243 if (fDebug > 1) cout<<"\tsectorsByRows = "<<sectorsByRows<<endl;
1244 for (Int_t i=0; i<sectorsByRows; i++) {
1245 if (!fTreeD->GetEvent(i)) continue;
1246 Int_t sec,row;
1247 fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
1248 if (fDebug > 1) cout<<sec<<' '<<row<<" \r";
1249 // here I expect that upper sectors follow lower sectors
1250 if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow();
1251 //
1252 digits->ExpandTrackBuffer();
1253 digits->First();
1254 do {
1255 Int_t iRow=digits->CurrentRow();
1256 Int_t iColumn=digits->CurrentColumn();
1257 Short_t digitValue = digits->CurrentDigit();
1258 if (digitValue >= zero) {
1259 Int_t label;
1260 for (Int_t j = 0; j<3; j++) {
1261 // label = digits->GetTrackID(iRow,iColumn,j);
1262 label = digits->GetTrackIDFast(iRow,iColumn,j)-2;
1263 if (label >= (Int_t)fNParticles) { //don't label from bakground event
1264 continue;
1265 }
1266 if (label >= 0 && label <= (Int_t)fNParticles) {
1267 if (fDebug > 6 ) {
1268 cout<<"Inner loop: sector, iRow, iColumn, label, value, row "
1269 <<sec<<" "
1270 <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
1271 <<" "<<row<<endl;
1272 }
1273 AliMCInfo * info = GetInfo(label);
1274 if (info){
1275 info->fTPCRow.SetRow(row+rowShift);
1276 }
1277 }
1278 }
1279 }
1280 } while (digits->Next());
1281 }
1282
1283 if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;
1284 return 0;
1285}
1286
1287
1288////////////////////////////////////////////////////////////////////////
1289Int_t AliGenInfoMaker::TreeTRLoop()
1290{
1291 //
1292 // loop over TrackReferences and store the first one for each track
1293 //
1294 TTree * treeTR = fTreeTR;
1295 Int_t nPrimaries = (Int_t) treeTR->GetEntries();
1296 if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl;
1297 //
1298 //
1299 //track references for TPC
1300 TClonesArray* TPCArrayTR = new TClonesArray("AliTrackReference");
1301 TClonesArray* ITSArrayTR = new TClonesArray("AliTrackReference");
1302 TClonesArray* TRDArrayTR = new TClonesArray("AliTrackReference");
51ad6848 1303 TClonesArray* TOFArrayTR = new TClonesArray("AliTrackReference");
ae7d73d2 1304 TClonesArray* RunArrayTR = new TClonesArray("AliTrackReference");
1305 //
1306 if (treeTR->GetBranch("TPC")) treeTR->GetBranch("TPC")->SetAddress(&TPCArrayTR);
1307 if (treeTR->GetBranch("ITS")) treeTR->GetBranch("ITS")->SetAddress(&ITSArrayTR);
1308 if (treeTR->GetBranch("TRD")) treeTR->GetBranch("TRD")->SetAddress(&TRDArrayTR);
51ad6848 1309 if (treeTR->GetBranch("TOF")) treeTR->GetBranch("TOF")->SetAddress(&TOFArrayTR);
ae7d73d2 1310 if (treeTR->GetBranch("AliRun")) treeTR->GetBranch("AliRun")->SetAddress(&RunArrayTR);
1311 //
1312 //
1313 //
1314 for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
1315 treeTR->GetEntry(iPrimPart);
1316 //
1317 // Loop over TPC references
1318 //
1319 for (Int_t iTrackRef = 0; iTrackRef < TPCArrayTR->GetEntriesFast(); iTrackRef++) {
1320 AliTrackReference *trackRef = (AliTrackReference*)TPCArrayTR->At(iTrackRef);
1321 //
51ad6848 1322 if (trackRef->P()<fgTPCPtCut) continue;
ae7d73d2 1323 Int_t label = trackRef->GetTrack();
1324 AliMCInfo * info = GetInfo(label);
1325 if (!info) info = MakeInfo(label);
1326 if (!info) continue;
51ad6848 1327 info->fPrimPart = iPrimPart;
ae7d73d2 1328 TClonesArray & arr = *(info->fTPCReferences);
1329 new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
1330 }
1331 //
1332 // Loop over ITS references
1333 //
1334 for (Int_t iTrackRef = 0; iTrackRef < ITSArrayTR->GetEntriesFast(); iTrackRef++) {
1335 AliTrackReference *trackRef = (AliTrackReference*)ITSArrayTR->At(iTrackRef);
1336 //
51ad6848 1337 if (trackRef->P()<fgTPCPtCut) continue;
ae7d73d2 1338 Int_t label = trackRef->GetTrack();
1339 AliMCInfo * info = GetInfo(label);
1340 if ( (!info) && trackRef->Pt()<fgITSPtCut) continue;
1341 if (!info) info = MakeInfo(label);
1342 if (!info) continue;
51ad6848 1343 info->fPrimPart = iPrimPart;
ae7d73d2 1344 TClonesArray & arr = *(info->fITSReferences);
1345 new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
1346 }
1347 //
1348 // Loop over TRD references
1349 //
1350 for (Int_t iTrackRef = 0; iTrackRef < TRDArrayTR->GetEntriesFast(); iTrackRef++) {
1351 AliTrackReference *trackRef = (AliTrackReference*)TRDArrayTR->At(iTrackRef);
1352 //
51ad6848 1353 if (trackRef->P()<fgTPCPtCut) continue;
ae7d73d2 1354 Int_t label = trackRef->GetTrack();
1355 AliMCInfo * info = GetInfo(label);
1356 if ( (!info) && trackRef->Pt()<fgTRDPtCut) continue;
1357 if (!info) info = MakeInfo(label);
1358 if (!info) continue;
51ad6848 1359 info->fPrimPart = iPrimPart;
ae7d73d2 1360 TClonesArray & arr = *(info->fTRDReferences);
1361 new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
1362 }
1363 //
51ad6848 1364 // Loop over TOF references
1365 //
1366 for (Int_t iTrackRef = 0; iTrackRef < TOFArrayTR->GetEntriesFast(); iTrackRef++) {
1367 AliTrackReference *trackRef = (AliTrackReference*)TOFArrayTR->At(iTrackRef);
1368 Int_t label = trackRef->GetTrack();
1369 AliMCInfo * info = GetInfo(label);
1370 if (!info){
1371 if (trackRef->Pt()<fgTPCPtCut) continue;
1372 if ( (!info) && trackRef->Pt()<fgTOFPtCut) continue;
1373 }
1374 if (!info) info = MakeInfo(label);
1375 if (!info) continue;
1376 info->fPrimPart = iPrimPart;
1377 TClonesArray & arr = *(info->fTOFReferences);
1378 new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
1379 }
1380 //
ae7d73d2 1381 // get dacay position
1382 for (Int_t iTrackRef = 0; iTrackRef < RunArrayTR->GetEntriesFast(); iTrackRef++) {
1383 AliTrackReference *trackRef = (AliTrackReference*)RunArrayTR->At(iTrackRef);
1384 //
1385 Int_t label = trackRef->GetTrack();
1386 AliMCInfo * info = GetInfo(label);
1387 if (!info) continue;
1388 if (!trackRef->TestBit(BIT(2))) continue; //if not decay
1389 // if (TMath::Abs(trackRef.X());
1390 info->fTRdecay = *trackRef;
1391 }
1392 }
1393 //
1394 TPCArrayTR->Delete();
1395 delete TPCArrayTR;
1396 TRDArrayTR->Delete();
1397 delete TRDArrayTR;
51ad6848 1398 TOFArrayTR->Delete();
1399 delete TOFArrayTR;
1400
ae7d73d2 1401 ITSArrayTR->Delete();
1402 delete ITSArrayTR;
1403 RunArrayTR->Delete();
1404 delete RunArrayTR;
1405 //
1406 return 0;
1407}
1408
1409////////////////////////////////////////////////////////////////////////
1410Float_t AliGenInfoMaker::TR2LocalX(AliTrackReference *trackRef,
1411 AliTPCParam *paramTPC) {
1412
1413 Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
1414 Int_t index[4];
1415 paramTPC->Transform0to1(x,index);
1416 paramTPC->Transform1to2(x,index);
1417 return x[0];
1418}
1419////////////////////////////////////////////////////////////////////////
1420
1421
1422
1423TH1F * AliComparisonDraw::DrawXY(const char * chx, const char *chy, const char* selection,
1424 const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy)
1425{
1426 //
1427 Double_t* bins = CreateLogBins(nbins, minx, maxx);
1428 Int_t nBinsRes = 30;
1429 TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, minx, maxx, nBinsRes, miny, maxy);
1430 char cut[1000];
1431 sprintf(cut,"%s&&%s",selection,quality);
1432 char expression[1000];
1433 sprintf(expression,"%s:%s>>hRes2",chy,chx);
1434 fTree->Draw(expression, cut, "groff");
1435 TH1F* hMean=0;
1436 TH1F* hRes = CreateResHisto(hRes2, &hMean);
1437 AliLabelAxes(hRes, chx, chy);
1438 //
1439 delete hRes2;
1440 delete[] bins;
1441 fRes = hRes;
1442 fMean = hMean;
1443 return hRes;
1444}
1445
1446TH1F * AliComparisonDraw::DrawLogXY(const char * chx, const char *chy, const char* selection,
1447 const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy)
1448{
1449 //
1450 Double_t* bins = CreateLogBins(nbins, minx, maxx);
1451 Int_t nBinsRes = 30;
1452 TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, bins, nBinsRes, miny, maxy);
1453 char cut[1000];
1454 sprintf(cut,"%s&&%s",selection,quality);
1455 char expression[1000];
1456 sprintf(expression,"%s:%s>>hRes2",chy,chx);
1457 fTree->Draw(expression, cut, "groff");
1458 TH1F* hMean=0;
1459 TH1F* hRes = CreateResHisto(hRes2, &hMean);
1460 AliLabelAxes(hRes, chx, chy);
1461 //
1462 delete hRes2;
1463 delete[] bins;
1464 fRes = hRes;
1465 fMean = hMean;
1466 return hRes;
1467}
1468
1469///////////////////////////////////////////////////////////////////////////////////
1470///////////////////////////////////////////////////////////////////////////////////
1471TH1F * AliComparisonDraw::Eff(const char *variable, const char* selection, const char * quality,
1472 Int_t nbins, Float_t min, Float_t max)
1473{
1474 //
1475 //
1476 TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, min, max);
1477 TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, min, max);
1478 char inputGen[1000];
1479 sprintf(inputGen,"%s>>hGen", variable);
1480 fTree->Draw(inputGen, selection, "groff");
1481 char selectionRec[256];
1482 sprintf(selectionRec, "%s && %s", selection, quality);
1483 char inputRec[1000];
1484 sprintf(inputRec,"%s>>hRec", variable);
1485 fTree->Draw(inputRec, selectionRec, "groff");
1486 //
1487 TH1F* hEff = CreateEffHisto(hGen, hRec);
1488 AliLabelAxes(hEff, variable, "#epsilon [%]");
1489 fRes = hEff;
1490 delete hRec;
1491 delete hGen;
1492 return hEff;
1493}
1494
1495
1496
1497///////////////////////////////////////////////////////////////////////////////////
1498///////////////////////////////////////////////////////////////////////////////////
1499TH1F * AliComparisonDraw::EffLog(const char *variable, const char* selection, const char * quality,
1500 Int_t nbins, Float_t min, Float_t max)
1501{
1502 //
1503 //
1504 Double_t* bins = CreateLogBins(nbins, min, max);
1505 TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, bins);
1506 TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, bins);
1507 char inputGen[1000];
1508 sprintf(inputGen,"%s>>hGen", variable);
1509 fTree->Draw(inputGen, selection, "groff");
1510 char selectionRec[256];
1511 sprintf(selectionRec, "%s && %s", selection, quality);
1512 char inputRec[1000];
1513 sprintf(inputRec,"%s>>hRec", variable);
1514 fTree->Draw(inputRec, selectionRec, "groff");
1515 //
1516 TH1F* hEff = CreateEffHisto(hGen, hRec);
1517 AliLabelAxes(hEff, variable, "#epsilon [%]");
1518 fRes = hEff;
1519 delete hRec;
1520 delete hGen;
1521 delete[] bins;
1522 return hEff;
1523}
1524
1525
1526///////////////////////////////////////////////////////////////////////////////////
1527///////////////////////////////////////////////////////////////////////////////////
1528
1529Double_t* AliComparisonDraw::CreateLogBins(Int_t nBins, Double_t xMin, Double_t xMax)
1530{
1531 Double_t* bins = new Double_t[nBins+1];
1532 bins[0] = xMin;
1533 Double_t factor = pow(xMax/xMin, 1./nBins);
1534 for (Int_t i = 1; i <= nBins; i++)
1535 bins[i] = factor * bins[i-1];
1536 return bins;
1537}
1538
1539
1540
1541
1542void AliComparisonDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
1543{
1544 //
1545 histo->GetXaxis()->SetTitle(xAxisTitle);
1546 histo->GetYaxis()->SetTitle(yAxisTitle);
1547}
1548
1549
1550TH1F* AliComparisonDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec)
1551{
1552 //
1553 Int_t nBins = hGen->GetNbinsX();
1554 TH1F* hEff = (TH1F*) hGen->Clone("hEff");
1555 hEff->SetTitle("");
1556 hEff->SetStats(kFALSE);
1557 hEff->SetMinimum(0.);
1558 hEff->SetMaximum(110.);
1559 //
1560 for (Int_t iBin = 0; iBin <= nBins; iBin++) {
1561 Double_t nGen = hGen->GetBinContent(iBin);
1562 Double_t nRec = hRec->GetBinContent(iBin);
1563 if (nGen > 0) {
1564 Double_t eff = nRec/nGen;
1565 hEff->SetBinContent(iBin, 100. * eff);
1566 Double_t error = sqrt((eff*(1.-eff)+0.01) / nGen);
1567 // if (error == 0) error = sqrt(0.1/nGen);
1568 //
1569 if (error == 0) error = 0.0001;
1570 hEff->SetBinError(iBin, 100. * error);
1571 } else {
1572 hEff->SetBinContent(iBin, 100. * 0.5);
1573 hEff->SetBinError(iBin, 100. * 0.5);
1574 }
1575 }
1576 return hEff;
1577}
1578
1579
1580
1581TH1F* AliComparisonDraw::CreateResHisto(TH2F* hRes2, TH1F **phMean, Bool_t drawBinFits,
1582 Bool_t overflowBinFits)
1583{
1584 TVirtualPad* currentPad = gPad;
1585 TAxis* axis = hRes2->GetXaxis();
1586 Int_t nBins = axis->GetNbins();
1587 TH1F* hRes, *hMean;
1588 if (axis->GetXbins()->GetSize()){
1589 hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
1590 hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
1591 }
1592 else{
1593 hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
1594 hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
1595
1596 }
1597 hRes->SetStats(false);
1598 hRes->SetOption("E");
1599 hRes->SetMinimum(0.);
1600 //
1601 hMean->SetStats(false);
1602 hMean->SetOption("E");
1603
1604 // create the fit function
1605 //TKFitGaus* fitFunc = new TKFitGaus("resFunc");
1606 // TF1 * fitFunc = new TF1("G","[3]+[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
1607 TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
1608
1609 fitFunc->SetLineWidth(2);
1610 fitFunc->SetFillStyle(0);
1611 // create canvas for fits
1612 TCanvas* canBinFits = NULL;
1613 Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
1614 Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
1615 Int_t ny = (nPads-1) / nx + 1;
1616 if (drawBinFits) {
1617 canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
1618 if (canBinFits) delete canBinFits;
1619 canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
1620 canBinFits->Divide(nx, ny);
1621 }
1622
1623 // loop over x bins and fit projection
1624 Int_t dBin = ((overflowBinFits) ? 1 : 0);
1625 for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
1626 if (drawBinFits) canBinFits->cd(bin + dBin);
1627 TH1D* hBin = hRes2->ProjectionY("hBin", bin, bin);
1628 //
1629 if (hBin->GetEntries() > 5) {
1630 fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
1631 hBin->Fit(fitFunc,"s");
1632 Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
1633
1634 if (sigma > 0.){
1635 hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
1636 hMean->SetBinContent(bin, fitFunc->GetParameter(1));
1637 }
1638 else{
1639 hRes->SetBinContent(bin, 0.);
1640 hMean->SetBinContent(bin,0);
1641 }
1642 hRes->SetBinError(bin, fitFunc->GetParError(2));
1643 hMean->SetBinError(bin, fitFunc->GetParError(1));
1644
1645 //
1646 //
1647
1648 } else {
1649 hRes->SetBinContent(bin, 0.);
1650 hRes->SetBinError(bin, 0.);
1651 hMean->SetBinContent(bin, 0.);
1652 hMean->SetBinError(bin, 0.);
1653 }
1654
1655
1656 if (drawBinFits) {
1657 char name[256];
1658 if (bin == 0) {
1659 sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
1660 } else if (bin == nBins+1) {
1661 sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
1662 } else {
1663 sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
1664 axis->GetTitle(), axis->GetBinUpEdge(bin));
1665 }
1666 canBinFits->cd(bin + dBin);
1667 hBin->SetTitle(name);
1668 hBin->SetStats(kTRUE);
1669 hBin->DrawCopy("E");
1670 canBinFits->Update();
1671 canBinFits->Modified();
1672 canBinFits->Update();
1673 }
1674
1675 delete hBin;
1676 }
1677
1678 delete fitFunc;
1679 currentPad->cd();
1680 *phMean = hMean;
1681 return hRes;
1682}
1683
1684
51ad6848 1685void AliComparisonDraw::DrawFriend2D(const char * chx, const char *chy, const char* selection, TTree * tfriend)
1686{
1687
1688}
ae7d73d2 1689
51ad6848 1690
1691void AliComparisonDraw::GetPoints3D(const char * label, const char * chpoints, const char* selection, TTree * tpoints, Int_t color,Float_t rmin){
1692 //
1693 //
1694 if (!fPoints) fPoints = new TObjArray;
1695 if (tpoints->GetIndex()==0) tpoints->BuildIndex("Label","Label");
1696 AliPointsMI * points = new AliPointsMI;
1697 TBranch * br = tpoints->GetBranch(chpoints);
1698 br->SetAddress(&points);
1699 Int_t npoints = fTree->Draw(label,selection);
1700 Float_t xyz[30000];
1701 rmin*=rmin;
1702 for (Int_t i=0;i<npoints;i++){
1703 Int_t index = (Int_t)fTree->GetV1()[i];
1704 tpoints->GetEntryWithIndex(index,index);
1705 if (points->fN<2) continue;
1706 Int_t goodpoints=0;
1707 for (Int_t i=0;i<points->fN; i++){
1708 xyz[goodpoints*3] = points->fX[i];
1709 xyz[goodpoints*3+1] = points->fY[i];
1710 xyz[goodpoints*3+2] = points->fZ[i];
1711 if ( points->fX[i]*points->fX[i]+points->fY[i]*points->fY[i]>rmin) goodpoints++;
1712 }
1713 TPolyMarker3D * marker = new TPolyMarker3D(goodpoints,xyz);
1714 // marker->SeWidth(1);
1715 marker->SetMarkerColor(color);
1716 marker->SetMarkerStyle(1);
1717 fPoints->AddLast(marker);
1718 }
1719}
1720
1721void AliComparisonDraw::Draw3D(Int_t min, Int_t max){
1722 if (!fPoints) return;
1723 Int_t npoints = fPoints->GetEntriesFast();
1724 max = TMath::Min(npoints,max);
1725 min = TMath::Min(npoints,min);
1726
1727 for (Int_t i =min; i<max;i++){
1728 TObject * obj = fPoints->At(i);
1729 if (obj) obj->Draw();
1730 }
1731}
1732
1733void AliComparisonDraw:: SavePoints(const char* name)
1734{
1735 char fname[25];
1736 sprintf(fname,"TH_%s.root",name);
1737 TFile f(fname,"new");
1738 fPoints->Write("Tracks",TObject::kSingleKey);
1739 f.Close();
1740 sprintf(fname,"TH%s.gif",name);
1741 fCanvas->SaveAs(fname);
1742}
1743
1744void AliComparisonDraw::InitView(){
1745 //
1746 //
1747 AliRunLoader* rl = AliRunLoader::Open();
1748 rl->GetEvent(0);
1749 rl->CdGAFile();
1750 //
1751 fCanvas = new TCanvas("cdisplay", "Cluster display",0,0,700,730);
1752 fView = new TView(1);
1753 fView->SetRange(-330,-360,-330, 360,360, 330);
1754 fCanvas->Clear();
1755 fCanvas->SetFillColor(1);
1756 fCanvas->SetTheta(90.);
1757 fCanvas->SetPhi(0.);
1758 //
1759 TGeometry *geom =(TGeometry*)gDirectory->Get("AliceGeom");
1760 geom->Draw("same");
1761
1762}