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