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