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