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