Bug fixes
[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 "AliPMDhit.h"
37 #include "AliHit.h"
38 #include "AliDetector.h"
39 #include "AliRunLoader.h"
40 #include "AliLoader.h"
41 #include "AliConfig.h"
42 #include "AliMagF.h"
43 #include "AliRunDigitizer.h"
44 #include "AliDigitizer.h"
45 #include "AliHeader.h"
46
47 #include "AliPMDcell.h"
48 #include "AliPMDsdigit.h"
49 #include "AliPMDdigit.h"
50 #include "AliPMDDigitizer.h"
51 #include "AliMC.h"
52
53 ClassImp(AliPMDDigitizer)
54
55 AliPMDDigitizer::AliPMDDigitizer() :
56   fRunLoader(0),
57   fPMDHit(0),
58   fPMD(0),
59   fPMDLoader(0),
60   fSDigits(0),
61   fDigits(0),
62   fCell(0),
63   fDebug(0),
64   fNsdigit(0),
65   fNdigit(0),
66   fDetNo(0),
67   fZPos(361.5)// in units of cm, This is the default position of PMD
68 {
69   // Default Constructor
70   //
71   for (Int_t i = 0; i < fgkTotUM; i++)
72     {
73       for (Int_t j = 0; j < fgkRow; j++)
74         {
75           for (Int_t k = 0; k < fgkCol; k++)
76             {
77               fCPV[i][j][k] = 0.; 
78               fPRE[i][j][k] = 0.; 
79               fPRECounter[i][j][k] =  0; 
80               fPRETrackNo[i][j][k] = -1; 
81               fCPVTrackNo[i][j][k] = -1; 
82             }
83         }
84     }
85
86 }
87 //____________________________________________________________________________
88 AliPMDDigitizer::AliPMDDigitizer(AliRunDigitizer* manager) 
89   :AliDigitizer(manager),
90   fRunLoader(0),
91   fPMDHit(0),
92   fPMD(0),
93   fPMDLoader(0),
94   fSDigits(new TClonesArray("AliPMDsdigit", 1000)),
95   fDigits(new TClonesArray("AliPMDdigit", 1000)),
96   fCell(0),
97   fDebug(0),
98   fNsdigit(0),
99   fNdigit(0),
100   fDetNo(0),
101   fZPos(361.5)// in units of cm, This is the default position of PMD
102 {
103   // ctor which should be used
104
105   for (Int_t i = 0; i < fgkTotUM; i++)
106     {
107       for (Int_t j = 0; j < fgkRow; j++)
108         {
109           for (Int_t k = 0; k < fgkCol; k++)
110             {
111               fCPV[i][j][k] = 0.; 
112               fPRE[i][j][k] = 0.; 
113               fPRECounter[i][j][k] =  0; 
114               fPRETrackNo[i][j][k] = -1; 
115               fCPVTrackNo[i][j][k] = -1; 
116             }
117         }
118     }
119 }
120 //____________________________________________________________________________
121 AliPMDDigitizer::~AliPMDDigitizer()
122 {
123   // Default Destructor
124   //
125   if (fSDigits) {
126     fSDigits->Delete();
127     delete fSDigits;
128     fSDigits=0;
129   }
130   if (fDigits) {
131     fDigits->Delete();
132     delete fDigits;
133     fDigits=0;
134   }
135   fCell.Delete();
136 }
137 //
138 // Member functions
139 //
140 //____________________________________________________________________________
141 void AliPMDDigitizer::OpengAliceFile(const char *file, Option_t *option)
142 {
143   // Loads galice.root file and corresponding header, kinematics
144   // hits and sdigits or digits depending on the option
145   //
146
147   TString evfoldname = AliConfig::GetDefaultEventFolderName();
148   fRunLoader = AliRunLoader::GetRunLoader(evfoldname);
149   if (!fRunLoader)
150     fRunLoader = AliRunLoader::Open(file,AliConfig::GetDefaultEventFolderName(),
151                                     "UPDATE");
152   
153   if (!fRunLoader)
154    {
155      Error("Open","Can not open session for file %s.",file);
156    }
157   
158   if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
159   if (!fRunLoader->TreeE()) fRunLoader->LoadHeader();
160   if (!fRunLoader->TreeK()) fRunLoader->LoadKinematics();
161
162   gAlice = fRunLoader->GetAliRun();
163   
164   if (fDebug) {
165     if (gAlice)
166       {
167         printf("<AliPMDdigitizer::Open> ");
168         printf("AliRun object found on file.\n");
169       }
170     else
171       {
172         printf("<AliPMDdigitizer::Open> ");
173         printf("Could not find AliRun object.\n");
174       }
175   }
176   
177   fPMD  = (AliPMD*)gAlice->GetDetector("PMD");
178   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
179   if (fPMDLoader == 0x0)
180     {
181       cerr<<"Hits2Digits : Can not find PMD or PMDLoader\n";
182     }
183
184   const char *cHS = strstr(option,"HS");
185   const char *cHD = strstr(option,"HD");
186   const char *cSD = strstr(option,"SD");
187
188   if (cHS)
189     {
190       fPMDLoader->LoadHits("READ");
191       fPMDLoader->LoadSDigits("recreate");
192     }
193   else if (cHD)
194     {
195       fPMDLoader->LoadHits("READ");
196       fPMDLoader->LoadDigits("recreate");
197     }
198   else if (cSD)
199     {
200       fPMDLoader->LoadSDigits("READ");
201       fPMDLoader->LoadDigits("recreate");
202     }
203 }
204 //____________________________________________________________________________
205 void AliPMDDigitizer::Hits2SDigits(Int_t ievt)
206 {
207   // This reads the PMD Hits tree and assigns the right track number
208   // to a cell and stores in the summable digits tree
209   //
210   // cout << " -------- Beginning of Hits2SDigits ----------- " << endl;
211
212   const Int_t kPi0 = 111;
213   const Int_t kGamma = 22;
214   Int_t npmd;
215   Int_t trackno;
216   Int_t smnumber;
217   Int_t trackpid;
218   Int_t mtrackno;
219   Int_t mtrackpid;
220
221   Float_t xPos, yPos, zPos;
222   Int_t xpad = -1, ypad = -1;
223   Float_t edep;
224   Float_t vx = -999.0, vy = -999.0, vz = -999.0;
225
226
227   if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
228   ResetSDigit();
229
230   if (fDebug) printf("Event Number =  %d \n",ievt); 
231   Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
232   if (fDebug) printf("Number of Particles = %d \n", nparticles);
233   fRunLoader->GetEvent(ievt);
234   // ------------------------------------------------------- //
235   // Pointer to specific detector hits.
236   // Get pointers to Alice detectors and Hits containers
237
238   TTree* treeH = fPMDLoader->TreeH();
239   
240   Int_t ntracks    = (Int_t) treeH->GetEntries();
241   if (fDebug) printf("Number of Tracks in the TreeH = %d \n", ntracks);
242
243   TTree* treeS = fPMDLoader->TreeS();
244   if (treeS == 0x0)
245     {
246       fPMDLoader->MakeTree("S");
247       treeS = fPMDLoader->TreeS();
248     }
249   Int_t bufsize = 16000;
250   treeS->Branch("PMDSDigit", &fSDigits, bufsize); 
251   
252   TClonesArray* hits = 0;
253   if (fPMD) hits = fPMD->Hits();
254
255   // Start loop on tracks in the hits containers
256
257   for (Int_t track=0; track<ntracks;track++) 
258     {
259       gAlice->ResetHits();
260       treeH->GetEvent(track);
261       if (fPMD) 
262         {
263           npmd = hits->GetEntriesFast();
264           for (int ipmd = 0; ipmd < npmd; ipmd++) 
265             {
266               fPMDHit = (AliPMDhit*) hits->UncheckedAt(ipmd);
267               trackno = fPMDHit->GetTrack();
268               //  get kinematics of the particles
269
270               TParticle* mparticle = gAlice->GetMCApp()->Particle(trackno);
271               trackpid  = mparticle->GetPdgCode();
272
273               Int_t igatr = -999;
274               Int_t ichtr = -999;
275               Int_t igapid = -999;
276               Int_t imo;
277               Int_t igen = 0;
278               Int_t idmo = -999;
279
280               Int_t tracknoOld=0, trackpidOld=0, statusOld = 0;
281               if (mparticle->GetFirstMother() == -1)
282                 {
283                   tracknoOld  = trackno;
284                   trackpidOld = trackpid;
285                   statusOld   = -1;
286                 }
287               Int_t igstatus = 0;
288               while((imo = mparticle->GetFirstMother()) >= 0)
289                 {
290                   igen++;
291
292                   mparticle =  gAlice->GetMCApp()->Particle(imo);
293                   idmo = mparticle->GetPdgCode();
294                   
295                   vx = mparticle->Vx();
296                   vy = mparticle->Vy();
297                   vz = mparticle->Vz();
298                 
299                   //printf("==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, idmo, vx, vy, vz);
300                   //fprintf(fpw1,"==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, idmo, vx, vy, vz);
301                   if ((idmo == kGamma || idmo == -11 || idmo == 11) && vx == 0. && vy == 0. && vz == 0.)
302                     {
303                       igatr = imo;
304                       igapid = idmo;
305                       igstatus = 1;
306                     }
307                   if(igstatus == 0)
308                     {
309                       if (idmo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
310                         {
311                           igatr = imo;
312                           igapid = idmo;
313                         }
314                     }
315                   ichtr = imo;
316                 }
317
318               if (idmo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
319                 {
320                   mtrackno = igatr;
321                   mtrackpid = igapid;
322                 }
323               else
324                 {
325                   mtrackno  = ichtr;
326                   mtrackpid = idmo;
327                 }
328               if (statusOld == -1)
329                 {
330                   mtrackno  = tracknoOld;
331                   mtrackpid = trackpidOld;
332                 }
333               xPos = fPMDHit->X();
334               yPos = fPMDHit->Y();
335               zPos = fPMDHit->Z();
336               
337               edep       = fPMDHit->GetEnergy();
338               Int_t vol1 = fPMDHit->GetVolume(1); // Column
339               Int_t vol2 = fPMDHit->GetVolume(2); // Row
340               Int_t vol3 = fPMDHit->GetVolume(3); // UnitModule
341               Int_t vol6 = fPMDHit->GetVolume(6); // SuperModule
342
343               // -----------------------------------------//
344               // For Super Module 1 & 2                   //
345               //  nrow = 96, ncol = 48                    //
346               // For Super Module 3 & 4                   //
347               //  nrow = 48, ncol = 96                    //
348               // -----------------------------------------//
349               
350               smnumber = (vol6-1)*6 + vol3;
351
352               if (vol6 == 1 || vol6 == 2)
353                 {
354                   xpad = vol1;
355                   ypad = vol2;
356                 }
357               else if (vol6 == 3 || vol6 == 4)
358                 {
359                   xpad = vol2;
360                   ypad = vol1;
361                 }
362
363               //cout << "zpos = " << zPos << " edep = " << edep << endl;
364
365               Float_t zposition = TMath::Abs(zPos);
366               if (zposition < fZPos)
367                 {
368                   // CPV
369                   fDetNo = 1;
370                 }
371               else if (zposition > fZPos)
372                 {
373                   // PMD
374                   fDetNo = 0;
375                 }
376               Int_t smn = smnumber - 1;
377               Int_t ixx = xpad     - 1;
378               Int_t iyy = ypad     - 1;
379               if (fDetNo == 0)
380                 {
381                   fPRE[smn][ixx][iyy] += edep;
382                   fPRECounter[smn][ixx][iyy]++;
383
384                   AliPMDcell* cell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
385
386                   fCell.Add(cell);
387                 }
388               else if(fDetNo == 1)
389                 {
390                   fCPV[smn][ixx][iyy] += edep;
391                   fCPVTrackNo[smn][ixx][iyy] = mtrackno;
392                 }
393             }
394         }
395     } // Track Loop ended
396
397   TrackAssignment2Cell();
398   ResetCell();
399
400   Float_t deltaE      = 0.;
401   Int_t   detno       = 0;
402   Int_t   trno        = -1;
403
404   for (Int_t idet = 0; idet < 2; idet++)
405     {
406       for (Int_t ism = 0; ism < fgkTotUM; ism++)
407         {
408           for (Int_t jrow = 0; jrow < fgkRow; jrow++)
409             {
410               for (Int_t kcol = 0; kcol < fgkCol; kcol++)
411                 {
412                   if (idet == 0)
413                     {
414                       deltaE = fPRE[ism][jrow][kcol];
415                       trno   = fPRETrackNo[ism][jrow][kcol];
416                       detno = 0;
417                     }
418                   else if (idet == 1)
419                     {
420                       deltaE = fCPV[ism][jrow][kcol];
421                       trno   = fCPVTrackNo[ism][jrow][kcol];
422                       detno = 1;
423                     }
424                   if (deltaE > 0.)
425                     {
426                       AddSDigit(trno,detno,ism,jrow,kcol,deltaE);
427                     }
428                 }
429             }
430           treeS->Fill();
431           ResetSDigit();
432         }
433     }
434   fPMDLoader->WriteSDigits("OVERWRITE");
435   ResetCellADC();
436
437   //  cout << " -------- End of Hits2SDigit ----------- " << endl;
438 }
439 //____________________________________________________________________________
440
441 void AliPMDDigitizer::Hits2Digits(Int_t ievt)
442 {
443   // This reads the PMD Hits tree and assigns the right track number
444   // to a cell and stores in the digits tree
445   //
446   const Int_t kPi0 = 111;
447   const Int_t kGamma = 22;
448   Int_t npmd;
449   Int_t trackno;
450   Int_t smnumber;
451   Int_t trackpid;
452   Int_t mtrackno;
453   Int_t mtrackpid;
454
455   Float_t xPos, yPos, zPos;
456   Int_t xpad = -1, ypad = -1;
457   Float_t edep;
458   Float_t vx = -999.0, vy = -999.0, vz = -999.0;
459
460   if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
461   ResetDigit();
462
463   if (fDebug) printf("Event Number =  %d \n",ievt); 
464
465   Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
466   if (fDebug) printf("Number of Particles = %d \n", nparticles);
467   fRunLoader->GetEvent(ievt);
468   // ------------------------------------------------------- //
469   // Pointer to specific detector hits.
470   // Get pointers to Alice detectors and Hits containers
471
472   fPMD  = (AliPMD*)gAlice->GetDetector("PMD");
473   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
474
475   if (fPMDLoader == 0x0)
476     {
477       cerr<<"Hits2Digits method : Can not find PMD or PMDLoader\n";
478     }
479   TTree* treeH = fPMDLoader->TreeH();
480   Int_t ntracks    = (Int_t) treeH->GetEntries();
481   if (fDebug) printf("Number of Tracks in the TreeH = %d \n", ntracks);
482   fPMDLoader->LoadDigits("recreate");
483   TTree* treeD = fPMDLoader->TreeD();
484   if (treeD == 0x0)
485     {
486       fPMDLoader->MakeTree("D");
487       treeD = fPMDLoader->TreeD();
488     }
489   Int_t bufsize = 16000;
490   treeD->Branch("PMDDigit", &fDigits, bufsize); 
491   
492   TClonesArray* hits = 0;
493   if (fPMD) hits = fPMD->Hits();
494
495   // Start loop on tracks in the hits containers
496
497   for (Int_t track=0; track<ntracks;track++) 
498     {
499       gAlice->ResetHits();
500       treeH->GetEvent(track);
501       
502       if (fPMD) 
503         {
504           npmd = hits->GetEntriesFast();
505           for (int ipmd = 0; ipmd < npmd; ipmd++) 
506             {
507               fPMDHit = (AliPMDhit*) hits->UncheckedAt(ipmd);
508               trackno = fPMDHit->GetTrack();
509               
510               //  get kinematics of the particles
511               
512               TParticle* mparticle = gAlice->GetMCApp()->Particle(trackno);
513               trackpid  = mparticle->GetPdgCode();
514
515               Int_t igatr = -999;
516               Int_t ichtr = -999;
517               Int_t igapid = -999;
518               Int_t imo;
519               Int_t igen = 0;
520               Int_t idmo = -999;
521
522               Int_t tracknoOld=0, trackpidOld=0, statusOld = 0;
523               if (mparticle->GetFirstMother() == -1)
524                 {
525                   tracknoOld  = trackno;
526                   trackpidOld = trackpid;
527                   statusOld   = -1;
528                 }
529
530               Int_t igstatus = 0;
531               while((imo = mparticle->GetFirstMother()) >= 0)
532                 {
533                   igen++;
534
535                   mparticle =  gAlice->GetMCApp()->Particle(imo);
536                   idmo = mparticle->GetPdgCode();
537                   
538                   vx = mparticle->Vx();
539                   vy = mparticle->Vy();
540                   vz = mparticle->Vz();
541                 
542                   //printf("==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, idmo, vx, vy, vz);
543                   //fprintf(fpw1,"==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, idmo, vx, vy, vz);
544                   if ((idmo == kGamma || idmo == -11 || idmo == 11) && vx == 0. && vy == 0. && vz == 0.)
545                     {
546                       igatr = imo;
547                       igapid = idmo;
548                       igstatus = 1;
549                     }
550                   if(igstatus == 0)
551                     {
552                       if (idmo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
553                         {
554                           igatr = imo;
555                           igapid = idmo;
556                         }
557                     }
558                   ichtr = imo;
559                 }
560
561               if (idmo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
562                 {
563                   mtrackno = igatr;
564                   mtrackpid = igapid;
565                 }
566               else
567                 {
568                   mtrackno  = ichtr;
569                   mtrackpid = idmo;
570                 }
571               if (statusOld == -1)
572                 {
573                   mtrackno  = tracknoOld;
574                   mtrackpid = trackpidOld;
575                 }
576               
577               xPos = fPMDHit->X();
578               yPos = fPMDHit->Y();
579               zPos = fPMDHit->Z();
580
581               edep       = fPMDHit->GetEnergy();
582               Int_t vol1 = fPMDHit->GetVolume(1); // Column
583               Int_t vol2 = fPMDHit->GetVolume(2); // Row
584               Int_t vol3 = fPMDHit->GetVolume(3); // UnitModule
585               Int_t vol6 = fPMDHit->GetVolume(6); // SuperModule
586
587               // -----------------------------------------//
588               // For Super Module 1 & 2                   //
589               //  nrow = 96, ncol = 48                    //
590               // For Super Module 3 & 4                   //
591               //  nrow = 48, ncol = 96                    //
592               // -----------------------------------------//
593               
594               smnumber = (vol6-1)*6 + vol3;
595
596               if (vol6 == 1 || vol6 == 2)
597                 {
598                   xpad = vol1;
599                   ypad = vol2;
600                 }
601               else if (vol6 == 3 || vol6 == 4)
602                 {
603                   xpad = vol2;
604                   ypad = vol1;
605                 }
606
607               //cout << "-zpos = " << -zPos << endl;
608
609               Float_t zposition = TMath::Abs(zPos);
610
611               if (zposition < fZPos)
612                 {
613                   // CPV
614                   fDetNo = 1;
615                 }
616               else if (zposition > fZPos)
617                 {
618                   // PMD
619                   fDetNo = 0;
620                 }
621
622               Int_t smn = smnumber - 1;
623               Int_t ixx = xpad     - 1;
624               Int_t iyy = ypad     - 1;
625               if (fDetNo == 0)
626                 {
627                   fPRE[smn][ixx][iyy] += edep;
628                   fPRECounter[smn][ixx][iyy]++;
629
630                   AliPMDcell* cell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
631
632                   fCell.Add(cell);
633                 }
634               else if(fDetNo == 1)
635                 {
636                   fCPV[smn][ixx][iyy] += edep;
637                   fCPVTrackNo[smn][ixx][iyy] = mtrackno;
638                 }
639             }
640         }
641     } // Track Loop ended
642
643   TrackAssignment2Cell();
644   ResetCell();
645
646   Float_t adc;
647   Float_t deltaE = 0.;
648   Int_t detno = 0;
649   Int_t trno = 1;
650
651   for (Int_t idet = 0; idet < 2; idet++)
652     {
653       for (Int_t ism = 0; ism < fgkTotUM; ism++)
654         {
655           for (Int_t jrow = 0; jrow < fgkRow; jrow++)
656             {
657               for (Int_t kcol = 0; kcol < fgkCol; kcol++)
658                 {
659                   if (idet == 0)
660                     {
661                       deltaE = fPRE[ism][jrow][kcol];
662                       trno   = fPRETrackNo[ism][jrow][kcol];
663                       detno = 0;
664                     }
665                   else if (idet == 1)
666                     {
667                       deltaE = fCPV[ism][jrow][kcol];
668                       trno   = fCPVTrackNo[ism][jrow][kcol];
669                       detno = 1;
670                     }
671                   if (deltaE > 0.)
672                     {
673                       MeV2ADC(deltaE,adc);
674                       AddDigit(trno,detno,ism,jrow,kcol,adc);
675                     }
676                 } // column loop
677             } // row    loop
678         } // supermodule loop
679       treeD->Fill();
680       ResetDigit();
681     } // detector loop
682
683   fPMDLoader->WriteDigits("OVERWRITE");
684   ResetCellADC();
685   
686   //  cout << " -------- End of Hits2Digit ----------- " << endl;
687 }
688 //____________________________________________________________________________
689
690
691 void AliPMDDigitizer::SDigits2Digits(Int_t ievt)
692 {
693   // This reads the PMD sdigits tree and converts energy deposition
694   // in a cell to ADC and stores in the digits tree
695   //
696   //  cout << " -------- Beginning of SDigits2Digit ----------- " << endl;
697   fRunLoader->GetEvent(ievt);
698
699   TTree* treeS = fPMDLoader->TreeS();
700   AliPMDsdigit  *pmdsdigit;
701   TBranch *branch = treeS->GetBranch("PMDSDigit");
702   if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
703   branch->SetAddress(&fSDigits);
704
705   TTree* treeD = fPMDLoader->TreeD();
706   if (treeD == 0x0)
707     {
708       fPMDLoader->MakeTree("D");
709       treeD = fPMDLoader->TreeD();
710     }
711   Int_t bufsize = 16000;
712   if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
713   treeD->Branch("PMDDigit", &fDigits, bufsize); 
714
715   Int_t   trno, det, smn;
716   Int_t   irow, icol;
717   Float_t edep, adc;
718
719   Int_t nmodules = (Int_t) treeS->GetEntries();
720
721   for (Int_t imodule = 0; imodule < nmodules; imodule++)
722     {
723       treeS->GetEntry(imodule); 
724       Int_t nentries = fSDigits->GetLast();
725       //cout << " nentries = " << nentries << endl;
726       for (Int_t ient = 0; ient < nentries+1; ient++)
727         {
728           pmdsdigit = (AliPMDsdigit*)fSDigits->UncheckedAt(ient);
729           trno   = pmdsdigit->GetTrackNumber();
730           det    = pmdsdigit->GetDetector();
731           smn    = pmdsdigit->GetSMNumber();
732           irow   = pmdsdigit->GetRow();
733           icol   = pmdsdigit->GetColumn();
734           edep   = pmdsdigit->GetCellEdep();
735
736           MeV2ADC(edep,adc);
737           AddDigit(trno,det,smn,irow,icol,adc);      
738         }
739       treeD->Fill();
740       ResetDigit();
741     }
742   fPMDLoader->WriteDigits("OVERWRITE");
743   //  cout << " -------- End of SDigits2Digit ----------- " << endl;
744 }
745 //____________________________________________________________________________
746 void AliPMDDigitizer::Exec(Option_t *option)
747 {
748   // Does the event merging and digitization
749
750   fDebug = 0;
751   const char *cdeb = strstr(option,"deb");
752   if(cdeb)
753     {
754       cout << "**************** PMD Exec *************** " << endl;
755       fDebug = 1;
756     }
757   
758   Int_t ninputs = fManager->GetNinputs();
759   if(fDebug)
760     {
761       cout << " Number of files = " << ninputs << endl;
762     }
763   ResetCellADC();
764
765   for (Int_t i = 0; i < ninputs; i++)
766     {
767       Int_t troffset = fManager->GetMask(i);
768       MergeSDigits(i, troffset);
769     }
770
771   fRunLoader = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
772   fPMD  = (AliPMD*)gAlice->GetDetector("PMD");
773   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
774   if (fPMDLoader == 0x0)
775     {
776       cerr<<"AliPMDDigitizer::Exec : Can not find PMD or PMDLoader\n";
777     }
778   fPMDLoader->LoadDigits("update");
779   TTree* treeD = fPMDLoader->TreeD();
780   if (treeD == 0x0)
781     {
782       fPMDLoader->MakeTree("D");
783       treeD = fPMDLoader->TreeD();
784     }
785   Int_t bufsize = 16000;
786   if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
787   treeD->Branch("PMDDigit", &fDigits, bufsize); 
788
789   Float_t adc;
790   Float_t deltaE = 0.;
791   Int_t detno = 0;
792   Int_t trno = 1;
793
794   for (Int_t idet = 0; idet < 2; idet++)
795     {
796       for (Int_t ism = 0; ism < fgkTotUM; ism++)
797         {
798           for (Int_t jrow = 0; jrow < fgkRow; jrow++)
799             {
800               for (Int_t kcol = 0; kcol < fgkCol; kcol++)
801                 {
802                   if (idet == 0)
803                     {
804                       deltaE = fPRE[ism][jrow][kcol];
805                       trno   = fPRETrackNo[ism][jrow][kcol];
806                       detno = 0;
807                     }
808                   else if (idet == 1)
809                     {
810                       deltaE = fCPV[ism][jrow][kcol];
811                       trno   = fCPVTrackNo[ism][jrow][kcol];
812                       detno = 1;
813                     }
814                   if (deltaE > 0.)
815                     {
816                       MeV2ADC(deltaE,adc);
817                       AddDigit(trno,detno,ism,jrow,kcol,adc);
818                     }
819                 } // column loop
820             } // row    loop
821           treeD->Fill();
822           ResetDigit();
823         } // supermodule loop
824     } // detector loop
825   fPMDLoader->WriteDigits("OVERWRITE");  
826   fPMDLoader->UnloadDigits();
827   ResetCellADC();
828 }
829 //____________________________________________________________________________
830
831 void AliPMDDigitizer::MergeSDigits(Int_t filenumber, Int_t troffset)
832 {
833   // merging sdigits
834   fRunLoader = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(filenumber));
835   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
836   fPMDLoader->LoadSDigits("read");
837   TTree* treeS = fPMDLoader->TreeS();
838   AliPMDsdigit  *pmdsdigit;
839   TBranch *branch = treeS->GetBranch("PMDSDigit");
840   if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
841   branch->SetAddress(&fSDigits);
842
843   Int_t   itrackno, idet, ism;
844   Int_t   ixp, iyp;
845   Float_t edep;
846   
847   Int_t nmodules = (Int_t) treeS->GetEntries();
848   if(fDebug)
849     {
850       cout << " nmodules = " << nmodules << endl;
851       cout << " tr offset = " << troffset << endl;
852     }
853   for (Int_t imodule = 0; imodule < nmodules; imodule++)
854     {
855       treeS->GetEntry(imodule); 
856       Int_t nentries = fSDigits->GetLast();
857       if(fDebug)
858         {
859           cout << " nentries = " << nentries << endl;
860         }
861       for (Int_t ient = 0; ient < nentries+1; ient++)
862         {
863           pmdsdigit = (AliPMDsdigit*)fSDigits->UncheckedAt(ient);
864           itrackno  = pmdsdigit->GetTrackNumber();
865           idet      = pmdsdigit->GetDetector();
866           ism       = pmdsdigit->GetSMNumber();
867           ixp       = pmdsdigit->GetRow();
868           iyp       = pmdsdigit->GetColumn();
869           edep      = pmdsdigit->GetCellEdep();
870
871           if (idet == 0)
872             {
873               if (fPRE[ism][ixp][iyp] < edep)
874                 {
875                   fPRETrackNo[ism][ixp][iyp] = troffset + itrackno;
876                 }
877               fPRE[ism][ixp][iyp] += edep;
878             }
879           else if (idet == 1)
880             {
881               if (fCPV[ism][ixp][iyp] < edep)
882                 {
883                   fCPVTrackNo[ism][ixp][iyp] = troffset + itrackno;
884                 }
885               fCPV[ism][ixp][iyp] += edep;
886             }
887         }
888     }
889
890 }
891 // ----------------------------------------------------------------------
892 void AliPMDDigitizer::TrackAssignment2Cell()
893 {
894   // 
895   // This block assigns the cell id when there are
896   // multiple tracks in a cell according to the
897   // energy deposition
898   //
899   Bool_t jsort = false;
900
901   Int_t i, j, k;
902
903   Float_t *fracEdp;
904   Float_t *trEdp;
905   Int_t *status1;
906   Int_t *status2;
907   Int_t *trnarray;
908   Int_t   ****pmdTrack;
909   Float_t ****pmdEdep;
910
911   pmdTrack = new Int_t ***[fgkTotUM];
912   pmdEdep  = new Float_t ***[fgkTotUM];
913   for (i=0; i<fgkTotUM; i++)
914     {
915       pmdTrack[i] = new Int_t **[fgkRow];
916       pmdEdep[i]  = new Float_t **[fgkRow];
917     }
918
919   for (i = 0; i < fgkTotUM; i++)
920     {
921       for (j = 0; j < fgkRow; j++)
922         {
923           pmdTrack[i][j] = new Int_t *[fgkCol];
924           pmdEdep[i][j]  = new Float_t *[fgkCol];
925         }
926     }
927   
928   for (i = 0; i < fgkTotUM; i++)
929     {
930       for (j = 0; j < fgkRow; j++)
931         {
932           for (k = 0; k < fgkCol; k++)
933             {
934               Int_t nn = fPRECounter[i][j][k];
935               if(nn > 0)
936                 {
937                   pmdTrack[i][j][k] = new Int_t[nn];
938                   pmdEdep[i][j][k] = new Float_t[nn];
939                 }
940               else
941                 {
942                   nn = 1;
943                   pmdTrack[i][j][k] = new Int_t[nn];
944                   pmdEdep[i][j][k] = new Float_t[nn];
945                 }                     
946               fPRECounter[i][j][k] = 0;
947             }
948         }
949     }
950
951
952   Int_t nentries = fCell.GetEntries();
953
954   Int_t   mtrackno, ism, ixp, iyp;
955   Float_t edep;
956
957   for (i = 0; i < nentries; i++)
958     {
959       AliPMDcell* cell = (AliPMDcell*)fCell.UncheckedAt(i);
960       
961       mtrackno = cell->GetTrackNumber();
962       ism      = cell->GetSMNumber();
963       ixp      = cell->GetX();
964       iyp      = cell->GetY();
965       edep     = cell->GetEdep();
966       Int_t nn = fPRECounter[ism][ixp][iyp];
967       //      cout << " nn = " << nn << endl;
968       pmdTrack[ism][ixp][iyp][nn] = (Int_t) mtrackno;
969       pmdEdep[ism][ixp][iyp][nn] = edep;
970       fPRECounter[ism][ixp][iyp]++;
971     }
972   
973   Int_t iz, il;
974   Int_t im, ix, iy;
975   Int_t nn;
976   
977   for (im=0; im<fgkTotUM; im++)
978     {
979       for (ix=0; ix<fgkRow; ix++)
980         {
981           for (iy=0; iy<fgkCol; iy++)
982             {
983               nn = fPRECounter[im][ix][iy];
984               if (nn > 1)
985                 {
986                   // This block handles if a cell is fired
987                   // many times by many tracks
988                   status1  = new Int_t[nn];
989                   status2  = new Int_t[nn];
990                   trnarray = new Int_t[nn];
991                   for (iz = 0; iz < nn; iz++)
992                     {
993                       status1[iz] = pmdTrack[im][ix][iy][iz];
994                     }
995                   TMath::Sort(nn,status1,status2,jsort);
996                   Int_t trackOld = -99999;
997                   Int_t track, trCount = 0;
998                   for (iz = 0; iz < nn; iz++)
999                     {
1000                       track = status1[status2[iz]];
1001                       if (trackOld != track)
1002                         {
1003                           trnarray[trCount] = track;
1004                           trCount++;
1005                         }                             
1006                       trackOld = track;
1007                     }
1008                   delete [] status1;
1009                   delete [] status2;
1010                   Float_t totEdp = 0.;
1011                   trEdp = new Float_t[trCount];
1012                   fracEdp = new Float_t[trCount];
1013                   for (il = 0; il < trCount; il++)
1014                     {
1015                       trEdp[il] = 0.;
1016                       track = trnarray[il];
1017                       for (iz = 0; iz < nn; iz++)
1018                         {
1019                           if (track == pmdTrack[im][ix][iy][iz])
1020                             {
1021                               trEdp[il] += pmdEdep[im][ix][iy][iz];
1022                             }
1023                         }
1024                       totEdp += trEdp[il];
1025                     }
1026                   Int_t ilOld = 0;
1027                   Float_t fracOld = 0.;
1028                   
1029                   for (il = 0; il < trCount; il++)
1030                     {
1031                       fracEdp[il] = trEdp[il]/totEdp;
1032                       if (fracOld < fracEdp[il])
1033                         {
1034                           fracOld = fracEdp[il];
1035                           ilOld = il;
1036                         }
1037                     }
1038                   fPRETrackNo[im][ix][iy] = trnarray[ilOld];
1039                   delete [] fracEdp;
1040                   delete [] trEdp;
1041                   delete [] trnarray;
1042                 }
1043               else if (nn == 1)
1044                 {
1045                   // This only handles if a cell is fired
1046                   // by only one track
1047                   
1048                   fPRETrackNo[im][ix][iy] = pmdTrack[im][ix][iy][0];
1049                   
1050                 }
1051               else if (nn ==0)
1052                 {
1053                   // This is if no cell is fired
1054                   fPRETrackNo[im][ix][iy] = -999;
1055                 }
1056             } // end of iy
1057         } // end of ix
1058     } // end of im
1059   
1060   // Delete all the pointers
1061   
1062   for (i = 0; i < fgkTotUM; i++)
1063     {
1064       for (j = 0; j < fgkRow; j++)
1065         {
1066           for (k = 0; k < fgkCol; k++)
1067             {
1068               delete [] pmdTrack[i][j][k];
1069               delete [] pmdEdep[i][j][k];
1070             }
1071         }
1072     }
1073   
1074   for (i = 0; i < fgkTotUM; i++)
1075     {
1076       for (j = 0; j < fgkRow; j++)
1077         {
1078           delete [] pmdTrack[i][j];
1079           delete [] pmdEdep[i][j];
1080         }
1081     }
1082   
1083   for (i = 0; i < fgkTotUM; i++)
1084     {
1085       delete [] pmdTrack[i];
1086       delete [] pmdEdep[i];
1087     }
1088   delete [] pmdTrack;
1089   delete [] pmdEdep;
1090   // 
1091   // End of the cell id assignment
1092   //
1093 }
1094 //____________________________________________________________________________
1095 void AliPMDDigitizer::MeV2ADC(Float_t mev, Float_t & adc) const
1096 {
1097   // This converts the simulated edep to ADC according to the
1098   // Test Beam Data
1099   // To be done
1100   //
1101
1102   //  adc = mev*1.;
1103
1104   
1105   // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% //
1106   // PS Test in September 2003
1107   // MeV - ADC conversion for 10bit ADC
1108
1109   const Float_t kConstant   = 7.181;
1110   const Float_t kErConstant = 0.6899;
1111   const Float_t kSlope      = 35.93;
1112   const Float_t kErSlope    = 0.306;
1113
1114   //gRandom->SetSeed();
1115
1116   Float_t cons   = gRandom->Gaus(kConstant,kErConstant);
1117   Float_t slop   = gRandom->Gaus(kSlope,kErSlope);
1118
1119   Float_t adc10bit = slop*mev*0.001 + cons;
1120
1121   // 12 bit adc
1122
1123   Int_t adc12bit = (Int_t) (4.0*adc10bit);
1124
1125   if(adc12bit < 3000)
1126     {
1127       adc = (Float_t) adc12bit;
1128     }
1129   else if (adc12bit >= 3000)
1130     {
1131       adc = 3000.0;
1132     }
1133
1134   // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% //
1135
1136 }
1137 //____________________________________________________________________________
1138 void AliPMDDigitizer::AddSDigit(Int_t trnumber, Int_t det, Int_t smnumber, 
1139   Int_t irow, Int_t icol, Float_t adc)
1140 {
1141   // Add SDigit
1142   //
1143   if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
1144   TClonesArray &lsdigits = *fSDigits;
1145   new(lsdigits[fNsdigit++])  AliPMDsdigit(trnumber,det,smnumber,irow,icol,adc);
1146 }
1147 //____________________________________________________________________________
1148
1149 void AliPMDDigitizer::AddDigit(Int_t trnumber, Int_t det, Int_t smnumber, 
1150   Int_t irow, Int_t icol, Float_t adc)
1151 {
1152   // Add Digit
1153   //
1154   if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
1155   TClonesArray &ldigits = *fDigits;
1156   new(ldigits[fNdigit++]) AliPMDdigit(trnumber,det,smnumber,irow,icol,adc);
1157 }
1158 //____________________________________________________________________________
1159
1160 void AliPMDDigitizer::SetZPosition(Float_t zpos)
1161 {
1162   fZPos = zpos;
1163 }
1164 //____________________________________________________________________________
1165 Float_t AliPMDDigitizer::GetZPosition() const
1166 {
1167   return fZPos;
1168 }
1169 //____________________________________________________________________________
1170
1171 void AliPMDDigitizer::ResetCell()
1172 {
1173   // clears the cell array and also the counter
1174   //  for each cell
1175   //
1176   fCell.Delete();
1177   for (Int_t i = 0; i < fgkTotUM; i++)
1178     {
1179       for (Int_t j = 0; j < fgkRow; j++)
1180         {
1181           for (Int_t k = 0; k < fgkCol; k++)
1182             {
1183               fPRECounter[i][j][k] = 0; 
1184             }
1185         }
1186     }
1187 }
1188 //____________________________________________________________________________
1189 void AliPMDDigitizer::ResetSDigit()
1190 {
1191   // Clears SDigits
1192   fNsdigit = 0;
1193   if (fSDigits) fSDigits->Delete();
1194 }
1195 //____________________________________________________________________________
1196 void AliPMDDigitizer::ResetDigit()
1197 {
1198   // Clears Digits
1199   fNdigit = 0;
1200   if (fDigits) fDigits->Delete();
1201 }
1202 //____________________________________________________________________________
1203
1204 void AliPMDDigitizer::ResetCellADC()
1205 {
1206   // Clears individual cells edep and track number
1207   for (Int_t i = 0; i < fgkTotUM; i++)
1208     {
1209       for (Int_t j = 0; j < fgkRow; j++)
1210         {
1211           for (Int_t k = 0; k < fgkCol; k++)
1212             {
1213               fCPV[i][j][k] = 0.; 
1214               fPRE[i][j][k] = 0.; 
1215               fPRETrackNo[i][j][k] = 0;
1216               fCPVTrackNo[i][j][k] = 0;
1217             }
1218         }
1219     }
1220 }
1221 //____________________________________________________________________________
1222
1223 void AliPMDDigitizer::UnLoad(Option_t *option)
1224 {
1225   // Unloads all the root files
1226   //
1227   const char *cS = strstr(option,"S");
1228   const char *cD = strstr(option,"D");
1229
1230   fRunLoader->UnloadgAlice();
1231   fRunLoader->UnloadHeader();
1232   fRunLoader->UnloadKinematics();
1233
1234   if (cS)
1235     {
1236       fPMDLoader->UnloadHits();
1237     }
1238   if (cD)
1239     {
1240       fPMDLoader->UnloadHits();
1241       fPMDLoader->UnloadSDigits();
1242     }
1243 }