]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/AliGenInfo.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Obsolet classes removed -renamed (Marian)
[u/mrichter/AliRoot.git] / PWG1 / AliGenInfo.cxx
CommitLineData
c92725b7 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
22Generate complex MC information - used for Comparison later on
23How to use it?
24
022044bf 25gSystem->Load("libPWG1.so")
c92725b7 26AliGenInfoMaker *t = new AliGenInfoMaker("galice.root","genTracks.root",0,0)
27t->Exec();
28
29*/
30
31#if !defined(__CINT__) || defined(__MAKECINT__)
32#include <stdio.h>
33#include <string.h>
34//ROOT includes
35#include "TROOT.h"
36#include "Rtypes.h"
37#include "TFile.h"
38#include "TTree.h"
39#include "TChain.h"
40#include "TCut.h"
41#include "TString.h"
c92725b7 42#include "TStopwatch.h"
43#include "TParticle.h"
44#include "TSystem.h"
c92725b7 45#include "TCanvas.h"
c92725b7 46#include "TGeometry.h"
47#include "TPolyLine3D.h"
48
49//ALIROOT includes
50#include "AliRun.h"
51#include "AliStack.h"
52#include "AliSimDigits.h"
53#include "AliTPCParam.h"
54#include "AliTPC.h"
55#include "AliTPCLoader.h"
56#include "AliDetector.h"
57#include "AliTrackReference.h"
58#include "AliTPCParamSR.h"
59#include "AliTracker.h"
60#include "AliMagF.h"
61#include "AliHelix.h"
c92725b7 62#include "AliTrackPointArray.h"
63
64#endif
65#include "AliGenInfo.h"
66//
67//
68
69ClassImp(AliTPCdigitRow);
70ClassImp(AliMCInfo);
71ClassImp(AliGenV0Info)
72ClassImp(AliGenKinkInfo)
73ClassImp(AliGenInfoMaker)
74
75
76
77AliTPCParam * GetTPCParam(){
78 AliTPCParamSR * par = new AliTPCParamSR;
79 par->Update();
80 return par;
81}
82
83
84//_____________________________________________________________________________
85Float_t TPCBetheBloch(Float_t bg)
86{
87 //
88 // Bethe-Bloch energy loss formula
89 //
90 const Double_t kp1=0.76176e-1;
91 const Double_t kp2=10.632;
92 const Double_t kp3=0.13279e-4;
93 const Double_t kp4=1.8631;
94 const Double_t kp5=1.9479;
95
96 Double_t dbg = (Double_t) bg;
97
98 Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
99
100 Double_t aa = TMath::Power(beta,kp4);
101 Double_t bb = TMath::Power(1./dbg,kp5);
102
103 bb=TMath::Log(kp3+bb);
104
105 return ((Float_t)((kp2-aa-bb)*kp1/aa));
106}
107
108
109
110
111
112////////////////////////////////////////////////////////////////////////
022044bf 113AliMCInfo::AliMCInfo():
114 fTrackRef(),
115 fTrackRefOut(),
116 fTRdecay(),
117 fPrimPart(0),
118 fParticle(),
119 fMass(0),
120 fCharge(0),
121 fLabel(0),
122 fEventNr(),
123 fMCtracks(),
124 fPdg(0),
125 fTPCdecay(0),
126 fRowsWithDigitsInn(0),
127 fRowsWithDigits(0),
128 fRowsTrackLength(0),
129 fPrim(0),
130 fTPCRow(),
131 fNTPCRef(0), // tpc references counter
132 fNITSRef(0), // ITS references counter
133 fNTRDRef(0), // TRD references counter
134 fNTOFRef(0), // TOF references counter
135 fTPCReferences(0),
136 fITSReferences(0),
137 fTRDReferences(0),
138 fTOFReferences(0)
c92725b7 139{
140 fTPCReferences = new TClonesArray("AliTrackReference",10);
141 fITSReferences = new TClonesArray("AliTrackReference",10);
142 fTRDReferences = new TClonesArray("AliTrackReference",10);
143 fTOFReferences = new TClonesArray("AliTrackReference",10);
144 fTRdecay.SetTrack(-1);
145 fCharge = 0;
146}
147
022044bf 148AliMCInfo::AliMCInfo(const AliMCInfo& info):
149 TObject(),
150 fTrackRef(info.fTrackRef),
151 fTrackRefOut(info.fTrackRefOut),
152 fTRdecay(info.fTRdecay),
153 fPrimPart(info.fPrimPart),
154 fParticle(info.fParticle),
155 fMass(info.fMass),
156 fCharge(info.fCharge),
157 fLabel(info.fLabel),
158 fEventNr(info.fEventNr),
159 fMCtracks(info.fMCtracks),
160 fPdg(info.fPdg),
161 fTPCdecay(info.fTPCdecay),
162 fRowsWithDigitsInn(info.fRowsWithDigitsInn),
163 fRowsWithDigits(info.fRowsWithDigits),
164 fRowsTrackLength(info.fRowsTrackLength),
165 fPrim(info.fPrim),
166 fTPCRow(info.fTPCRow),
167 fNTPCRef(info.fNTPCRef), // tpc references counter
168 fNITSRef(info.fNITSRef), // ITS references counter
169 fNTRDRef(info.fNTRDRef), // TRD references counter
170 fNTOFRef(info.fNTOFRef), // TOF references counter
171 fTPCReferences(0),
172 fITSReferences(0),
173 fTRDReferences(0),
174 fTOFReferences(0)
175{
176 fTPCReferences = (TClonesArray*)info.fTPCReferences->Clone();
177 fITSReferences = (TClonesArray*)info.fITSReferences->Clone();
178 fTRDReferences = (TClonesArray*)info.fTRDReferences->Clone();
179 fTOFReferences = (TClonesArray*)info.fTOFReferences->Clone();
180}
181
182
c92725b7 183AliMCInfo::~AliMCInfo()
184{
185 if (fTPCReferences) {
186 delete fTPCReferences;
187 }
188 if (fITSReferences){
189 delete fITSReferences;
190 }
191 if (fTRDReferences){
192 delete fTRDReferences;
193 }
194 if (fTOFReferences){
195 delete fTOFReferences;
196 }
197
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();
215 fNTRDRef = fTRDReferences->GetEntriesFast();
216 fNTOFRef = fTOFReferences->GetEntriesFast();
217
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;
267 fPrim = TPCBetheBloch(betagama);
268 }else fPrim=0;
269 }
270 }else{
271 fMass =0;
272 fPrim =0;
273 }
274}
275
276/////////////////////////////////////////////////////////////////////////////////
022044bf 277AliGenV0Info::AliGenV0Info():
278 fMCd(), //info about daughter particle -
279 fMCm(), //info about mother particle - first particle for V0
280 fMotherP(), //particle info about mother particle
281 fMCDist1(0), //info about closest distance according closest MC - linear DCA
282 fMCDist2(0), //info about closest distance parabolic DCA
283 fMCRr(0), // rec position of the vertex
284 fMCR(0), //exact r position of the vertex
285 fInvMass(0), //reconstructed invariant mass -
286 fPointAngleFi(0), //point angle fi
287 fPointAngleTh(0), //point angle theta
288 fPointAngle(0) //point angle full
c92725b7 289{
022044bf 290 for (Int_t i=0;i<3; i++){
291 fMCPdr[i]=0;
292 fMCX[i]=0;
293 fMCXr[i]=0;
294 fMCPm[i]=0;
295 fMCAngle[i]=0;
296 fMCPd[i]=0;
297 }
298 fMCPd[3]=0;
299 for (Int_t i=0; i<2;i++){
300 fPdg[i]=0;
301 fLab[i]=0;
302 }
c92725b7 303}
c92725b7 304
305void AliGenV0Info::Update(Float_t vertex[3])
306{
022044bf 307 //
308 // Update information - calculates derived variables
309 //
310
c92725b7 311 fMCPd[0] = fMCd.fParticle.Px();
312 fMCPd[1] = fMCd.fParticle.Py();
313 fMCPd[2] = fMCd.fParticle.Pz();
314 fMCPd[3] = fMCd.fParticle.P();
315 //
316 fMCX[0] = fMCd.fParticle.Vx();
317 fMCX[1] = fMCd.fParticle.Vy();
318 fMCX[2] = fMCd.fParticle.Vz();
319 fMCR = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
320 //
321 fPdg[0] = fMCd.fParticle.GetPdgCode();
322 fPdg[1] = fMCm.fParticle.GetPdgCode();
323 //
324 fLab[0] = fMCd.fParticle.GetUniqueID();
325 fLab[1] = fMCm.fParticle.GetUniqueID();
326 //
327 //
328 //
329 Double_t x1[3],p1[3];
330 TParticle & pdaughter = fMCd.fParticle;
331 x1[0] = pdaughter.Vx();
332 x1[1] = pdaughter.Vy();
333 x1[2] = pdaughter.Vz();
334 p1[0] = pdaughter.Px();
335 p1[1] = pdaughter.Py();
336 p1[2] = pdaughter.Pz();
337 Double_t sign = (pdaughter.GetPDG()->Charge()>0)? -1:1;
338 AliHelix dhelix1(x1,p1,sign);
339 //
340 //
341 Double_t x2[3],p2[3];
342 //
343 TParticle & pmother = fMCm.fParticle;
344 x2[0] = pmother.Vx();
345 x2[1] = pmother.Vy();
346 x2[2] = pmother.Vz();
347 p2[0] = pmother.Px();
348 p2[1] = pmother.Py();
349 p2[2] = pmother.Pz();
350 //
351 //
352 sign = (pmother.GetPDG()->Charge()>0) ? -1:1;
353 AliHelix mhelix(x2,p2,sign);
354
355 //
356 //
357 //
358 //find intersection linear
359 //
360 Double_t distance1, distance2;
361 Double_t phase[2][2],radius[2];
362 Int_t points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
363 Double_t delta1=10000,delta2=10000;
364 if (points>0){
365 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
366 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
367 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
368 }
369 else{
370 fInvMass=-1;
371 return;
372 }
373 if (points==2){
374 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
375 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
376 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
377 }
378 distance1 = TMath::Min(delta1,delta2);
379 //
380 //find intersection parabolic
381 //
382 points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
383 delta1=10000,delta2=10000;
384
385 if (points>0){
386 dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
387 }
388 if (points==2){
389 dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
390 }
391
392 distance2 = TMath::Min(delta1,delta2);
393 //
394 if (delta1<delta2){
395 //get V0 info
396 dhelix1.Evaluate(phase[0][0],fMCXr);
397 dhelix1.GetMomentum(phase[0][0],fMCPdr);
398 mhelix.GetMomentum(phase[0][1],fMCPm);
399 dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fMCAngle);
400 fMCRr = TMath::Sqrt(radius[0]);
401 }
402 else{
403 dhelix1.Evaluate(phase[1][0],fMCXr);
404 dhelix1.GetMomentum(phase[1][0],fMCPdr);
405 mhelix.GetMomentum(phase[1][1],fMCPm);
406 dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fMCAngle);
407 fMCRr = TMath::Sqrt(radius[1]);
408 }
409 //
410 //
411 fMCDist1 = TMath::Sqrt(distance1);
412 fMCDist2 = TMath::Sqrt(distance2);
413 //
414 //
415 //
416 Float_t v[3] = {fMCXr[0]-vertex[0],fMCXr[1]-vertex[1],fMCXr[2]-vertex[2]};
417 //Float_t v[3] = {fMCXr[0],fMCXr[1],fMCXr[2]};
418 Float_t p[3] = {fMCPdr[0]+fMCPm[0], fMCPdr[1]+fMCPm[1],fMCPdr[2]+fMCPm[2]};
419 Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
420 Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
421 vnorm2 = TMath::Sqrt(vnorm2);
422 Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
423 Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
424 pnorm2 = TMath::Sqrt(pnorm2);
425 //
426 fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
427 fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);
428 fPointAngle = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
429 Double_t mass1 = fMCd.fMass;
430 Double_t mass2 = fMCm.fMass;
431
432
433 Double_t e1 = TMath::Sqrt(mass1*mass1+
434 fMCPd[0]*fMCPd[0]+
435 fMCPd[1]*fMCPd[1]+
436 fMCPd[2]*fMCPd[2]);
437 Double_t e2 = TMath::Sqrt(mass2*mass2+
438 fMCPm[0]*fMCPm[0]+
439 fMCPm[1]*fMCPm[1]+
440 fMCPm[2]*fMCPm[2]);
441
442 fInvMass =
443 (fMCPm[0]+fMCPd[0])*(fMCPm[0]+fMCPd[0])+
444 (fMCPm[1]+fMCPd[1])*(fMCPm[1]+fMCPd[1])+
445 (fMCPm[2]+fMCPd[2])*(fMCPm[2]+fMCPd[2]);
446
447 // fInvMass = TMath::Sqrt((e1+e2)*(e1+e2)-fInvMass);
448 fInvMass = (e1+e2)*(e1+e2)-fInvMass;
022044bf 449 if (fInvMass>0) fInvMass = TMath::Sqrt(fInvMass);
c92725b7 450}
451
452
453
454/////////////////////////////////////////////////////////////////////////////////
022044bf 455
456AliGenKinkInfo::AliGenKinkInfo():
457 fMCd(), //info about daughter particle - second particle for V0
458 fMCm(), //info about mother particle - first particle for V0
459 fMCDist1(0), //info about closest distance according closest MC - linear DCA
460 fMCDist2(0), //info about closest distance parabolic DCA
461 fMCRr(0), // rec position of the vertex
462 fMCR(0) //exact r position of the vertex
463{
464 //
465 // default constructor
466 //
467 for (Int_t i=0;i<3;i++){
468 fMCPdr[i]=0;
469 fMCPd[i]=0;
470 fMCX[i]=0;
471 fMCPm[i]=0;
472 fMCAngle[i]=0;
473 }
474 for (Int_t i=0; i<2; i++) {
475 fPdg[i]= 0;
476 fLab[i]=0;
477 }
478}
479
c92725b7 480void AliGenKinkInfo::Update()
481{
022044bf 482 //
483 // Update information
484 // some redundancy - faster acces to the values in analysis code
485 //
c92725b7 486 fMCPd[0] = fMCd.fParticle.Px();
487 fMCPd[1] = fMCd.fParticle.Py();
488 fMCPd[2] = fMCd.fParticle.Pz();
489 fMCPd[3] = fMCd.fParticle.P();
490 //
491 fMCX[0] = fMCd.fParticle.Vx();
492 fMCX[1] = fMCd.fParticle.Vy();
493 fMCX[2] = fMCd.fParticle.Vz();
494 fMCR = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
495 //
496 fPdg[0] = fMCd.fParticle.GetPdgCode();
497 fPdg[1] = fMCm.fParticle.GetPdgCode();
498 //
499 fLab[0] = fMCd.fParticle.GetUniqueID();
500 fLab[1] = fMCm.fParticle.GetUniqueID();
501 //
502 //
503 //
504 Double_t x1[3],p1[3];
505 TParticle & pdaughter = fMCd.fParticle;
506 x1[0] = pdaughter.Vx();
507 x1[1] = pdaughter.Vy();
508 x1[2] = pdaughter.Vz();
509 p1[0] = pdaughter.Px();
510 p1[1] = pdaughter.Py();
511 p1[2] = pdaughter.Pz();
512 Double_t sign = (pdaughter.GetPDG()->Charge()>0)? -1:1;
513 AliHelix dhelix1(x1,p1,sign);
514 //
515 //
516 Double_t x2[3],p2[3];
517 //
518 TParticle & pmother = fMCm.fParticle;
519 x2[0] = pmother.Vx();
520 x2[1] = pmother.Vy();
521 x2[2] = pmother.Vz();
522 p2[0] = pmother.Px();
523 p2[1] = pmother.Py();
524 p2[2] = pmother.Pz();
525 //
526 AliTrackReference & pdecay = fMCm.fTRdecay;
527 x2[0] = pdecay.X();
528 x2[1] = pdecay.Y();
529 x2[2] = pdecay.Z();
530 p2[0] = pdecay.Px();
531 p2[1] = pdecay.Py();
532 p2[2] = pdecay.Pz();
533 //
534 sign = (pmother.GetPDG()->Charge()>0) ? -1:1;
535 AliHelix mhelix(x2,p2,sign);
536
537 //
538 //
539 //
540 //find intersection linear
541 //
542 Double_t distance1, distance2;
543 Double_t phase[2][2],radius[2];
544 Int_t points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
545 Double_t delta1=10000,delta2=10000;
546 if (points>0){
547 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
548 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
549 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
550 }
551 if (points==2){
552 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
553 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
554 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
555 }
556 distance1 = TMath::Min(delta1,delta2);
557 //
558 //find intersection parabolic
559 //
560 points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
561 delta1=10000,delta2=10000;
562
563 if (points>0){
564 dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
565 }
566 if (points==2){
567 dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
568 }
569
570 distance2 = TMath::Min(delta1,delta2);
571 //
572 if (delta1<delta2){
573 //get V0 info
574 dhelix1.Evaluate(phase[0][0],fMCXr);
575 dhelix1.GetMomentum(phase[0][0],fMCPdr);
576 mhelix.GetMomentum(phase[0][1],fMCPm);
577 dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fMCAngle);
578 fMCRr = TMath::Sqrt(radius[0]);
579 }
580 else{
581 dhelix1.Evaluate(phase[1][0],fMCXr);
582 dhelix1.GetMomentum(phase[1][0],fMCPdr);
583 mhelix.GetMomentum(phase[1][1],fMCPm);
584 dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fMCAngle);
585 fMCRr = TMath::Sqrt(radius[1]);
586 }
587 //
588 //
589 fMCDist1 = TMath::Sqrt(distance1);
590 fMCDist2 = TMath::Sqrt(distance2);
591
592}
593
594
595Float_t AliGenKinkInfo::GetQt(){
596 //
597 //
598 Float_t momentumd = TMath::Sqrt(fMCPd[0]*fMCPd[0]+fMCPd[1]*fMCPd[1]+fMCPd[2]*fMCPd[2]);
599 return TMath::Sin(fMCAngle[2])*momentumd;
600}
601
602
603
c92725b7 604
605////////////////////////////////////////////////////////////////////////
606AliTPCdigitRow::AliTPCdigitRow()
607{
608 Reset();
609}
610////////////////////////////////////////////////////////////////////////
611AliTPCdigitRow & AliTPCdigitRow::operator=(const AliTPCdigitRow &digOld)
612{
613 for (Int_t i = 0; i<kgRowBytes; i++) fDig[i] = digOld.fDig[i];
614 return (*this);
615}
616////////////////////////////////////////////////////////////////////////
617void AliTPCdigitRow::SetRow(Int_t row)
618{
619 if (row >= 8*kgRowBytes) {
620 cerr<<"AliTPCdigitRow::SetRow: index "<<row<<" out of bounds."<<endl;
621 return;
622 }
623 Int_t iC = row/8;
624 Int_t iB = row%8;
625 SETBIT(fDig[iC],iB);
626}
627
628////////////////////////////////////////////////////////////////////////
022044bf 629Bool_t AliTPCdigitRow::TestRow(Int_t row) const
c92725b7 630{
631//
632// return kTRUE if row is on
633//
634 Int_t iC = row/8;
635 Int_t iB = row%8;
636 return TESTBIT(fDig[iC],iB);
637}
638////////////////////////////////////////////////////////////////////////
022044bf 639Int_t AliTPCdigitRow::RowsOn(Int_t upto) const
c92725b7 640{
641//
642// returns number of rows with a digit
643// count only rows less equal row number upto
644//
645 Int_t total = 0;
646 for (Int_t i = 0; i<kgRowBytes; i++) {
647 for (Int_t j = 0; j < 8; j++) {
648 if (i*8+j > upto) return total;
649 if (TESTBIT(fDig[i],j)) total++;
650 }
651 }
652 return total;
653}
654////////////////////////////////////////////////////////////////////////
655void AliTPCdigitRow::Reset()
656{
657//
658// resets all rows to zero
659//
660 for (Int_t i = 0; i<kgRowBytes; i++) {
661 fDig[i] <<= 8;
662 }
663}
664////////////////////////////////////////////////////////////////////////
022044bf 665Int_t AliTPCdigitRow::Last() const
c92725b7 666{
667//
668// returns the last row number with a digit
669// returns -1 if now digits
670//
671 for (Int_t i = kgRowBytes-1; i>=0; i--) {
672 for (Int_t j = 7; j >= 0; j--) {
673 if TESTBIT(fDig[i],j) return i*8+j;
674 }
675 }
676 return -1;
677}
678////////////////////////////////////////////////////////////////////////
022044bf 679Int_t AliTPCdigitRow::First() const
c92725b7 680{
681//
682// returns the first row number with a digit
683// returns -1 if now digits
684//
685 for (Int_t i = 0; i<kgRowBytes; i++) {
686 for (Int_t j = 0; j < 8; j++) {
687 if (TESTBIT(fDig[i],j)) return i*8+j;
688 }
689 }
690 return -1;
691}
692
022044bf 693
c92725b7 694
695////////////////////////////////////////////////////////////////////////
022044bf 696AliGenInfoMaker::AliGenInfoMaker():
697 fDebug(0), //! debug flag
698 fEventNr(0), //! current event number
699 fLabel(0), //! track label
700 fNEvents(0), //! number of events to process
701 fFirstEventNr(0), //! first event to process
702 fNParticles(0), //! number of particles in TreeK
703 fTreeGenTracks(0), //! output tree with generated tracks
704 fTreeKinks(0), //! output tree with Kinks
705 fTreeV0(0), //! output tree with V0
706 fTreeHitLines(0), //! tree with hit lines
707 fFileGenTracks(0), //! output file with stored fTreeGenTracks
708 fLoader(0), //! pointer to the run loader
709 fTreeD(0), //! current tree with digits
710 fTreeTR(0), //! current tree with TR
711 fStack(0), //! current stack
712 fGenInfo(0), //! array with pointers to gen info
713 fNInfos(0), //! number of tracks with infos
714 fParamTPC(0), //! AliTPCParam
715 fTPCPtCut(0.03),
716 fITSPtCut(0.1),
717 fTRDPtCut(0.1),
718 fTOFPtCut(0.1)
719{
c92725b7 720}
721
722////////////////////////////////////////////////////////////////////////
723AliGenInfoMaker::AliGenInfoMaker(const char * fnGalice, const char* fnRes,
022044bf 724 Int_t nEvents, Int_t firstEvent):
725 fDebug(0), //! debug flag
726 fEventNr(0), //! current event number
727 fLabel(0), //! track label
728 fNEvents(0), //! number of events to process
729 fFirstEventNr(0), //! first event to process
730 fNParticles(0), //! number of particles in TreeK
731 fTreeGenTracks(0), //! output tree with generated tracks
732 fTreeKinks(0), //! output tree with Kinks
733 fTreeV0(0), //! output tree with V0
734 fTreeHitLines(0), //! tree with hit lines
735 fFileGenTracks(0), //! output file with stored fTreeGenTracks
736 fLoader(0), //! pointer to the run loader
737 fTreeD(0), //! current tree with digits
738 fTreeTR(0), //! current tree with TR
739 fStack(0), //! current stack
740 fGenInfo(0), //! array with pointers to gen info
741 fNInfos(0), //! number of tracks with infos
742 fParamTPC(0), //! AliTPCParam
743 fTPCPtCut(0.03),
744 fITSPtCut(0.1),
745 fTRDPtCut(0.1),
746 fTOFPtCut(0.1)
747
c92725b7 748{
022044bf 749 //
750 //
751 //
c92725b7 752 fFirstEventNr = firstEvent;
753 fEventNr = firstEvent;
754 fNEvents = nEvents;
c92725b7 755 sprintf(fFnRes,"%s",fnRes);
756 //
757 fLoader = AliRunLoader::Open(fnGalice);
758 if (gAlice){
759 delete gAlice->GetRunLoader();
760 delete gAlice;
761 gAlice = 0x0;
762 }
763 if (fLoader->LoadgAlice()){
764 cerr<<"Error occured while l"<<endl;
765 }
766 Int_t nall = fLoader->GetNumberOfEvents();
767 if (nEvents==0) {
768 nEvents =nall;
769 fNEvents=nall;
770 fFirstEventNr=0;
771 }
772
773 if (nall<=0){
774 cerr<<"no events available"<<endl;
775 fEventNr = 0;
776 return;
777 }
778 if (firstEvent+nEvents>nall) {
779 fEventNr = nall-firstEvent;
780 cerr<<"restricted number of events availaible"<<endl;
781 }
782 AliMagF * magf = gAlice->Field();
783 AliTracker::SetFieldMap(magf,0);
784}
785
786
787AliMCInfo * AliGenInfoMaker::MakeInfo(UInt_t i)
788{
789 //
790 if (i<fNParticles) {
791 if (fGenInfo[i]) return fGenInfo[i];
792 fGenInfo[i] = new AliMCInfo;
793 fNInfos++;
794 return fGenInfo[i];
795 }
796 else
797 return 0;
798}
799
c92725b7 800////////////////////////////////////////////////////////////////////////
801AliGenInfoMaker::~AliGenInfoMaker()
802{
803
804 if (fLoader){
805 fLoader->UnloadgAlice();
806 gAlice = 0;
807 delete fLoader;
808 }
809}
810
811Int_t AliGenInfoMaker::SetIO()
812{
813 //
814 //
815 CreateTreeGenTracks();
816 if (!fTreeGenTracks) return 1;
817 // AliTracker::SetFieldFactor();
818
819 fParamTPC = GetTPCParam();
820 //
821 return 0;
822}
823
824////////////////////////////////////////////////////////////////////////
825Int_t AliGenInfoMaker::SetIO(Int_t eventNr)
826{
827 //
828 //
829 // SET INPUT
830 fLoader->SetEventNumber(eventNr);
831 //
832 fLoader->LoadHeader();
833 fLoader->LoadKinematics();
834 fStack = fLoader->Stack();
835 //
836 fLoader->LoadTrackRefs();
837 fLoader->LoadHits();
838 fTreeTR = fLoader->TreeTR();
839 //
840 AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
841 tpcl->LoadDigits();
842 fTreeD = tpcl->TreeD();
843 return 0;
844}
845
846Int_t AliGenInfoMaker::CloseIOEvent()
847{
848 fLoader->UnloadHeader();
849 fLoader->UnloadKinematics();
850 fLoader->UnloadTrackRefs();
851 AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
852 tpcl->UnloadDigits();
853 return 0;
854}
855
856Int_t AliGenInfoMaker::CloseIO()
857{
858 fLoader->UnloadgAlice();
859 return 0;
860}
861
862
863
864////////////////////////////////////////////////////////////////////////
865Int_t AliGenInfoMaker::Exec(Int_t nEvents, Int_t firstEventNr)
866{
867 fNEvents = nEvents;
868 fFirstEventNr = firstEventNr;
869 return Exec();
870}
871
872////////////////////////////////////////////////////////////////////////
873Int_t AliGenInfoMaker::Exec()
874{
875 TStopwatch timer;
876 timer.Start();
877 Int_t status =SetIO();
878 if (status>0) return status;
879 //
880
881 for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
882 fEventNr++) {
883 SetIO(fEventNr);
884 fNParticles = fStack->GetNtrack();
885 //
886 fGenInfo = new AliMCInfo*[fNParticles];
887 for (UInt_t i = 0; i<fNParticles; i++) {
888 fGenInfo[i]=0;
889 }
890 //
891 cout<<"Start to process event "<<fEventNr<<endl;
892 cout<<"\tfNParticles = "<<fNParticles<<endl;
893 if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeTR"<<endl;
894 if (TreeTRLoop()>0) return 1;
895 //
896 if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeD"<<endl;
897 if (TreeDLoop()>0) return 1;
898 //
899 if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeK"<<endl;
900 if (TreeKLoop()>0) return 1;
901 if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
902 //
903 if (BuildKinkInfo()>0) return 1;
904 if (BuildV0Info()>0) return 1;
905 //if (BuildHitLines()>0) return 1;
906 if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
907 //
908 for (UInt_t i = 0; i<fNParticles; i++) {
909 if (fGenInfo[i]) delete fGenInfo[i];
910 }
911 delete []fGenInfo;
912 CloseIOEvent();
913 }
914 //
915 CloseIO();
916 CloseOutputFile();
917
918 cerr<<"Exec finished"<<endl;
919
920 timer.Stop();
921 timer.Print();
922 return 0;
923}
924
925////////////////////////////////////////////////////////////////////////
926void AliGenInfoMaker::CreateTreeGenTracks()
927{
928 fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
929 if (!fFileGenTracks) {
930 cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
931 return;
932 }
933 fTreeGenTracks = new TTree("genTracksTree","genTracksTree");
934 AliMCInfo * info = new AliMCInfo;
935 fTreeGenTracks->Branch("MC","AliMCInfo",&info);
936 delete info;
937 //
938 AliGenKinkInfo *kinkinfo = new AliGenKinkInfo;
939 fTreeKinks = new TTree("genKinksTree","genKinksTree");
940 fTreeKinks->Branch("MC","AliGenKinkInfo",&kinkinfo);
941 delete kinkinfo;
942 //
943 AliGenV0Info *v0info = new AliGenV0Info;
944 fTreeV0 = new TTree("genV0Tree","genV0Tree");
945 fTreeV0->Branch("MC","AliGenV0Info",&v0info);
946 delete v0info;
947 //
948 //
949 AliTrackPointArray * points0 = new AliTrackPointArray;
950 AliTrackPointArray * points1 = new AliTrackPointArray;
951 AliTrackPointArray * points2 = new AliTrackPointArray;
952 fTreeHitLines = new TTree("HitLines","HitLines");
953 fTreeHitLines->Branch("TPC.","AliTrackPointArray",&points0,32000,10);
954 fTreeHitLines->Branch("TRD.","AliTrackPointArray",&points1,32000,10);
955 fTreeHitLines->Branch("ITS.","AliTrackPointArray",&points2,32000,10);
956 Int_t label=0;
957 fTreeHitLines->Branch("Label",&label,"label/I");
958 //
959 fTreeGenTracks->AutoSave();
960 fTreeKinks->AutoSave();
961 fTreeV0->AutoSave();
962 fTreeHitLines->AutoSave();
963}
964////////////////////////////////////////////////////////////////////////
965void AliGenInfoMaker::CloseOutputFile()
966{
967 if (!fFileGenTracks) {
968 cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
969 return;
970 }
971 fFileGenTracks->cd();
972 fTreeGenTracks->Write();
973 delete fTreeGenTracks;
974 fTreeKinks->Write();
975 delete fTreeKinks;
976 fTreeV0->Write();
977 delete fTreeV0;
978
979 fFileGenTracks->Close();
980 delete fFileGenTracks;
981 return;
982}
983
984////////////////////////////////////////////////////////////////////////
985Int_t AliGenInfoMaker::TreeKLoop()
986{
987//
988// open the file with treeK
989// loop over all entries there and save information about some tracks
990//
991
992 AliStack * stack = fStack;
993 if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
994
995 if (fDebug > 0) {
996 cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
997 <<fEventNr<<endl;
998 }
999 Int_t ipdg = 0;
1000 TParticlePDG *ppdg = 0;
1001 // not all generators give primary vertex position. Take the vertex
1002 // of the particle 0 as primary vertex.
1003 TDatabasePDG pdg; //get pdg table
1004 //thank you very much root for this
1005 TBranch * br = fTreeGenTracks->GetBranch("MC");
1006 TParticle *particle = stack->ParticleFromTreeK(0);
1007 fVPrim[0] = particle->Vx();
1008 fVPrim[1] = particle->Vy();
1009 fVPrim[2] = particle->Vz();
1010 for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
1011 // load only particles with TR
1012 AliMCInfo * info = GetInfo(iParticle);
1013 if (!info) continue;
1014 //////////////////////////////////////////////////////////////////////
1015 info->fLabel = iParticle;
1016 //
1017 info->fParticle = *(stack->Particle(iParticle));
1018 info->fVDist[0] = info->fParticle.Vx()-fVPrim[0];
1019 info->fVDist[1] = info->fParticle.Vy()-fVPrim[1];
1020 info->fVDist[2] = info->fParticle.Vz()-fVPrim[2];
1021 info->fVDist[3] = TMath::Sqrt(info->fVDist[0]*info->fVDist[0]+
1022 info->fVDist[1]*info->fVDist[1]+info->fVDist[2]*info->fVDist[2]);
1023 //
1024 //
1025 ipdg = info->fParticle.GetPdgCode();
1026 info->fPdg = ipdg;
1027 ppdg = pdg.GetParticle(ipdg);
1028 info->fEventNr = fEventNr;
1029 info->Update();
1030 //////////////////////////////////////////////////////////////////////
1031 br->SetAddress(&info);
1032 fTreeGenTracks->Fill();
1033 }
1034 fTreeGenTracks->AutoSave();
1035 if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
1036 return 0;
1037}
1038
1039
1040////////////////////////////////////////////////////////////////////////////////////
1041Int_t AliGenInfoMaker::BuildKinkInfo()
1042{
1043 //
1044 // Build tree with Kink Information
1045 //
1046 AliStack * stack = fStack;
1047 if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
1048
1049 if (fDebug > 0) {
1050 cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
1051 <<fEventNr<<endl;
1052 }
1053 // Int_t ipdg = 0;
1054 //TParticlePDG *ppdg = 0;
1055 // not all generators give primary vertex position. Take the vertex
1056 // of the particle 0 as primary vertex.
1057 TDatabasePDG pdg; //get pdg table
1058 //thank you very much root for this
1059 TBranch * br = fTreeKinks->GetBranch("MC");
1060 //
1061 AliGenKinkInfo * kinkinfo = new AliGenKinkInfo;
1062 //
1063 for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
1064 // load only particles with TR
1065 AliMCInfo * info = GetInfo(iParticle);
1066 if (!info) continue;
1067 if (info->fCharge==0) continue;
1068 if (info->fTRdecay.P()<0.13) continue; //momenta cut
1069 if (info->fTRdecay.R()>500) continue; //R cut - decay outside barrel
1070 TParticle & particle = info->fParticle;
1071 Int_t first = particle.GetDaughter(0);
1072 if (first ==0) continue;
1073
1074 Int_t last = particle.GetDaughter(1);
1075 if (last ==0) last=first;
1076 AliMCInfo * info2 = 0;
1077 AliMCInfo * dinfo = 0;
1078
1079 for (Int_t id2=first;id2<=last;id2++){
1080 info2 = GetInfo(id2);
1081 if (!info2) continue;
1082 TParticle &p2 = info2->fParticle;
1083 Double_t r = TMath::Sqrt(p2.Vx()*p2.Vx()+p2.Vy()*p2.Vy());
1084 if ( TMath::Abs(info->fTRdecay.R()-r)>0.01) continue;
1085 if (!(p2.GetPDG())) continue;
1086 dinfo =info2;
1087
1088 kinkinfo->fMCm = (*info);
1089 kinkinfo->fMCd = (*dinfo);
1090 if (kinkinfo->fMCm.fParticle.GetPDG()==0 || kinkinfo->fMCd.fParticle.GetPDG()==0) continue;
1091 kinkinfo->Update();
1092 br->SetAddress(&kinkinfo);
1093 fTreeKinks->Fill();
1094 }
1095 /*
1096 if (dinfo){
1097 kinkinfo->fMCm = (*info);
1098 kinkinfo->fMCd = (*dinfo);
1099 kinkinfo->Update();
1100 br->SetAddress(&kinkinfo);
1101 fTreeKinks->Fill();
1102 }
1103 */
1104 }
1105 fTreeGenTracks->AutoSave();
1106 if (fDebug > 2) cerr<<"end of Kink Loop"<<endl;
1107 return 0;
1108}
1109
1110
1111////////////////////////////////////////////////////////////////////////////////////
1112Int_t AliGenInfoMaker::BuildV0Info()
1113{
1114 //
1115 // Build tree with V0 Information
1116 //
1117 AliStack * stack = fStack;
1118 if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
1119
1120 if (fDebug > 0) {
1121 cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
1122 <<fEventNr<<endl;
1123 }
1124 // Int_t ipdg = 0;
1125 //TParticlePDG *ppdg = 0;
1126 // not all generators give primary vertex position. Take the vertex
1127 // of the particle 0 as primary vertex.
1128 TDatabasePDG pdg; //get pdg table
1129 //thank you very much root for this
1130 TBranch * br = fTreeV0->GetBranch("MC");
1131 //
1132 AliGenV0Info * v0info = new AliGenV0Info;
1133 //
1134 //
1135 for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
1136 // load only particles with TR
1137 AliMCInfo * info = GetInfo(iParticle);
1138 if (!info) continue;
1139 if (info->fCharge==0) continue;
1140 //
1141 //
1142 TParticle & particle = info->fParticle;
1143 Int_t mother = particle.GetMother(0);
1144 if (mother <=0) continue;
1145 //
1146 TParticle * motherparticle = stack->Particle(mother);
1147 if (!motherparticle) continue;
1148 //
1149 Int_t last = motherparticle->GetDaughter(1);
1150 if (last==(int)iParticle) continue;
1151 AliMCInfo * info2 = info;
1152 AliMCInfo * dinfo = GetInfo(last);
1153 if (!dinfo) continue;
1154 if (!dinfo->fParticle.GetPDG()) continue;
1155 if (!info2->fParticle.GetPDG()) continue;
1156 if (dinfo){
1157 v0info->fMCm = (*info);
1158 v0info->fMCd = (*dinfo);
1159 v0info->fMotherP = (*motherparticle);
1160 v0info->Update(fVPrim);
1161 br->SetAddress(&v0info);
1162 fTreeV0->Fill();
1163 }
1164 }
1165 fTreeV0->AutoSave();
1166 if (fDebug > 2) cerr<<"end of V0 Loop"<<endl;
1167 return 0;
1168}
1169
1170
1171
1172
1173////////////////////////////////////////////////////////////////////////
1174Int_t AliGenInfoMaker::BuildHitLines()
1175{
1176
1177//
1178// open the file with treeK
1179// loop over all entries there and save information about some tracks
1180//
1181
1182 AliStack * stack = fStack;
1183 if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
1184
1185 if (fDebug > 0) {
1186 cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
1187 <<fEventNr<<endl;
1188 }
1189// Int_t ipdg = 0;
1190// // TParticlePDG *ppdg = 0;
1191// // not all generators give primary vertex position. Take the vertex
1192// // of the particle 0 as primary vertex.
1193// TDatabasePDG pdg; //get pdg table
1194// //thank you very much root for this
1195// AliTrackPointArray *tpcp = new AliTrackPointArray;
1196// AliTrackPointArray *trdp = new AliTrackPointArray;
1197// AliTrackPointArray *itsp = new AliTrackPointArray;
1198// Int_t label =0;
1199// //
1200// TBranch * brtpc = fTreeHitLines->GetBranch("TPC.");
1201// TBranch * brtrd = fTreeHitLines->GetBranch("TRD.");
1202// TBranch * brits = fTreeHitLines->GetBranch("ITS.");
1203// TBranch * brlabel = fTreeHitLines->GetBranch("Label");
1204// brlabel->SetAddress(&label);
1205// brtpc->SetAddress(&tpcp);
1206// brtrd->SetAddress(&trdp);
1207// brits->SetAddress(&itsp);
1208// //
1209// AliDetector *dtpc = gAlice->GetDetector("TPC");
1210// AliDetector *dtrd = gAlice->GetDetector("TRD");
1211// AliDetector *dits = gAlice->GetDetector("ITS");
1212
1213// for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
1214// // load only particles with TR
1215// AliMCInfo * info = GetInfo(iParticle);
1216// if (!info) continue;
1217// Int_t primpart = info->fPrimPart;
1218// ipdg = info->fParticle.GetPdgCode();
1219// label = iParticle;
1220// //
1221// gAlice->ResetHits();
1222// tpcp->Reset();
1223// itsp->Reset();
1224// trdp->Reset();
1225// tpcp->fLabel1 = ipdg;
1226// trdp->fLabel1 = ipdg;
1227// itsp->fLabel1 = ipdg;
1228// if (dtpc->TreeH()->GetEvent(primpart)){
1229// dtpc->LoadPoints(primpart);
1230// tpcp->Reset(dtpc,iParticle);
1231// }
1232// if (dtrd->TreeH()->GetEvent(primpart)){
1233// dtrd->LoadPoints(primpart);
1234// trdp->Reset(dtrd,iParticle);
1235// }
1236// if (dits->TreeH()->GetEvent(primpart)){
1237// dits->LoadPoints(primpart);
1238// itsp->Reset(dits,iParticle);
1239// }
1240// //
1241// fTreeHitLines->Fill();
1242// dtpc->ResetPoints();
1243// dtrd->ResetPoints();
1244// dits->ResetPoints();
1245// }
1246// delete tpcp;
1247// delete trdp;
1248// delete itsp;
1249// fTreeHitLines->AutoSave();
1250// if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
1251 return 0;
1252}
1253
1254
1255////////////////////////////////////////////////////////////////////////
1256Int_t AliGenInfoMaker::TreeDLoop()
1257{
1258 //
1259 // open the file with treeD
1260 // loop over all entries there and save information about some tracks
1261 //
1262
1263 Int_t nInnerSector = fParamTPC->GetNInnerSector();
1264 Int_t rowShift = 0;
1265 Int_t zero=fParamTPC->GetZeroSup()+6;
1266 //
1267 //
1268 AliSimDigits digitsAddress, *digits=&digitsAddress;
1269 fTreeD->GetBranch("Segment")->SetAddress(&digits);
1270
1271 Int_t sectorsByRows=(Int_t)fTreeD->GetEntries();
1272 if (fDebug > 1) cout<<"\tsectorsByRows = "<<sectorsByRows<<endl;
1273 for (Int_t i=0; i<sectorsByRows; i++) {
1274 if (!fTreeD->GetEvent(i)) continue;
1275 Int_t sec,row;
1276 fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
1277 if (fDebug > 1) cout<<sec<<' '<<row<<" \r";
1278 // here I expect that upper sectors follow lower sectors
1279 if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow();
1280 //
1281 digits->ExpandTrackBuffer();
1282 digits->First();
1283 do {
1284 Int_t iRow=digits->CurrentRow();
1285 Int_t iColumn=digits->CurrentColumn();
1286 Short_t digitValue = digits->CurrentDigit();
1287 if (digitValue >= zero) {
1288 Int_t label;
1289 for (Int_t j = 0; j<3; j++) {
1290 // label = digits->GetTrackID(iRow,iColumn,j);
1291 label = digits->GetTrackIDFast(iRow,iColumn,j)-2;
1292 if (label >= (Int_t)fNParticles) { //don't label from bakground event
1293 continue;
1294 }
1295 if (label >= 0 && label <= (Int_t)fNParticles) {
1296 if (fDebug > 6 ) {
1297 cout<<"Inner loop: sector, iRow, iColumn, label, value, row "
1298 <<sec<<" "
1299 <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
1300 <<" "<<row<<endl;
1301 }
1302 AliMCInfo * info = GetInfo(label);
1303 if (info){
1304 info->fTPCRow.SetRow(row+rowShift);
1305 }
1306 }
1307 }
1308 }
1309 } while (digits->Next());
1310 }
1311
1312 if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;
1313 return 0;
1314}
1315
1316
1317////////////////////////////////////////////////////////////////////////
1318Int_t AliGenInfoMaker::TreeTRLoop()
1319{
1320 //
1321 // loop over TrackReferences and store the first one for each track
1322 //
1323 TTree * treeTR = fTreeTR;
1324 Int_t nPrimaries = (Int_t) treeTR->GetEntries();
1325 if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl;
1326 //
1327 //
1328 //track references for TPC
022044bf 1329 TClonesArray* tpcArrayTR = new TClonesArray("AliTrackReference");
1330 TClonesArray* itsArrayTR = new TClonesArray("AliTrackReference");
1331 TClonesArray* trdArrayTR = new TClonesArray("AliTrackReference");
1332 TClonesArray* tofArrayTR = new TClonesArray("AliTrackReference");
1333 TClonesArray* runArrayTR = new TClonesArray("AliTrackReference");
c92725b7 1334 //
022044bf 1335 if (treeTR->GetBranch("TPC")) treeTR->GetBranch("TPC")->SetAddress(&tpcArrayTR);
1336 if (treeTR->GetBranch("ITS")) treeTR->GetBranch("ITS")->SetAddress(&itsArrayTR);
1337 if (treeTR->GetBranch("TRD")) treeTR->GetBranch("TRD")->SetAddress(&trdArrayTR);
1338 if (treeTR->GetBranch("TOF")) treeTR->GetBranch("TOF")->SetAddress(&tofArrayTR);
1339 if (treeTR->GetBranch("AliRun")) treeTR->GetBranch("AliRun")->SetAddress(&runArrayTR);
c92725b7 1340 //
1341 //
1342 //
1343 for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
1344 treeTR->GetEntry(iPrimPart);
1345 //
1346 // Loop over TPC references
1347 //
022044bf 1348 for (Int_t iTrackRef = 0; iTrackRef < tpcArrayTR->GetEntriesFast(); iTrackRef++) {
1349 AliTrackReference *trackRef = (AliTrackReference*)tpcArrayTR->At(iTrackRef);
c92725b7 1350 //
1351 if (trackRef->TestBit(BIT(2))){
1352 //if decay
022044bf 1353 if (trackRef->P()<fTPCPtCut) continue;
c92725b7 1354 Int_t label = trackRef->GetTrack();
1355 AliMCInfo * info = GetInfo(label);
1356 if (!info) info = MakeInfo(label);
1357 info->fTRdecay = *trackRef;
1358 }
1359 //
022044bf 1360 if (trackRef->P()<fTPCPtCut) continue;
c92725b7 1361 Int_t label = trackRef->GetTrack();
1362 AliMCInfo * info = GetInfo(label);
1363 if (!info) info = MakeInfo(label);
1364 if (!info) continue;
1365 info->fPrimPart = iPrimPart;
1366 TClonesArray & arr = *(info->fTPCReferences);
1367 new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
1368 }
1369 //
1370 // Loop over ITS references
1371 //
022044bf 1372 for (Int_t iTrackRef = 0; iTrackRef < itsArrayTR->GetEntriesFast(); iTrackRef++) {
1373 AliTrackReference *trackRef = (AliTrackReference*)itsArrayTR->At(iTrackRef);
c92725b7 1374 //
1375 //
022044bf 1376 if (trackRef->P()<fTPCPtCut) continue;
c92725b7 1377 Int_t label = trackRef->GetTrack();
1378 AliMCInfo * info = GetInfo(label);
022044bf 1379 if ( (!info) && trackRef->Pt()<fITSPtCut) continue;
c92725b7 1380 if (!info) info = MakeInfo(label);
1381 if (!info) continue;
1382 info->fPrimPart = iPrimPart;
1383 TClonesArray & arr = *(info->fITSReferences);
1384 new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
1385 }
1386 //
1387 // Loop over TRD references
1388 //
022044bf 1389 for (Int_t iTrackRef = 0; iTrackRef < trdArrayTR->GetEntriesFast(); iTrackRef++) {
1390 AliTrackReference *trackRef = (AliTrackReference*)trdArrayTR->At(iTrackRef);
c92725b7 1391 //
022044bf 1392 if (trackRef->P()<fTPCPtCut) continue;
c92725b7 1393 Int_t label = trackRef->GetTrack();
1394 AliMCInfo * info = GetInfo(label);
022044bf 1395 if ( (!info) && trackRef->Pt()<fTRDPtCut) continue;
c92725b7 1396 if (!info) info = MakeInfo(label);
1397 if (!info) continue;
1398 info->fPrimPart = iPrimPart;
1399 TClonesArray & arr = *(info->fTRDReferences);
1400 new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
1401 }
1402 //
1403 // Loop over TOF references
1404 //
022044bf 1405 for (Int_t iTrackRef = 0; iTrackRef < tofArrayTR->GetEntriesFast(); iTrackRef++) {
1406 AliTrackReference *trackRef = (AliTrackReference*)tofArrayTR->At(iTrackRef);
c92725b7 1407 Int_t label = trackRef->GetTrack();
1408 AliMCInfo * info = GetInfo(label);
1409 if (!info){
022044bf 1410 if (trackRef->Pt()<fTPCPtCut) continue;
1411 if ( (!info) && trackRef->Pt()<fTOFPtCut) continue;
c92725b7 1412 }
1413 if (!info) info = MakeInfo(label);
1414 if (!info) continue;
1415 info->fPrimPart = iPrimPart;
1416 TClonesArray & arr = *(info->fTOFReferences);
1417 new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
1418 }
1419 //
1420 // get dacay position
022044bf 1421 for (Int_t iTrackRef = 0; iTrackRef < runArrayTR->GetEntriesFast(); iTrackRef++) {
1422 AliTrackReference *trackRef = (AliTrackReference*)runArrayTR->At(iTrackRef);
c92725b7 1423 //
1424 Int_t label = trackRef->GetTrack();
1425 AliMCInfo * info = GetInfo(label);
1426 if (!info) continue;
1427 if (!trackRef->TestBit(BIT(2))) continue; //if not decay
1428 // if (TMath::Abs(trackRef.X());
1429 info->fTRdecay = *trackRef;
1430 }
1431 }
1432 //
022044bf 1433 tpcArrayTR->Delete();
1434 delete tpcArrayTR;
1435 trdArrayTR->Delete();
1436 delete trdArrayTR;
1437 tofArrayTR->Delete();
1438 delete tofArrayTR;
1439 itsArrayTR->Delete();
1440 delete itsArrayTR;
1441 runArrayTR->Delete();
1442 delete runArrayTR;
c92725b7 1443 //
1444 return 0;
1445}
1446
1447////////////////////////////////////////////////////////////////////////
1448Float_t AliGenInfoMaker::TR2LocalX(AliTrackReference *trackRef,
022044bf 1449 AliTPCParam *paramTPC) const {
c92725b7 1450
1451 Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
1452 Int_t index[4];
1453 paramTPC->Transform0to1(x,index);
1454 paramTPC->Transform1to2(x,index);
1455 return x[0];
1456}
1457////////////////////////////////////////////////////////////////////////
1458
1459
1460