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