c95006ed3a478c35aba7af2c9cafca77a4014222
[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   fpw = fopen("digits_test.dat","w");
751   fDebug = 0;
752   const char *cdeb = strstr(option,"deb");
753   if(cdeb)
754     {
755       cout << "**************** PMD Exec *************** " << endl;
756       fDebug = 1;
757     }
758   
759   Int_t ninputs = fManager->GetNinputs();
760   if(fDebug)
761     {
762       cout << " Number of files = " << ninputs << endl;
763     }
764   ResetCellADC();
765
766   for (Int_t i = 0; i < ninputs; i++)
767     {
768       Int_t troffset = fManager->GetMask(i);
769       MergeSDigits(i, troffset);
770     }
771
772   fRunLoader = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
773   fPMD  = (AliPMD*)gAlice->GetDetector("PMD");
774   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
775   if (fPMDLoader == 0x0)
776     {
777       cerr<<"AliPMDDigitizer::Exec : Can not find PMD or PMDLoader\n";
778     }
779   fPMDLoader->LoadDigits("update");
780   TTree* treeD = fPMDLoader->TreeD();
781   if (treeD == 0x0)
782     {
783       fPMDLoader->MakeTree("D");
784       treeD = fPMDLoader->TreeD();
785     }
786   Int_t bufsize = 16000;
787   if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
788   treeD->Branch("PMDDigit", &fDigits, bufsize); 
789
790   Float_t adc;
791   Float_t deltaE = 0.;
792   Int_t detno = 0;
793   Int_t trno = 1;
794
795   for (Int_t idet = 0; idet < 2; idet++)
796     {
797       for (Int_t ism = 0; ism < fgkTotUM; ism++)
798         {
799           for (Int_t jrow = 0; jrow < fgkRow; jrow++)
800             {
801               for (Int_t kcol = 0; kcol < fgkCol; kcol++)
802                 {
803                   if (idet == 0)
804                     {
805                       deltaE = fPRE[ism][jrow][kcol];
806                       trno   = fPRETrackNo[ism][jrow][kcol];
807                       detno = 0;
808                     }
809                   else if (idet == 1)
810                     {
811                       deltaE = fCPV[ism][jrow][kcol];
812                       trno   = fCPVTrackNo[ism][jrow][kcol];
813                       detno = 1;
814                     }
815                   if (deltaE > 0.)
816                     {
817                       MeV2ADC(deltaE,adc);
818                       AddDigit(trno,detno,ism,jrow,kcol,adc);
819                     }
820                 } // column loop
821             } // row    loop
822           treeD->Fill();
823           ResetDigit();
824         } // supermodule loop
825     } // detector loop
826   fPMDLoader->WriteDigits("OVERWRITE");  
827   fPMDLoader->UnloadDigits();
828   ResetCellADC();
829 }
830 //____________________________________________________________________________
831
832 void AliPMDDigitizer::MergeSDigits(Int_t filenumber, Int_t troffset)
833 {
834   // merging sdigits
835   fRunLoader = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(filenumber));
836   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
837   fPMDLoader->LoadSDigits("read");
838   TTree* treeS = fPMDLoader->TreeS();
839   AliPMDsdigit  *pmdsdigit;
840   TBranch *branch = treeS->GetBranch("PMDSDigit");
841   if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
842   branch->SetAddress(&fSDigits);
843
844   Int_t   itrackno, idet, ism;
845   Int_t   ixp, iyp;
846   Float_t edep;
847   
848   Int_t nmodules = (Int_t) treeS->GetEntries();
849   if(fDebug)
850     {
851       cout << " nmodules = " << nmodules << endl;
852       cout << " tr offset = " << troffset << endl;
853     }
854   for (Int_t imodule = 0; imodule < nmodules; imodule++)
855     {
856       treeS->GetEntry(imodule); 
857       Int_t nentries = fSDigits->GetLast();
858       if(fDebug)
859         {
860           cout << " nentries = " << nentries << endl;
861         }
862       for (Int_t ient = 0; ient < nentries+1; ient++)
863         {
864           pmdsdigit = (AliPMDsdigit*)fSDigits->UncheckedAt(ient);
865           itrackno  = pmdsdigit->GetTrackNumber();
866           idet      = pmdsdigit->GetDetector();
867           ism       = pmdsdigit->GetSMNumber();
868           ixp       = pmdsdigit->GetRow();
869           iyp       = pmdsdigit->GetColumn();
870           edep      = pmdsdigit->GetCellEdep();
871
872           if (idet == 0)
873             {
874               if (fPRE[ism][ixp][iyp] < edep)
875                 {
876                   fPRETrackNo[ism][ixp][iyp] = troffset + itrackno;
877                 }
878               fPRE[ism][ixp][iyp] += edep;
879             }
880           else if (idet == 1)
881             {
882               if (fCPV[ism][ixp][iyp] < edep)
883                 {
884                   fCPVTrackNo[ism][ixp][iyp] = troffset + itrackno;
885                 }
886               fCPV[ism][ixp][iyp] += edep;
887             }
888         }
889     }
890
891 }
892 // ----------------------------------------------------------------------
893 void AliPMDDigitizer::TrackAssignment2Cell()
894 {
895   // 
896   // This block assigns the cell id when there are
897   // multiple tracks in a cell according to the
898   // energy deposition
899   //
900   Bool_t jsort = false;
901
902   Int_t i, j, k;
903
904   Float_t *fracEdp;
905   Float_t *trEdp;
906   Int_t *status1;
907   Int_t *status2;
908   Int_t *trnarray;
909   Int_t   ****pmdTrack;
910   Float_t ****pmdEdep;
911
912   pmdTrack = new Int_t ***[fgkTotUM];
913   pmdEdep  = new Float_t ***[fgkTotUM];
914   for (i=0; i<fgkTotUM; i++)
915     {
916       pmdTrack[i] = new Int_t **[fgkRow];
917       pmdEdep[i]  = new Float_t **[fgkRow];
918     }
919
920   for (i = 0; i < fgkTotUM; i++)
921     {
922       for (j = 0; j < fgkRow; j++)
923         {
924           pmdTrack[i][j] = new Int_t *[fgkCol];
925           pmdEdep[i][j]  = new Float_t *[fgkCol];
926         }
927     }
928   
929   for (i = 0; i < fgkTotUM; i++)
930     {
931       for (j = 0; j < fgkRow; j++)
932         {
933           for (k = 0; k < fgkCol; k++)
934             {
935               Int_t nn = fPRECounter[i][j][k];
936               if(nn > 0)
937                 {
938                   pmdTrack[i][j][k] = new Int_t[nn];
939                   pmdEdep[i][j][k] = new Float_t[nn];
940                 }
941               else
942                 {
943                   nn = 1;
944                   pmdTrack[i][j][k] = new Int_t[nn];
945                   pmdEdep[i][j][k] = new Float_t[nn];
946                 }                     
947               fPRECounter[i][j][k] = 0;
948             }
949         }
950     }
951
952
953   Int_t nentries = fCell.GetEntries();
954
955   Int_t   mtrackno, ism, ixp, iyp;
956   Float_t edep;
957
958   for (i = 0; i < nentries; i++)
959     {
960       AliPMDcell* cell = (AliPMDcell*)fCell.UncheckedAt(i);
961       
962       mtrackno = cell->GetTrackNumber();
963       ism      = cell->GetSMNumber();
964       ixp      = cell->GetX();
965       iyp      = cell->GetY();
966       edep     = cell->GetEdep();
967       Int_t nn = fPRECounter[ism][ixp][iyp];
968       //      cout << " nn = " << nn << endl;
969       pmdTrack[ism][ixp][iyp][nn] = (Int_t) mtrackno;
970       pmdEdep[ism][ixp][iyp][nn] = edep;
971       fPRECounter[ism][ixp][iyp]++;
972     }
973   
974   Int_t iz, il;
975   Int_t im, ix, iy;
976   Int_t nn;
977   
978   for (im=0; im<fgkTotUM; im++)
979     {
980       for (ix=0; ix<fgkRow; ix++)
981         {
982           for (iy=0; iy<fgkCol; iy++)
983             {
984               nn = fPRECounter[im][ix][iy];
985               if (nn > 1)
986                 {
987                   // This block handles if a cell is fired
988                   // many times by many tracks
989                   status1  = new Int_t[nn];
990                   status2  = new Int_t[nn];
991                   trnarray = new Int_t[nn];
992                   for (iz = 0; iz < nn; iz++)
993                     {
994                       status1[iz] = pmdTrack[im][ix][iy][iz];
995                     }
996                   TMath::Sort(nn,status1,status2,jsort);
997                   Int_t trackOld = -99999;
998                   Int_t track, trCount = 0;
999                   for (iz = 0; iz < nn; iz++)
1000                     {
1001                       track = status1[status2[iz]];
1002                       if (trackOld != track)
1003                         {
1004                           trnarray[trCount] = track;
1005                           trCount++;
1006                         }                             
1007                       trackOld = track;
1008                     }
1009                   delete [] status1;
1010                   delete [] status2;
1011                   Float_t totEdp = 0.;
1012                   trEdp = new Float_t[trCount];
1013                   fracEdp = new Float_t[trCount];
1014                   for (il = 0; il < trCount; il++)
1015                     {
1016                       trEdp[il] = 0.;
1017                       track = trnarray[il];
1018                       for (iz = 0; iz < nn; iz++)
1019                         {
1020                           if (track == pmdTrack[im][ix][iy][iz])
1021                             {
1022                               trEdp[il] += pmdEdep[im][ix][iy][iz];
1023                             }
1024                         }
1025                       totEdp += trEdp[il];
1026                     }
1027                   Int_t ilOld = 0;
1028                   Float_t fracOld = 0.;
1029                   
1030                   for (il = 0; il < trCount; il++)
1031                     {
1032                       fracEdp[il] = trEdp[il]/totEdp;
1033                       if (fracOld < fracEdp[il])
1034                         {
1035                           fracOld = fracEdp[il];
1036                           ilOld = il;
1037                         }
1038                     }
1039                   fPRETrackNo[im][ix][iy] = trnarray[ilOld];
1040                   delete [] fracEdp;
1041                   delete [] trEdp;
1042                   delete [] trnarray;
1043                 }
1044               else if (nn == 1)
1045                 {
1046                   // This only handles if a cell is fired
1047                   // by only one track
1048                   
1049                   fPRETrackNo[im][ix][iy] = pmdTrack[im][ix][iy][0];
1050                   
1051                 }
1052               else if (nn ==0)
1053                 {
1054                   // This is if no cell is fired
1055                   fPRETrackNo[im][ix][iy] = -999;
1056                 }
1057             } // end of iy
1058         } // end of ix
1059     } // end of im
1060   
1061   // Delete all the pointers
1062   
1063   for (i = 0; i < fgkTotUM; i++)
1064     {
1065       for (j = 0; j < fgkRow; j++)
1066         {
1067           for (k = 0; k < fgkCol; k++)
1068             {
1069               delete [] pmdTrack[i][j][k];
1070               delete [] pmdEdep[i][j][k];
1071             }
1072         }
1073     }
1074   
1075   for (i = 0; i < fgkTotUM; i++)
1076     {
1077       for (j = 0; j < fgkRow; j++)
1078         {
1079           delete [] pmdTrack[i][j];
1080           delete [] pmdEdep[i][j];
1081         }
1082     }
1083   
1084   for (i = 0; i < fgkTotUM; i++)
1085     {
1086       delete [] pmdTrack[i];
1087       delete [] pmdEdep[i];
1088     }
1089   delete [] pmdTrack;
1090   delete [] pmdEdep;
1091   // 
1092   // End of the cell id assignment
1093   //
1094 }
1095 //____________________________________________________________________________
1096 void AliPMDDigitizer::MeV2ADC(Float_t mev, Float_t & adc) const
1097 {
1098   // This converts the simulated edep to ADC according to the
1099   // Test Beam Data
1100   // To be done
1101   //
1102
1103   //  adc = mev*1.;
1104
1105   
1106   // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% //
1107   // PS Test in September 2003
1108   // MeV - ADC conversion for 10bit ADC
1109
1110   const Float_t kConstant   = 7.181;
1111   const Float_t kErConstant = 0.6899;
1112   const Float_t kSlope      = 35.93;
1113   const Float_t kErSlope    = 0.306;
1114
1115   //gRandom->SetSeed();
1116
1117   Float_t cons   = gRandom->Gaus(kConstant,kErConstant);
1118   Float_t slop   = gRandom->Gaus(kSlope,kErSlope);
1119
1120   Float_t adc10bit = slop*mev*0.001 + cons;
1121
1122   // 12 bit adc
1123
1124   Int_t adc12bit = (Int_t) (4.0*adc10bit);
1125
1126   if(adc12bit < 3000)
1127     {
1128       adc = (Float_t) adc12bit;
1129     }
1130   else if (adc12bit >= 3000)
1131     {
1132       adc = 3000.0;
1133     }
1134
1135   // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% //
1136
1137 }
1138 //____________________________________________________________________________
1139 void AliPMDDigitizer::AddSDigit(Int_t trnumber, Int_t det, Int_t smnumber, 
1140   Int_t irow, Int_t icol, Float_t adc)
1141 {
1142   // Add SDigit
1143   //
1144   if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
1145   TClonesArray &lsdigits = *fSDigits;
1146   new(lsdigits[fNsdigit++])  AliPMDsdigit(trnumber,det,smnumber,irow,icol,adc);
1147 }
1148 //____________________________________________________________________________
1149
1150 void AliPMDDigitizer::AddDigit(Int_t trnumber, Int_t det, Int_t smnumber, 
1151   Int_t irow, Int_t icol, Float_t adc)
1152 {
1153   // Add Digit
1154   //
1155   if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
1156   TClonesArray &ldigits = *fDigits;
1157   new(ldigits[fNdigit++]) AliPMDdigit(trnumber,det,smnumber,irow,icol,adc);
1158 }
1159 //____________________________________________________________________________
1160
1161 void AliPMDDigitizer::SetZPosition(Float_t zpos)
1162 {
1163   fZPos = zpos;
1164 }
1165 //____________________________________________________________________________
1166 Float_t AliPMDDigitizer::GetZPosition() const
1167 {
1168   return fZPos;
1169 }
1170 //____________________________________________________________________________
1171
1172 void AliPMDDigitizer::ResetCell()
1173 {
1174   // clears the cell array and also the counter
1175   //  for each cell
1176   //
1177   fCell.Delete();
1178   for (Int_t i = 0; i < fgkTotUM; i++)
1179     {
1180       for (Int_t j = 0; j < fgkRow; j++)
1181         {
1182           for (Int_t k = 0; k < fgkCol; k++)
1183             {
1184               fPRECounter[i][j][k] = 0; 
1185             }
1186         }
1187     }
1188 }
1189 //____________________________________________________________________________
1190 void AliPMDDigitizer::ResetSDigit()
1191 {
1192   // Clears SDigits
1193   fNsdigit = 0;
1194   if (fSDigits) fSDigits->Delete();
1195 }
1196 //____________________________________________________________________________
1197 void AliPMDDigitizer::ResetDigit()
1198 {
1199   // Clears Digits
1200   fNdigit = 0;
1201   if (fDigits) fDigits->Delete();
1202 }
1203 //____________________________________________________________________________
1204
1205 void AliPMDDigitizer::ResetCellADC()
1206 {
1207   // Clears individual cells edep and track number
1208   for (Int_t i = 0; i < fgkTotUM; i++)
1209     {
1210       for (Int_t j = 0; j < fgkRow; j++)
1211         {
1212           for (Int_t k = 0; k < fgkCol; k++)
1213             {
1214               fCPV[i][j][k] = 0.; 
1215               fPRE[i][j][k] = 0.; 
1216               fPRETrackNo[i][j][k] = 0;
1217               fCPVTrackNo[i][j][k] = 0;
1218             }
1219         }
1220     }
1221 }
1222 //____________________________________________________________________________
1223
1224 void AliPMDDigitizer::UnLoad(Option_t *option)
1225 {
1226   // Unloads all the root files
1227   //
1228   const char *cS = strstr(option,"S");
1229   const char *cD = strstr(option,"D");
1230
1231   fRunLoader->UnloadgAlice();
1232   fRunLoader->UnloadHeader();
1233   fRunLoader->UnloadKinematics();
1234
1235   if (cS)
1236     {
1237       fPMDLoader->UnloadHits();
1238     }
1239   if (cD)
1240     {
1241       fPMDLoader->UnloadHits();
1242       fPMDLoader->UnloadSDigits();
1243     }
1244 }