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