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