]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PMD/AliPMDDigitizer.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Added explicit base class declaration
[u/mrichter/AliRoot.git] / PMD / AliPMDDigitizer.cxx
CommitLineData
b2b13ee2 1//-----------------------------------------------------//
2// //
8599803d 3// Source File : PMDDigitizer.cxx, Version 00 //
b2b13ee2 4// //
5// Date : September 20 2002 //
6// //
7//-----------------------------------------------------//
8
9#include <TBRIK.h>
10#include <TNode.h>
11#include <TTree.h>
12#include <TGeometry.h>
13#include <TObjArray.h>
14#include <TClonesArray.h>
15#include <TFile.h>
16#include <TNtuple.h>
17#include <TParticle.h>
18
19#include "AliRun.h"
20#include "AliPMD.h"
21#include "AliHit.h"
22#include "AliDetector.h"
23#include "AliRunLoader.h"
24#include "AliLoader.h"
25#include "AliConfig.h"
26#include "AliMagF.h"
27#include "AliRunDigitizer.h"
28#include "AliHeader.h"
29
30#include "AliPMDcell.h"
31#include "AliPMDsdigit.h"
32#include "AliPMDdigit.h"
33#include "AliPMDDigitizer.h"
34#include "AliPMDClustering.h"
b2b13ee2 35
36ClassImp(AliPMDDigitizer)
37//
38// Constructor
39//
40AliPMDDigitizer::AliPMDDigitizer()
41{
42 if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
43 fNsdigit = 0;
44 if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
45 fNdigit = 0;
46
8599803d 47 for (Int_t i = 0; i < fTotUM; i++)
b2b13ee2 48 {
8599803d 49 for (Int_t j = 0; j < fRow; j++)
b2b13ee2 50 {
8599803d 51 for (Int_t k = 0; k < fCol; k++)
b2b13ee2 52 {
53 fCPV[i][j][k] = 0.;
54 fPMD[i][j][k] = 0.;
55 }
56 }
57 }
58
59 if (!fCell) fCell = new TObjArray();
60
61 fZPos = 361.5; // in units of cm, This is the default position of PMD
b2b13ee2 62}
63AliPMDDigitizer::~AliPMDDigitizer()
64{
65 delete fSDigits;
66 delete fDigits;
67 delete fCell;
68}
69//
70// Member functions
71//
72void AliPMDDigitizer::OpengAliceFile(Char_t *file, Option_t *option)
73{
74
75 fRunLoader = AliRunLoader::Open(file,AliConfig::fgkDefaultEventFolderName,
76 "UPDATE");
77
78 if (!fRunLoader)
79 {
80 Error("Open","Can not open session for file %s.",file);
81 }
82
83 fRunLoader->LoadgAlice();
84 fRunLoader->LoadHeader();
85 fRunLoader->LoadKinematics();
86
87 gAlice = fRunLoader->GetAliRun();
88
89 if (gAlice)
90 {
91 printf("<AliPMDdigitizer::Open> ");
92 printf("AliRun object found on file.\n");
93 }
94 else
95 {
96 printf("<AliPMDdigitizer::Open> ");
97 printf("Could not find AliRun object.\n");
98 }
99
100 PMD = (AliPMD*)gAlice->GetDetector("PMD");
101 pmdloader = fRunLoader->GetLoader("PMDLoader");
102 if (pmdloader == 0x0)
103 {
104 cerr<<"Hits2Digits : Can not find PMD or PMDLoader\n";
105 }
106
107 const char *cHS = strstr(option,"HS");
108 const char *cHD = strstr(option,"HD");
109 const char *cSD = strstr(option,"SD");
110
111 if (cHS)
112 {
113 pmdloader->LoadHits("READ");
114 pmdloader->LoadSDigits("recreate");
115 }
116 else if (cHD)
117 {
118 pmdloader->LoadHits("READ");
119 pmdloader->LoadDigits("recreate");
120 }
121 else if (cSD)
122 {
123 pmdloader->LoadSDigits("READ");
124 pmdloader->LoadDigits("recreate");
125 }
126
127}
128void AliPMDDigitizer::Hits2SDigits(Int_t ievt)
129{
130 cout << " -------- Beginning of Hits2SDigits ----------- " << endl;
131
132 Int_t kPi0 = 111;
133 Int_t kGamma = 22;
134 Int_t npmd;
135 Int_t trackno;
b2b13ee2 136 Int_t smnumber;
137 Int_t trackpid;
138 Int_t mtrackno;
139 Int_t mtrackpid;
140
141 Float_t xPos, yPos, zPos;
8599803d 142 Int_t xpad = -1, ypad = -1;
b2b13ee2 143 Float_t edep;
144 Float_t vx = -999.0, vy = -999.0, vz = -999.0;
145
146
147 ResetSDigit();
148
149 printf("Event Number = %d \n",ievt);
150 Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
151 printf("Number of Particles = %d \n", nparticles);
152 fRunLoader->GetEvent(ievt);
153 Particles = gAlice->Particles();
154 // ------------------------------------------------------- //
155 // Pointer to specific detector hits.
156 // Get pointers to Alice detectors and Hits containers
157
158 treeH = pmdloader->TreeH();
159
160 Int_t ntracks = (Int_t) treeH->GetEntries();
161 printf("Number of Tracks in the TreeH = %d \n", ntracks);
162
163 treeS = pmdloader->TreeS();
164 if (treeS == 0x0)
165 {
166 pmdloader->MakeTree("S");
167 treeS = pmdloader->TreeS();
168 }
169 Int_t bufsize = 16000;
170 treeS->Branch("PMDSDigit", &fSDigits, bufsize);
171
172 if (PMD) PMDhits = PMD->Hits();
173
174 // Start loop on tracks in the hits containers
175
176
177 for (Int_t track=0; track<ntracks;track++)
178 {
179 gAlice->ResetHits();
180 treeH->GetEvent(track);
181
182 if (PMD)
183 {
184 npmd = PMDhits->GetEntriesFast();
185 for (int ipmd = 0; ipmd < npmd; ipmd++)
186 {
187 pmdHit = (AliPMDhit*) PMDhits->UncheckedAt(ipmd);
188 trackno = pmdHit->GetTrack();
189
190 // get kinematics of the particles
191
192 particle = gAlice->Particle(trackno);
193 trackpid = particle->GetPdgCode();
194
195 Int_t igatr = -999;
196 Int_t ichtr = -999;
8599803d 197 Int_t igapid = -999;
b2b13ee2 198 Int_t imo;
199 Int_t igen = 0;
200 Int_t id_mo = -999;
8599803d 201
b2b13ee2 202 TParticle* mparticle = particle;
8599803d 203 Int_t trackno_old=0, trackpid_old=0, status_old = 0;
204 if (mparticle->GetFirstMother() == -1)
205 {
206 trackno_old = trackno;
207 trackpid_old = trackpid;
208 status_old = -1;
209 }
210
211 Int_t igstatus = 0;
b2b13ee2 212 while((imo = mparticle->GetFirstMother()) >= 0)
213 {
214 igen++;
215 mparticle = gAlice->Particle(imo);
216 id_mo = mparticle->GetPdgCode();
217
218 vx = mparticle->Vx();
219 vy = mparticle->Vy();
220 vz = mparticle->Vz();
221
222 //printf("==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, id_mo, vx, vy, vz);
8599803d 223 //fprintf(fpw1,"==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, id_mo, vx, vy, vz);
224 if ((id_mo == kGamma || id_mo == -11 || id_mo == 11) && vx == 0. && vy == 0. && vz == 0.)
b2b13ee2 225 {
226 igatr = imo;
8599803d 227 igapid = id_mo;
228 igstatus = 1;
229 }
230 if(igstatus == 0)
231 {
232 if (id_mo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
233 {
234 igatr = imo;
235 igapid = id_mo;
236 }
b2b13ee2 237 }
238 ichtr = imo;
239 }
240
241 if (id_mo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
242 {
243 mtrackno = igatr;
8599803d 244 mtrackpid = igapid;
b2b13ee2 245 }
246 else
247 {
248 mtrackno = ichtr;
249 mtrackpid = id_mo;
250 }
8599803d 251 if (status_old == -1)
252 {
253 mtrackno = trackno_old;
254 mtrackpid = trackpid_old;
255 }
256
b2b13ee2 257 xPos = pmdHit->X();
258 yPos = pmdHit->Y();
259 zPos = pmdHit->Z();
b2b13ee2 260 edep = pmdHit->fEnergy;
8599803d 261 Int_t Vol1 = pmdHit->fVolume[1]; // Column
262 Int_t Vol2 = pmdHit->fVolume[2]; // Row
263 Int_t Vol3 = pmdHit->fVolume[3]; // UnitModule
264 Int_t Vol6 = pmdHit->fVolume[6]; // SuperModule
265 // -----------------------------------------//
266 // For Super Module 1 & 2 //
267 // nrow = 96, ncol = 48 //
268 // For Super Module 3 & 4 //
269 // nrow = 48, ncol = 96 //
270 // -----------------------------------------//
b2b13ee2 271
8599803d 272 smnumber = (Vol6-1)*6 + Vol3;
273
274 if (Vol6 == 1 || Vol6 == 2)
275 {
276 xpad = Vol1;
277 ypad = Vol2;
278 }
279 else if (Vol6 == 3 || Vol6 == 4)
b2b13ee2 280 {
8599803d 281 xpad = Vol2;
282 ypad = Vol1;
b2b13ee2 283 }
b2b13ee2 284
75bb524e 285 //cout << "zpos = " << zPos << " edep = " << edep << endl;
b2b13ee2 286
8599803d 287 Float_t zposition = TMath::Abs(zPos);
288 if (zposition < fZPos)
b2b13ee2 289 {
290 // CPV
291 fDetNo = 1;
292 }
8599803d 293 else if (zposition > fZPos)
b2b13ee2 294 {
295 // PMD
296 fDetNo = 0;
297 }
8599803d 298 Int_t smn = smnumber - 1;
299 Int_t ixx = xpad - 1;
300 Int_t iyy = ypad - 1;
b2b13ee2 301 if (fDetNo == 0)
302 {
8599803d 303 fPMD[smn][ixx][iyy] += edep;
304 fPMDCounter[smn][ixx][iyy]++;
b2b13ee2 305
306 pmdcell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
307
308 fCell->Add(pmdcell);
309 }
310 else if(fDetNo == 1)
311 {
8599803d 312 fCPV[smn][ixx][iyy] += edep;
313 fCPVTrackNo[smn][ixx][iyy] = mtrackno;
b2b13ee2 314 }
315 }
316 }
317 } // Track Loop ended
8599803d 318
b2b13ee2 319 TrackAssignment2Cell();
b2b13ee2 320 ResetCell();
321
322 Float_t deltaE = 0.;
323 Int_t detno = 0;
b2b13ee2 324 Int_t trno = -1;
325 Int_t cellno = 0;
326
327 for (Int_t idet = 0; idet < 2; idet++)
328 {
8599803d 329 for (Int_t ism = 0; ism < fTotUM; ism++)
b2b13ee2 330 {
8599803d 331 for (Int_t jrow = 0; jrow < fRow; jrow++)
b2b13ee2 332 {
8599803d 333 for (Int_t kcol = 0; kcol < fCol; kcol++)
b2b13ee2 334 {
8599803d 335 cellno = jrow*fCol + kcol;
b2b13ee2 336 if (idet == 0)
337 {
338 deltaE = fPMD[ism][jrow][kcol];
339 trno = fPMDTrackNo[ism][jrow][kcol];
340 detno = 0;
341 }
342 else if (idet == 1)
343 {
344 deltaE = fCPV[ism][jrow][kcol];
345 trno = fCPVTrackNo[ism][jrow][kcol];
346 detno = 1;
347 }
348 if (deltaE > 0.)
349 {
b2b13ee2 350 AddSDigit(trno,detno,ism,cellno,deltaE);
351 }
352 }
353 }
354 treeS->Fill();
355 ResetSDigit();
356 }
357 }
358 pmdloader->WriteSDigits("OVERWRITE");
b2b13ee2 359 ResetCellADC();
360
361 // cout << " -------- End of Hits2SDigit ----------- " << endl;
362}
363
364void AliPMDDigitizer::Hits2Digits(Int_t ievt)
365{
366 Int_t kPi0 = 111;
367 Int_t kGamma = 22;
368 Int_t npmd;
369 Int_t trackno;
b2b13ee2 370 Int_t smnumber;
371 Int_t trackpid;
372 Int_t mtrackno;
373 Int_t mtrackpid;
374
375 Float_t xPos, yPos, zPos;
8599803d 376 Int_t xpad = -1, ypad = -1;
b2b13ee2 377 Float_t edep;
378 Float_t vx = -999.0, vy = -999.0, vz = -999.0;
379
380
381 ResetDigit();
382
383 printf("Event Number = %d \n",ievt);
384
385 Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
386 printf("Number of Particles = %d \n", nparticles);
387 fRunLoader->GetEvent(ievt);
388 Particles = gAlice->Particles();
389 // ------------------------------------------------------- //
390 // Pointer to specific detector hits.
391 // Get pointers to Alice detectors and Hits containers
392
393 PMD = (AliPMD*)gAlice->GetDetector("PMD");
394 pmdloader = fRunLoader->GetLoader("PMDLoader");
395
396 if (pmdloader == 0x0)
397 {
398 cerr<<"Hits2Digits method : Can not find PMD or PMDLoader\n";
399 }
400 treeH = pmdloader->TreeH();
401 Int_t ntracks = (Int_t) treeH->GetEntries();
402 printf("Number of Tracks in the TreeH = %d \n", ntracks);
403 pmdloader->LoadDigits("recreate");
404 treeD = pmdloader->TreeD();
405 if (treeD == 0x0)
406 {
407 pmdloader->MakeTree("D");
408 treeD = pmdloader->TreeD();
409 }
410 Int_t bufsize = 16000;
411 treeD->Branch("PMDDigit", &fDigits, bufsize);
412
413 if (PMD) PMDhits = PMD->Hits();
414
415 // Start loop on tracks in the hits containers
416
417 for (Int_t track=0; track<ntracks;track++)
418 {
419 gAlice->ResetHits();
420 treeH->GetEvent(track);
421
422 if (PMD)
423 {
424 npmd = PMDhits->GetEntriesFast();
425 for (int ipmd = 0; ipmd < npmd; ipmd++)
426 {
427 pmdHit = (AliPMDhit*) PMDhits->UncheckedAt(ipmd);
428 trackno = pmdHit->GetTrack();
429
430 // get kinematics of the particles
431
432 particle = gAlice->Particle(trackno);
433 trackpid = particle->GetPdgCode();
434
435 Int_t igatr = -999;
436 Int_t ichtr = -999;
8599803d 437 Int_t igapid = -999;
b2b13ee2 438 Int_t imo;
439 Int_t igen = 0;
440 Int_t id_mo = -999;
8599803d 441
b2b13ee2 442 TParticle* mparticle = particle;
8599803d 443 Int_t trackno_old=0, trackpid_old=0, status_old = 0;
444 if (mparticle->GetFirstMother() == -1)
445 {
446 trackno_old = trackno;
447 trackpid_old = trackpid;
448 status_old = -1;
449 }
450
451 Int_t igstatus = 0;
b2b13ee2 452 while((imo = mparticle->GetFirstMother()) >= 0)
453 {
454 igen++;
455 mparticle = gAlice->Particle(imo);
456 id_mo = mparticle->GetPdgCode();
457
458 vx = mparticle->Vx();
459 vy = mparticle->Vy();
460 vz = mparticle->Vz();
461
462 //printf("==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, id_mo, vx, vy, vz);
8599803d 463 //fprintf(fpw1,"==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, id_mo, vx, vy, vz);
464 if ((id_mo == kGamma || id_mo == -11 || id_mo == 11) && vx == 0. && vy == 0. && vz == 0.)
b2b13ee2 465 {
466 igatr = imo;
8599803d 467 igapid = id_mo;
468 igstatus = 1;
469 }
470 if(igstatus == 0)
471 {
472 if (id_mo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
473 {
474 igatr = imo;
475 igapid = id_mo;
476 }
b2b13ee2 477 }
478 ichtr = imo;
479 }
480
481 if (id_mo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
482 {
483 mtrackno = igatr;
8599803d 484 mtrackpid = igapid;
b2b13ee2 485 }
486 else
487 {
488 mtrackno = ichtr;
489 mtrackpid = id_mo;
490 }
8599803d 491 if (status_old == -1)
492 {
493 mtrackno = trackno_old;
494 mtrackpid = trackpid_old;
495 }
b2b13ee2 496
497 xPos = pmdHit->X();
498 yPos = pmdHit->Y();
499 zPos = pmdHit->Z();
b2b13ee2 500 edep = pmdHit->fEnergy;
8599803d 501 Int_t Vol1 = pmdHit->fVolume[1]; // Column
502 Int_t Vol2 = pmdHit->fVolume[2]; // Row
503 Int_t Vol3 = pmdHit->fVolume[3]; // UnitModule
504 Int_t Vol6 = pmdHit->fVolume[6]; // SuperModule
505
506 // -----------------------------------------//
507 // For Super Module 1 & 2 //
508 // nrow = 96, ncol = 48 //
509 // For Super Module 3 & 4 //
510 // nrow = 48, ncol = 96 //
511 // -----------------------------------------//
b2b13ee2 512
8599803d 513 smnumber = (Vol6-1)*6 + Vol3;
514
515 if (Vol6 == 1 || Vol6 == 2)
b2b13ee2 516 {
8599803d 517 xpad = Vol1;
518 ypad = Vol2;
519 }
520 else if (Vol6 == 3 || Vol6 == 4)
521 {
522 xpad = Vol2;
523 ypad = Vol1;
b2b13ee2 524 }
b2b13ee2 525
526 //cout << "-zpos = " << -zPos << endl;
8599803d 527
528 Float_t zposition = TMath::Abs(zPos);
529
530 if (zposition < fZPos)
b2b13ee2 531 {
532 // CPV
533 fDetNo = 1;
534 }
8599803d 535 else if (zposition > fZPos)
b2b13ee2 536 {
537 // PMD
538 fDetNo = 0;
539 }
8599803d 540
541 Int_t smn = smnumber - 1;
542 Int_t ixx = xpad - 1;
543 Int_t iyy = ypad - 1;
544 if (fDetNo == 0)
b2b13ee2 545 {
8599803d 546 fPMD[smn][ixx][iyy] += edep;
547 fPMDCounter[smn][ixx][iyy]++;
548
549 pmdcell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
550
551 fCell->Add(pmdcell);
b2b13ee2 552 }
8599803d 553 else if(fDetNo == 1)
b2b13ee2 554 {
8599803d 555 fCPV[smn][ixx][iyy] += edep;
556 fCPVTrackNo[smn][ixx][iyy] = mtrackno;
b2b13ee2 557 }
558 }
559 }
560 } // Track Loop ended
561
562 TrackAssignment2Cell();
563 ResetCell();
564
565 Float_t deltaE = 0.;
566 Int_t detno = 0;
b2b13ee2 567 Int_t trno = 1;
568 Int_t cellno;
569
570 for (Int_t idet = 0; idet < 2; idet++)
571 {
8599803d 572 for (Int_t ism = 0; ism < fTotUM; ism++)
b2b13ee2 573 {
8599803d 574 for (Int_t jrow = 0; jrow < fRow; jrow++)
b2b13ee2 575 {
8599803d 576 for (Int_t kcol = 0; kcol < fCol; kcol++)
b2b13ee2 577 {
8599803d 578 cellno = jrow*fCol + kcol;
b2b13ee2 579 if (idet == 0)
580 {
581 deltaE = fPMD[ism][jrow][kcol];
8599803d 582 trno = fPMDTrackNo[ism][jrow][kcol];
b2b13ee2 583 detno = 0;
584 }
585 else if (idet == 1)
586 {
587 deltaE = fCPV[ism][jrow][kcol];
8599803d 588 trno = fCPVTrackNo[ism][jrow][kcol];
b2b13ee2 589 detno = 1;
590 }
591 if (deltaE > 0.)
592 {
b2b13ee2 593 AddDigit(trno,detno,ism,cellno,deltaE);
594 }
595 } // column loop
596 } // row loop
597 } // supermodule loop
598 treeD->Fill();
599 ResetDigit();
600 } // detector loop
601
602 pmdloader->WriteDigits("OVERWRITE");
b2b13ee2 603 ResetCellADC();
604
605 // cout << " -------- End of Hits2Digit ----------- " << endl;
606}
607
608
609void AliPMDDigitizer::SDigits2Digits(Int_t ievt)
610{
611 // cout << " -------- Beginning of SDigits2Digit ----------- " << endl;
612 fRunLoader->GetEvent(ievt);
613
614 treeS = pmdloader->TreeS();
615 AliPMDsdigit *pmdsdigit;
616 TBranch *branch = treeS->GetBranch("PMDSDigit");
617 branch->SetAddress(&fSDigits);
618
619 treeD = pmdloader->TreeD();
620 if (treeD == 0x0)
621 {
622 pmdloader->MakeTree("D");
623 treeD = pmdloader->TreeD();
624 }
625 Int_t bufsize = 16000;
626 treeD->Branch("PMDDigit", &fDigits, bufsize);
627
628 Int_t trno, det, smn;
629 Int_t cellno;
630 Float_t edep, adc;
631
632 Int_t nmodules = (Int_t) treeS->GetEntries();
633
634 for (Int_t imodule = 0; imodule < nmodules; imodule++)
635 {
636 treeS->GetEntry(imodule);
637 Int_t nentries = fSDigits->GetLast();
638 //cout << " nentries = " << nentries << endl;
639 for (Int_t ient = 0; ient < nentries+1; ient++)
640 {
641 pmdsdigit = (AliPMDsdigit*)fSDigits->UncheckedAt(ient);
642 trno = pmdsdigit->GetTrackNumber();
643 det = pmdsdigit->GetDetector();
644 smn = pmdsdigit->GetSMNumber();
645 cellno = pmdsdigit->GetCellNumber();
646 edep = pmdsdigit->GetCellEdep();
647
648 MeV2ADC(edep,adc);
b2b13ee2 649 AddDigit(trno,det,smn,cellno,adc);
650 }
651 treeD->Fill();
652 ResetDigit();
653 }
654 pmdloader->WriteDigits("OVERWRITE");
655 // cout << " -------- End of SDigits2Digit ----------- " << endl;
656}
657
658void AliPMDDigitizer::TrackAssignment2Cell()
659{
b2b13ee2 660 //
75bb524e 661 // This block assigns the cell id when there are
b2b13ee2 662 // multiple tracks in a cell according to the
663 // energy deposition
664 //
75bb524e 665 Bool_t jsort = false;
8599803d 666
b2b13ee2 667 Int_t i, j, k;
668
669 Float_t *frac_edp;
670 Float_t *tr_edp;
75bb524e 671 Int_t *status1;
672 Int_t *status2;
673 Int_t *trnarray;
b2b13ee2 674 Int_t ****PMDTrack;
675 Float_t ****PMDEdep;
676
8599803d 677 PMDTrack = new Int_t ***[fTotUM];
678 PMDEdep = new Float_t ***[fTotUM];
679 for (i=0; i<fTotUM; i++)
b2b13ee2 680 {
8599803d 681 PMDTrack[i] = new Int_t **[fRow];
682 PMDEdep[i] = new Float_t **[fRow];
b2b13ee2 683 }
684
8599803d 685 for (i = 0; i < fTotUM; i++)
b2b13ee2 686 {
8599803d 687 for (j = 0; j < fRow; j++)
b2b13ee2 688 {
8599803d 689 PMDTrack[i][j] = new Int_t *[fCol];
690 PMDEdep[i][j] = new Float_t *[fCol];
b2b13ee2 691 }
692 }
693
8599803d 694 for (i = 0; i < fTotUM; i++)
b2b13ee2 695 {
8599803d 696 for (j = 0; j < fRow; j++)
b2b13ee2 697 {
8599803d 698 for (k = 0; k < fCol; k++)
b2b13ee2 699 {
700 Int_t nn = fPMDCounter[i][j][k];
701 if(nn > 0)
702 {
703 PMDTrack[i][j][k] = new Int_t[nn];
704 PMDEdep[i][j][k] = new Float_t[nn];
705 }
706 else
707 {
708 nn = 1;
709 PMDTrack[i][j][k] = new Int_t[nn];
710 PMDEdep[i][j][k] = new Float_t[nn];
711 }
712 fPMDCounter[i][j][k] = 0;
713 }
714 }
715 }
716
717
718 Int_t nentries = fCell->GetEntries();
719
720 Int_t mtrackno, ism, ixp, iyp;
721 Float_t edep;
722
723 for (i = 0; i < nentries; i++)
724 {
725 pmdcell = (AliPMDcell*)fCell->UncheckedAt(i);
726
727 mtrackno = pmdcell->GetTrackNumber();
728 ism = pmdcell->GetSMNumber();
729 ixp = pmdcell->GetX();
730 iyp = pmdcell->GetY();
731 edep = pmdcell->GetEdep();
b2b13ee2 732 Int_t nn = fPMDCounter[ism][ixp][iyp];
b2b13ee2 733 // cout << " nn = " << nn << endl;
b2b13ee2 734 PMDTrack[ism][ixp][iyp][nn] = (Int_t) mtrackno;
735 PMDEdep[ism][ixp][iyp][nn] = edep;
736 fPMDCounter[ism][ixp][iyp]++;
737 }
738
739 Int_t iz, il;
740 Int_t im, ix, iy;
741 Int_t nn;
742
8599803d 743 for (im=0; im<fTotUM; im++)
b2b13ee2 744 {
8599803d 745 for (ix=0; ix<fRow; ix++)
b2b13ee2 746 {
8599803d 747 for (iy=0; iy<fCol; iy++)
b2b13ee2 748 {
749 nn = fPMDCounter[im][ix][iy];
750 if (nn > 1)
751 {
752 // This block handles if a cell is fired
753 // many times by many tracks
75bb524e 754 status1 = new Int_t[nn];
755 status2 = new Int_t[nn];
756 trnarray = new Int_t[nn];
b2b13ee2 757 for (iz = 0; iz < nn; iz++)
758 {
75bb524e 759 status1[iz] = PMDTrack[im][ix][iy][iz];
b2b13ee2 760 }
75bb524e 761 TMath::Sort(nn,status1,status2,jsort);
b2b13ee2 762 Int_t track_old = -99999;
763 Int_t track, tr_count = 0;
764 for (iz = 0; iz < nn; iz++)
765 {
75bb524e 766 track = status1[status2[iz]];
b2b13ee2 767 if (track_old != track)
768 {
75bb524e 769 trnarray[tr_count] = track;
b2b13ee2 770 tr_count++;
b2b13ee2 771 }
772 track_old = track;
773 }
75bb524e 774 delete status1;
775 delete status2;
b2b13ee2 776 Float_t tot_edp = 0.;
777 tr_edp = new Float_t[tr_count];
778 frac_edp = new Float_t[tr_count];
779 for (il = 0; il < tr_count; il++)
780 {
781 tr_edp[il] = 0.;
75bb524e 782 track = trnarray[il];
b2b13ee2 783 for (iz = 0; iz < nn; iz++)
784 {
785 if (track == PMDTrack[im][ix][iy][iz])
786 {
787 tr_edp[il] += PMDEdep[im][ix][iy][iz];
788 }
789 }
790 tot_edp += tr_edp[il];
791 }
b2b13ee2 792 Int_t il_old = 0;
793 Float_t frac_old = 0.;
794
795 for (il = 0; il < tr_count; il++)
796 {
797 frac_edp[il] = tr_edp[il]/tot_edp;
798 if (frac_old < frac_edp[il])
799 {
800 frac_old = frac_edp[il];
801 il_old = il;
802 }
803 }
75bb524e 804 fPMDTrackNo[im][ix][iy] = trnarray[il_old];
b2b13ee2 805 delete frac_edp;
806 delete tr_edp;
75bb524e 807 delete trnarray;
b2b13ee2 808 }
809 else if (nn == 1)
810 {
811 // This only handles if a cell is fired
812 // by only one track
813
814 fPMDTrackNo[im][ix][iy] = PMDTrack[im][ix][iy][0];
815
816 }
817 else if (nn ==0)
818 {
819 // This is if no cell is fired
820 fPMDTrackNo[im][ix][iy] = -999;
821 }
822 } // end of iy
823 } // end of ix
824 } // end of im
825
826 // Delete all the pointers
827
8599803d 828 for (i = 0; i < fTotUM; i++)
b2b13ee2 829 {
8599803d 830 for (j = 0; j < fRow; j++)
b2b13ee2 831 {
8599803d 832 for (k = 0; k < fCol; k++)
b2b13ee2 833 {
834 delete [] PMDTrack[i][j][k];
835 delete [] PMDEdep[i][j][k];
836 }
837 }
838 }
839
8599803d 840 for (i = 0; i < fTotUM; i++)
b2b13ee2 841 {
8599803d 842 for (j = 0; j < fRow; j++)
b2b13ee2 843 {
844 delete [] PMDTrack[i][j];
845 delete [] PMDEdep[i][j];
846 }
847 }
848
8599803d 849 for (i = 0; i < fTotUM; i++)
b2b13ee2 850 {
851 delete [] PMDTrack[i];
852 delete [] PMDEdep[i];
853 }
854 delete PMDTrack;
855 delete PMDEdep;
856 //
857 // End of the cell id assignment
858 //
859}
860
861
862void AliPMDDigitizer::MeV2ADC(Float_t mev, Float_t & adc)
863{
864 // To be done
865
866 adc = mev*1.;
867}
868void AliPMDDigitizer::AddSDigit(Int_t trnumber, Int_t det, Int_t smnumber,
869 Int_t cellnumber, Float_t adc)
870{
871 TClonesArray &lsdigits = *fSDigits;
872 AliPMDsdigit *newcell;
873 newcell = new AliPMDsdigit(trnumber,det,smnumber,cellnumber,adc);
874 new(lsdigits[fNsdigit++]) AliPMDsdigit(newcell);
875 delete newcell;
876}
877
878void AliPMDDigitizer::AddDigit(Int_t trnumber, Int_t det, Int_t smnumber,
879 Int_t cellnumber, Float_t adc)
880{
881 TClonesArray &ldigits = *fDigits;
882 AliPMDdigit *newcell;
883 newcell = new AliPMDdigit(trnumber,det,smnumber,cellnumber,adc);
884 new(ldigits[fNdigit++]) AliPMDdigit(newcell);
885 delete newcell;
886}
887
888Int_t AliPMDDigitizer::Convert2RealSMNumber(Int_t smnumber1)
889{
890 Int_t smnumber = -999;
891
892 if (smnumber1==1) smnumber = 1;
893 if (smnumber1==2) smnumber = 10;
894 if (smnumber1==3) smnumber = 19;
895 if (smnumber1==4) smnumber = 1;
896 if (smnumber1==5) smnumber = 10;
897 if (smnumber1==6) smnumber = 19;
898 if (smnumber1==7) smnumber = 2;
899 if (smnumber1==8) smnumber = 3;
900 if (smnumber1==9) smnumber = 4;
901 if (smnumber1==10) smnumber = 5;
902 if (smnumber1==11) smnumber = 6;
903 if (smnumber1==12) smnumber = 7;
904 if (smnumber1==13) smnumber = 8;
905 if (smnumber1==14) smnumber = 9;
906 if (smnumber1==15) smnumber = 11;
907 if (smnumber1==16) smnumber = 12;
908 if (smnumber1==17) smnumber = 13;
909 if (smnumber1==18) smnumber = 14;
910 if (smnumber1==19) smnumber = 15;
911 if (smnumber1==20) smnumber = 16;
912 if (smnumber1==21) smnumber = 17;
913 if (smnumber1==22) smnumber = 18;
914 if (smnumber1==23) smnumber = 20;
915 if (smnumber1==24) smnumber = 21;
916 if (smnumber1==25) smnumber = 22;
917 if (smnumber1==26) smnumber = 23;
918 if (smnumber1==27) smnumber = 24;
919 if (smnumber1==28) smnumber = 25;
920 if (smnumber1==29) smnumber = 26;
921 if (smnumber1==30) smnumber = 27;
922
923 return smnumber;
924}
925void AliPMDDigitizer::SetZPosition(Float_t zpos)
926{
927 fZPos = zpos;
928}
929Float_t AliPMDDigitizer::GetZPosition() const
930{
931 return fZPos;
932}
933
934void AliPMDDigitizer::ResetCell()
935{
936 fCell->Clear();
8599803d 937 for (Int_t i = 0; i < fTotUM; i++)
b2b13ee2 938 {
8599803d 939 for (Int_t j = 0; j < fRow; j++)
b2b13ee2 940 {
8599803d 941 for (Int_t k = 0; k < fCol; k++)
b2b13ee2 942 {
943 fPMDCounter[i][j][k] = 0;
944 }
945 }
946 }
947}
948void AliPMDDigitizer::ResetSDigit()
949{
950 fNsdigit = 0;
951 if (fSDigits) fSDigits->Clear();
952}
953void AliPMDDigitizer::ResetDigit()
954{
955 fNdigit = 0;
956 if (fDigits) fDigits->Clear();
957}
958
959void AliPMDDigitizer::ResetCellADC()
960{
8599803d 961 for (Int_t i = 0; i < fTotUM; i++)
b2b13ee2 962 {
8599803d 963 for (Int_t j = 0; j < fRow; j++)
b2b13ee2 964 {
8599803d 965 for (Int_t k = 0; k < fCol; k++)
b2b13ee2 966 {
967 fCPV[i][j][k] = 0.;
968 fPMD[i][j][k] = 0.;
969 }
970 }
971 }
972}
973
974void AliPMDDigitizer::UnLoad(Option_t *option)
975{
976 const char *cS = strstr(option,"S");
977 const char *cD = strstr(option,"D");
978
979 fRunLoader->UnloadgAlice();
980 fRunLoader->UnloadHeader();
981 fRunLoader->UnloadKinematics();
982
983 if (cS)
984 {
985 pmdloader->UnloadHits();
986 }
987 if (cD)
988 {
989 pmdloader->UnloadHits();
990 pmdloader->UnloadSDigits();
991 }
992}