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