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