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