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