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