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