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