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