CPV or PREshower plane is according to Module number instead of zposition
[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 #include <TRandom.h>
34
35 #include "AliLog.h"
36 #include "AliRun.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 #include "AliCDBManager.h"
47 #include "AliCDBStorage.h"
48 #include "AliCDBEntry.h"
49 #include "AliMC.h"
50
51 #include "AliPMD.h"
52 #include "AliPMDhit.h"
53 #include "AliPMDcell.h"
54 #include "AliPMDsdigit.h"
55 #include "AliPMDdigit.h"
56 #include "AliPMDCalibData.h"
57 #include "AliPMDPedestal.h"
58 #include "AliPMDDigitizer.h"
59
60
61 ClassImp(AliPMDDigitizer)
62
63 AliPMDDigitizer::AliPMDDigitizer() :
64   fRunLoader(0),
65   fPMDHit(0),
66   fPMD(0),
67   fPMDLoader(0),
68   fCalibGain(GetCalibGain()),
69   fCalibPed(GetCalibPed()),
70   fSDigits(0),
71   fDigits(0),
72   fCPVCell(0),
73   fCell(0),
74   fNsdigit(0),
75   fNdigit(0),
76   fDetNo(0),
77   fZPos(361.5)   // in units of cm, default position of PMD
78 {
79   // Default Constructor
80   //
81   for (Int_t i = 0; i < fgkTotUM; i++)
82     {
83       for (Int_t j = 0; j < fgkRow; j++)
84         {
85           for (Int_t k = 0; k < fgkCol; k++)
86             {
87               fCPV[i][j][k]         = 0.;
88               fPRE[i][j][k]         = 0.;
89               fCPVCounter[i][j][k]  =  0; 
90               fPRECounter[i][j][k]  =  0;
91               fCPVTrackNo[i][j][k]  = -1;
92               fPRETrackNo[i][j][k]  = -1;
93               fCPVTrackPid[i][j][k] = -1;
94               fPRETrackPid[i][j][k] = -1;
95             }
96         }
97     }
98
99 }
100 //____________________________________________________________________________
101 AliPMDDigitizer::AliPMDDigitizer(const AliPMDDigitizer& digitizer):
102   AliDigitizer(digitizer),
103   fRunLoader(0),
104   fPMDHit(0),
105   fPMD(0),
106   fPMDLoader(0),
107   fCalibGain(GetCalibGain()),
108   fCalibPed(GetCalibPed()),
109   fSDigits(0),
110   fDigits(0),
111   fCPVCell(0),
112   fCell(0),
113   fNsdigit(0),
114   fNdigit(0),
115   fDetNo(0),
116   fZPos(361.5)   // in units of cm, default position of PMD
117 {
118   // copy constructor
119   AliError("Copy constructor not allowed ");
120   
121 }
122 //____________________________________________________________________________
123 AliPMDDigitizer & AliPMDDigitizer::operator=(const AliPMDDigitizer& /*digitizer*/)
124 {
125   // Assignment operator
126   AliError("Assignement operator not allowed ");
127
128   return *this;
129 }
130 //____________________________________________________________________________
131 AliPMDDigitizer::AliPMDDigitizer(AliRunDigitizer* manager):
132   AliDigitizer(manager),
133   fRunLoader(0),
134   fPMDHit(0),
135   fPMD(0),
136   fPMDLoader(0),
137   fCalibGain(GetCalibGain()),
138   fCalibPed(GetCalibPed()),
139   fSDigits(new TClonesArray("AliPMDsdigit", 1000)),
140   fDigits(new TClonesArray("AliPMDdigit", 1000)),
141   fCPVCell(0),
142   fCell(0),
143   fNsdigit(0),
144   fNdigit(0),
145   fDetNo(0),
146   fZPos(361.5)// in units of cm, This is the default position of PMD
147 {
148   // ctor which should be used
149
150
151   for (Int_t i = 0; i < fgkTotUM; i++)
152     {
153       for (Int_t j = 0; j < fgkRow; j++)
154         {
155           for (Int_t k = 0; k < fgkCol; k++)
156             {
157               fCPV[i][j][k]         = 0.;
158               fPRE[i][j][k]         = 0.;
159               fCPVCounter[i][j][k]  =  0; 
160               fPRECounter[i][j][k]  =  0;
161               fCPVTrackNo[i][j][k]  = -1;
162               fPRETrackNo[i][j][k]  = -1;
163               fCPVTrackPid[i][j][k] = -1;
164               fPRETrackPid[i][j][k] = -1;
165             }
166         }
167     }
168 }
169
170 //____________________________________________________________________________
171 AliPMDDigitizer::~AliPMDDigitizer()
172 {
173   // Default Destructor
174   //
175   if (fSDigits) {
176     fSDigits->Delete();
177     delete fSDigits;
178     fSDigits=0;
179   }
180   if (fDigits) {
181     fDigits->Delete();
182     delete fDigits;
183     fDigits=0;
184   }
185   fCPVCell.Delete();
186   fCell.Delete();
187 }
188 //
189 // Member functions
190 //
191 //____________________________________________________________________________
192 void AliPMDDigitizer::OpengAliceFile(const char *file, Option_t *option)
193 {
194   // Loads galice.root file and corresponding header, kinematics
195   // hits and sdigits or digits depending on the option
196   //
197
198   TString evfoldname = AliConfig::GetDefaultEventFolderName();
199   fRunLoader = AliRunLoader::GetRunLoader(evfoldname);
200   if (!fRunLoader)
201       fRunLoader = AliRunLoader::Open(file,AliConfig::GetDefaultEventFolderName(), "UPDATE");
202   
203   if (!fRunLoader)
204    {
205      AliError(Form("Can not open session for file %s.",file));
206    }
207
208   const char *cHS = strstr(option,"HS");
209   const char *cHD = strstr(option,"HD");
210   const char *cSD = strstr(option,"SD");
211   
212   if(cHS || cHD)
213     {
214       if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
215       if (!fRunLoader->TreeE()) fRunLoader->LoadHeader();
216       if (!fRunLoader->TreeK()) fRunLoader->LoadKinematics();
217   
218       gAlice = fRunLoader->GetAliRun();
219   
220       if (gAlice)
221         {
222           AliDebug(1,"Alirun object found");
223         }
224       else
225         {
226           AliError("Could not found Alirun object");
227         }
228   
229       fPMD  = (AliPMD*)gAlice->GetDetector("PMD");
230     }
231
232   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
233   if (fPMDLoader == 0x0)
234     {
235       AliError("Can not find PMDLoader");
236     }
237
238
239   if (cHS)
240     {
241       fPMDLoader->LoadHits("READ");
242       fPMDLoader->LoadSDigits("recreate");
243     }
244   else if (cHD)
245     {
246       fPMDLoader->LoadHits("READ");
247       fPMDLoader->LoadDigits("recreate");
248     }
249   else if (cSD)
250     {
251       fPMDLoader->LoadSDigits("READ");
252       fPMDLoader->LoadDigits("recreate");
253     }
254 }
255 //____________________________________________________________________________
256 void AliPMDDigitizer::Hits2SDigits(Int_t ievt)
257 {
258   // This reads the PMD Hits tree and assigns the right track number
259   // to a cell and stores in the summable digits tree
260   //
261
262   const Int_t kPi0 = 111;
263   const Int_t kGamma = 22;
264   Int_t npmd;
265   Int_t trackno;
266   Int_t smnumber;
267   Int_t trackpid;
268   Int_t mtrackno;
269   Int_t mtrackpid;
270
271   Float_t xPos, yPos, zPos;
272   Int_t xpad = -1, ypad = -1;
273   Float_t edep;
274   Float_t vx = -999.0, vy = -999.0, vz = -999.0;
275
276
277   if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
278   ResetSDigit();
279
280   AliDebug(1,Form("Event Number = %d",ievt));
281   Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
282   AliDebug(1,Form("Number of Particles = %d",nparticles));
283
284   //
285
286   fRunLoader->GetEvent(ievt);
287   // ------------------------------------------------------- //
288   // Pointer to specific detector hits.
289   // Get pointers to Alice detectors and Hits containers
290
291   TTree* treeH = fPMDLoader->TreeH();
292   
293   Int_t ntracks    = (Int_t) treeH->GetEntries();
294   AliDebug(1,Form("Number of Tracks in the TreeH = %d", ntracks));
295   TTree* treeS = fPMDLoader->TreeS();
296   if (treeS == 0x0)
297     {
298       fPMDLoader->MakeTree("S");
299       treeS = fPMDLoader->TreeS();
300     }
301   Int_t bufsize = 16000;
302   treeS->Branch("PMDSDigit", &fSDigits, bufsize);
303   
304   TClonesArray* hits = 0;
305   if (fPMD) hits = fPMD->Hits();
306
307   // Start loop on tracks in the hits containers
308
309   for (Int_t track=0; track<ntracks;track++)
310     {
311       gAlice->GetMCApp()->ResetHits();
312       treeH->GetEvent(track);
313       if (fPMD)
314         {
315           npmd = hits->GetEntriesFast();
316           for (Int_t ipmd = 0; ipmd < npmd; ipmd++)
317             {
318               fPMDHit = (AliPMDhit*) hits->UncheckedAt(ipmd);
319               trackno = fPMDHit->GetTrack();
320               //  get kinematics of the particles
321
322               TParticle* mparticle = gAlice->GetMCApp()->Particle(trackno);
323               trackpid  = mparticle->GetPdgCode();
324               Int_t  ks = mparticle->GetStatusCode();
325               Int_t imo;
326               Int_t tracknoOld=0, trackpidOld=0, statusOld = 0;
327               
328               if (mparticle->GetFirstMother() == -1)
329                 {
330                   tracknoOld  = trackno;
331                   trackpidOld = trackpid;
332                   statusOld   = -1;
333                 }
334               Int_t igstatus = 0;
335
336               Int_t trnotemp = trackno;    // Modified on 25th Nov 2009
337               if(ks==1||(imo = mparticle->GetFirstMother())<0 ){
338                 vx = mparticle->Vx();
339                 vy = mparticle->Vy();
340                 vz = mparticle->Vz();
341                 
342                 if(trackpid==kGamma||trackpid==11||trackpid==-11||
343                    trackpid==kPi0)igstatus=1;
344               }
345               
346               
347               while(((imo = mparticle->GetFirstMother()) >= 0)&& 
348                     (ks = mparticle->GetStatusCode() <1) )
349                 {
350                   mparticle =  gAlice->GetMCApp()->Particle(imo);
351                   trackpid = mparticle->GetPdgCode();
352                   ks = mparticle->GetStatusCode();
353                   vx = mparticle->Vx();
354                   vy = mparticle->Vy();
355                   vz = mparticle->Vz();
356                   
357                   // Modified on 25th Nov 2009
358
359                   trnotemp = trackno;
360                   if(trackpid == 111)
361                     {
362                       trackno = trnotemp;
363                     }
364                   if(trackpid != 111)
365                     {
366                       trackno=imo;
367                     }
368                   // end of modification on 25th Nov 2009
369                 }
370          
371               if(trackpid==kGamma||trackpid==11||trackpid==-11||
372                  trackpid==kPi0)igstatus=1;
373               mtrackpid=trackpid;
374               mtrackno=trackno;
375               trackpid=trackpidOld;
376               trackno=tracknoOld;
377               
378               //-----------------end of modification----------------
379               xPos = fPMDHit->X();
380               yPos = fPMDHit->Y();
381               zPos = fPMDHit->Z();
382               
383               edep       = fPMDHit->GetEnergy();
384               Int_t vol1 = fPMDHit->GetVolume(1); // Column
385               Int_t vol2 = fPMDHit->GetVolume(2); // Row
386               Int_t vol7 = fPMDHit->GetVolume(4); // Serial Module No
387
388
389               // -----------------------------------------//
390               // In new geometry after adding electronics //
391               // For Super Module 1 & 2                   //
392               //  nrow = 48, ncol = 96                    //
393               // For Super Module 3 & 4                   //
394               //  nrow = 96, ncol = 48                    //
395               // -----------------------------------------//
396
397               if (vol7 < 24)
398                 {
399                   smnumber = vol7;
400                 }
401               else
402                 {
403                   smnumber = vol7 - 24;
404                 }
405               Int_t vol8 = smnumber/6 + 1;   // fake supermodule
406
407               if (vol8 == 1 || vol8 == 2)
408                 {
409                   xpad = vol2;
410                   ypad = vol1;
411                 }
412               else if (vol8 == 3 || vol8 == 4)
413                 {
414                   xpad = vol1;
415                   ypad = vol2;
416                 }
417
418               AliDebug(2,Form("Zposition = %f Edeposition = %f",zPos,edep));
419
420               if (vol7 < 24)
421                 {
422                   // PRE
423                   fDetNo = 0;
424                 }
425               else
426                 {
427                   // CPV
428                   fDetNo = 1;
429                 }
430
431               Int_t smn = smnumber;
432               Int_t ixx = xpad     - 1;
433               Int_t iyy = ypad     - 1;
434               if (fDetNo == 0)
435                 {
436                   fPRE[smn][ixx][iyy] += edep;
437                   fPRECounter[smn][ixx][iyy]++;
438
439                   AliPMDcell* cell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
440                   fCell.Add(cell);
441                 }
442               else if(fDetNo == 1)
443                 {
444                   fCPV[smn][ixx][iyy] += edep;
445                   fCPVCounter[smn][ixx][iyy]++;
446                   AliPMDcell* cpvcell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep); 
447                   fCPVCell.Add(cpvcell);
448                 }
449             }
450         }
451     } // Track Loop ended
452   TrackAssignment2CPVCell();
453   TrackAssignment2Cell();
454   ResetCell();
455
456   Float_t deltaE      = 0.;
457   Int_t   detno       = 0;
458   Int_t   trno        = -1;
459   Int_t   trpid       = -99;
460
461   for (Int_t idet = 0; idet < 2; idet++)
462     {
463       for (Int_t ism = 0; ism < fgkTotUM; ism++)
464         {
465           for (Int_t jrow = 0; jrow < fgkRow; jrow++)
466             {
467               for (Int_t kcol = 0; kcol < fgkCol; kcol++)
468                 {
469                   if (idet == 0)
470                     {
471                       deltaE = fPRE[ism][jrow][kcol];
472                       trno   = fPRETrackNo[ism][jrow][kcol];
473                       detno = 0;
474                     }
475                   else if (idet == 1)
476                     {
477                       deltaE = fCPV[ism][jrow][kcol];
478                       trno   = fCPVTrackNo[ism][jrow][kcol];
479                       detno = 1;
480                     }
481                   if (deltaE > 0.)
482                     {
483                       // Natasha
484                       TParticle *mparticle = gAlice->GetMCApp()->Particle(trno);
485                       trpid = mparticle->GetPdgCode();
486                       AddSDigit(trno,trpid,detno,ism,jrow,kcol,deltaE);
487                     }
488                 }
489             }
490           treeS->Fill();
491           ResetSDigit();
492         }
493     }
494   fPMDLoader->WriteSDigits("OVERWRITE");
495   ResetCellADC();
496 }
497 //____________________________________________________________________________
498
499 void AliPMDDigitizer::Hits2Digits(Int_t ievt)
500 {
501   // This reads the PMD Hits tree and assigns the right track number
502   // to a cell and stores in the digits tree
503   //
504   const Int_t kPi0 = 111;
505   const Int_t kGamma = 22;
506   Int_t npmd;
507   Int_t trackno;
508   Int_t smnumber;
509   Int_t trackpid;
510   Int_t mtrackno;
511   Int_t mtrackpid;
512
513   Float_t xPos, yPos, zPos;
514   Int_t xpad = -1, ypad = -1;
515   Float_t edep;
516   Float_t vx = -999.0, vy = -999.0, vz = -999.0;
517
518   if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
519   ResetDigit();
520
521   AliDebug(1,Form("Event Number =  %d",ievt));
522   Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
523   AliDebug(1,Form("Number of Particles = %d", nparticles));
524
525   fRunLoader->GetEvent(ievt);
526   // ------------------------------------------------------- //
527   // Pointer to specific detector hits.
528   // Get pointers to Alice detectors and Hits containers
529
530   fPMD  = (AliPMD*)gAlice->GetDetector("PMD");
531   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
532
533   if (fPMDLoader == 0x0)
534     {
535       AliError("Can not find PMD or PMDLoader");
536     }
537   TTree* treeH = fPMDLoader->TreeH();
538   Int_t ntracks    = (Int_t) treeH->GetEntries();
539   AliDebug(1,Form("Number of Tracks in the TreeH = %d", ntracks));
540   fPMDLoader->LoadDigits("recreate");
541   TTree* treeD = fPMDLoader->TreeD();
542   if (treeD == 0x0)
543     {
544       fPMDLoader->MakeTree("D");
545       treeD = fPMDLoader->TreeD();
546     }
547   Int_t bufsize = 16000;
548   treeD->Branch("PMDDigit", &fDigits, bufsize);
549   
550   TClonesArray* hits = 0;
551   if (fPMD) hits = fPMD->Hits();
552
553   // Start loop on tracks in the hits containers
554
555   for (Int_t track=0; track<ntracks;track++)
556     {
557       gAlice->GetMCApp()->ResetHits();
558       treeH->GetEvent(track);
559       
560       if (fPMD)
561         {
562           npmd = hits->GetEntriesFast();
563           for (Int_t ipmd = 0; ipmd < npmd; ipmd++)
564             {
565               fPMDHit = (AliPMDhit*) hits->UncheckedAt(ipmd);
566               trackno = fPMDHit->GetTrack();
567               
568               //  get kinematics of the particles
569               
570               TParticle* mparticle = gAlice->GetMCApp()->Particle(trackno);
571               trackpid  = mparticle->GetPdgCode();
572               Int_t  ks = mparticle->GetStatusCode();
573               Int_t imo;
574               Int_t tracknoOld=0, trackpidOld=0, statusOld = 0;
575               if (mparticle->GetFirstMother() == -1)
576                 {
577                   tracknoOld  = trackno;
578                   trackpidOld = trackpid;
579                   statusOld   = -1;
580                 }
581
582               Int_t igstatus = 0;
583
584               Int_t trnotemp = trackno;   // modified on 25th Nov 2009
585               if(ks==1||(imo = mparticle->GetFirstMother())<0 ){
586                 vx = mparticle->Vx();
587                 vy = mparticle->Vy();
588                 vz = mparticle->Vz();
589                 
590                 if(trackpid==kGamma||trackpid==11||trackpid==-11||trackpid==kPi0)
591                   igstatus=1;
592               }
593               
594               
595               while(((imo = mparticle->GetFirstMother()) >= 0)&& 
596                     (ks = mparticle->GetStatusCode() <1) )
597                 {
598                   mparticle =  gAlice->GetMCApp()->Particle(imo);
599                   trackpid = mparticle->GetPdgCode();
600                   ks = mparticle->GetStatusCode();
601                   vx = mparticle->Vx();
602                   vy = mparticle->Vy();
603                   vz = mparticle->Vz();
604                   
605                   // Modified on 25th Nov 2009
606                   
607                   trnotemp = trackno;
608                   if(trackpid == 111)
609                     {
610                       trackno = trnotemp;
611                     }
612                   if(trackpid != 111)
613                     {
614                       trackno=imo;
615                     }
616                 }
617          
618               if(trackpid==kGamma||trackpid==11||trackpid==-11||trackpid==kPi0)
619                 igstatus=1;
620               mtrackpid=trackpid;
621               mtrackno=trackno;
622               trackpid=trackpidOld;
623               trackno=tracknoOld;
624               
625               xPos = fPMDHit->X();
626               yPos = fPMDHit->Y();
627               zPos = fPMDHit->Z();
628               edep       = fPMDHit->GetEnergy();
629               Int_t vol1 = fPMDHit->GetVolume(1); // Column
630               Int_t vol2 = fPMDHit->GetVolume(2); // Row
631               Int_t vol7 = fPMDHit->GetVolume(4); // Serial Module No
632
633               // -----------------------------------------//
634               // In new geometry after adding electronics //
635               // For Super Module 1 & 2                   //
636               //  nrow = 48, ncol = 96                    //
637               // For Super Module 3 & 4                   //
638               //  nrow = 96, ncol = 48                    //
639               // -----------------------------------------//
640               
641               if (vol7 < 24)
642                 {
643                   smnumber = vol7;
644                 }
645               else
646                 {
647                   smnumber = vol7 - 24;
648                 }
649               Int_t vol8 = smnumber/6 + 1;    // fake supermodule
650
651               if (vol8 == 1 || vol8 == 2)
652                 {
653                   xpad = vol2;
654                   ypad = vol1;
655                 }
656               else if (vol8 == 3 || vol8 == 4)
657                 {
658                   xpad = vol1;
659                   ypad = vol2;
660                 }
661
662               AliDebug(2,Form("ZPosition = %f Edeposition = %d",zPos,edep));
663
664               if (vol7 < 24)
665                 {
666                   // PRE
667                   fDetNo = 0;
668                 }
669               else
670                 {
671                   fDetNo = 1;
672                 }
673               /*
674                 if (zPos < fZPos)
675                 {
676                 // CPV
677                 fDetNo = 1;
678                 }
679                 else if (zPos > fZPos)
680                 {
681                 // PMD
682                 fDetNo = 0;
683                 }
684               */
685               Int_t smn = smnumber;
686               Int_t ixx = xpad     - 1;
687               Int_t iyy = ypad     - 1;
688               if (fDetNo == 0)
689                 {
690                   fPRE[smn][ixx][iyy] += edep;
691                   fPRECounter[smn][ixx][iyy]++;
692
693                   AliPMDcell* cell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
694
695                   fCell.Add(cell);
696                 }
697               else if(fDetNo == 1)
698                 {
699                   fCPV[smn][ixx][iyy] += edep;
700                   fCPVCounter[smn][ixx][iyy]++;
701                   AliPMDcell* cpvcell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep); 
702                   fCPVCell.Add(cpvcell);
703                 }
704             }
705         }
706     } // Track Loop ended
707   TrackAssignment2CPVCell();
708   TrackAssignment2Cell();
709   ResetCell();
710
711   Float_t gain1;
712   Float_t adc;
713   Float_t deltaE = 0.;
714   Int_t detno = 0;
715   Int_t trno = 1;
716   Int_t trpid = -99;
717
718   for (Int_t idet = 0; idet < 2; idet++)
719   {
720       for (Int_t ism = 0; ism < fgkTotUM; ism++)
721       {
722           for (Int_t jrow = 0; jrow < fgkRow; jrow++)
723           {
724               for (Int_t kcol = 0; kcol < fgkCol; kcol++)
725               {
726                   if (idet == 0)
727                   {
728                       deltaE = fPRE[ism][jrow][kcol];
729                       trno   = fPRETrackNo[ism][jrow][kcol];
730                       detno  = 0;
731                   }
732                   else if (idet == 1)
733                   {
734                       deltaE = fCPV[ism][jrow][kcol];
735                       trno   = fCPVTrackNo[ism][jrow][kcol];
736                       detno  = 1;
737                   }
738                   if (deltaE > 0.)
739                   {
740                       MeV2ADC(deltaE,adc);
741
742                       // To decalibrate the adc values
743                       //
744                       gain1 = Gain(idet,ism,jrow,kcol);
745                       if (gain1 != 0.)
746                       {
747                           Int_t adcDecalib = (Int_t)(adc/gain1);
748                           adc = (Float_t) adcDecalib;
749                       }
750                       else if(gain1 == 0.)
751                       {
752                           adc = 0.;
753                       }
754
755                       // Pedestal Decalibration
756                       Int_t   pedmeanrms = 
757                           fCalibPed->GetPedMeanRms(idet,ism,jrow,kcol);
758                       Int_t   pedrms1    = (Int_t) pedmeanrms%100;
759                       Float_t pedrms     = (Float_t)pedrms1/10.;
760                       Float_t pedmean    = 
761                           (Float_t) (pedmeanrms - pedrms1)/1000.0;
762                       if (adc > 0.)
763                       {
764                           adc += (pedmean + 3.0*pedrms);
765                           TParticle *mparticle
766                             = gAlice->GetMCApp()->Particle(trno);
767                           trpid = mparticle->GetPdgCode();
768                           
769                           AddDigit(trno,trpid,detno,ism,jrow,kcol,adc);
770                       }
771                   }
772               } // column loop
773           } // row    loop
774           treeD->Fill();
775           ResetDigit();
776       } // supermodule loop
777   } // detector loop
778   
779   fPMDLoader->WriteDigits("OVERWRITE");
780   ResetCellADC();
781
782 }
783 //____________________________________________________________________________
784
785
786 void AliPMDDigitizer::SDigits2Digits(Int_t ievt)
787 {
788   // This reads the PMD sdigits tree and converts energy deposition
789   // in a cell to ADC and stores in the digits tree
790   //
791
792   fRunLoader->GetEvent(ievt);
793
794   TTree* treeS = fPMDLoader->TreeS();
795   AliPMDsdigit  *pmdsdigit;
796   TBranch *branch = treeS->GetBranch("PMDSDigit");
797   if(!branch)
798     {
799       AliError("PMD Sdigit branch does not exist");
800       return;
801     }
802   if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
803   branch->SetAddress(&fSDigits);
804
805   TTree* treeD = fPMDLoader->TreeD();
806   if (treeD == 0x0)
807     {
808       fPMDLoader->MakeTree("D");
809       treeD = fPMDLoader->TreeD();
810     }
811   Int_t bufsize = 16000;
812   if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
813   treeD->Branch("PMDDigit", &fDigits, bufsize);
814
815   Int_t   trno, trpid, det, smn;
816   Int_t   irow, icol;
817   Float_t edep, adc;
818
819   Int_t nmodules = (Int_t) treeS->GetEntries();
820   AliDebug(1,Form("Number of modules = %d",nmodules));
821
822   for (Int_t imodule = 0; imodule < nmodules; imodule++)
823     {
824       treeS->GetEntry(imodule);
825       Int_t nentries = fSDigits->GetLast();
826       AliDebug(2,Form("Number of entries per module = %d",nentries+1));
827       for (Int_t ient = 0; ient < nentries+1; ient++)
828         {
829           pmdsdigit = (AliPMDsdigit*)fSDigits->UncheckedAt(ient);
830           trno   = pmdsdigit->GetTrackNumber();
831           trpid  = pmdsdigit->GetTrackPid();
832           det    = pmdsdigit->GetDetector();
833           smn    = pmdsdigit->GetSMNumber();
834           irow   = pmdsdigit->GetRow();
835           icol   = pmdsdigit->GetColumn();
836           edep   = pmdsdigit->GetCellEdep();
837
838           MeV2ADC(edep,adc);
839
840           // To decalibrte the adc values
841           //
842           Float_t gain1 = Gain(det,smn,irow,icol);
843           if (gain1 != 0.)
844           {
845             Int_t adcDecalib = (Int_t)(adc/gain1);
846             adc = (Float_t) adcDecalib;
847           }
848           else if(gain1 == 0.)
849           {
850               adc = 0.;
851           }
852           // Pedestal Decalibration
853           Int_t   pedmeanrms = fCalibPed->GetPedMeanRms(det,smn,irow,icol);
854           Int_t   pedrms1    = (Int_t) pedmeanrms%100;
855           Float_t pedrms     = (Float_t)pedrms1/10.;
856           Float_t pedmean    = (Float_t) (pedmeanrms - pedrms1)/1000.0;
857           if(adc > 0.)
858           {
859               adc += (pedmean + 3.0*pedrms);
860               AddDigit(trno,trpid,det,smn,irow,icol,adc);
861           }
862
863         }
864       treeD->Fill();
865       ResetDigit();
866     }
867   fPMDLoader->WriteDigits("OVERWRITE");
868
869 }
870 //____________________________________________________________________________
871 void AliPMDDigitizer::Exec(Option_t *option)
872 {
873   // Does the event merging and digitization
874   const char *cdeb = strstr(option,"deb");
875   if(cdeb)
876     {
877       AliDebug(100," *** PMD Exec is called ***");
878     }
879
880   Int_t ninputs = fManager->GetNinputs();
881   AliDebug(1,Form("Number of files to be processed = %d",ninputs));
882   ResetCellADC();
883
884   for (Int_t i = 0; i < ninputs; i++)
885     {
886       Int_t troffset = fManager->GetMask(i);
887       MergeSDigits(i, troffset);
888     }
889
890   fRunLoader = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
891   fPMD  = (AliPMD*)gAlice->GetDetector("PMD");
892   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
893   if (fPMDLoader == 0x0)
894     {
895       AliError("Can not find PMD or PMDLoader");
896     }
897   fPMDLoader->LoadDigits("update");
898   TTree* treeD = fPMDLoader->TreeD();
899   if (treeD == 0x0)
900     {
901       fPMDLoader->MakeTree("D");
902       treeD = fPMDLoader->TreeD();
903     }
904   Int_t bufsize = 16000;
905   if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
906   treeD->Branch("PMDDigit", &fDigits, bufsize);
907
908   Float_t adc;
909   Float_t deltaE = 0.;
910   Int_t detno    = 0;
911   Int_t trno     = 1;
912   Int_t trpid    = -99;
913
914   for (Int_t idet = 0; idet < 2; idet++)
915     {
916       for (Int_t ism = 0; ism < fgkTotUM; ism++)
917         {
918           for (Int_t jrow = 0; jrow < fgkRow; jrow++)
919             {
920               for (Int_t kcol = 0; kcol < fgkCol; kcol++)
921                 {
922                   if (idet == 0)
923                     {
924                       deltaE = fPRE[ism][jrow][kcol];
925                       trno   = fPRETrackNo[ism][jrow][kcol];
926                       trpid  = fPRETrackPid[ism][jrow][kcol];
927                       detno  = 0;
928                     }
929                   else if (idet == 1)
930                     {
931                       deltaE = fCPV[ism][jrow][kcol];
932                       trno   = fCPVTrackNo[ism][jrow][kcol];
933                       trpid  = fCPVTrackPid[ism][jrow][kcol];
934                       detno  = 1;
935                     }
936                   if (deltaE > 0.)
937                     {
938                       MeV2ADC(deltaE,adc);
939
940                       //
941                       // Gain decalibration
942                       //
943                       Float_t gain1 = Gain(idet,ism,jrow,kcol);
944
945                       if (gain1 != 0.)
946                       {
947                           Int_t adcDecalib = (Int_t)(adc/gain1);
948                           adc = (Float_t) adcDecalib;
949                       }
950                       else if(gain1 == 0.)
951                       {
952                           adc = 0.;
953                       }
954                       // Pedestal Decalibration
955                       Int_t   pedmeanrms = 
956                           fCalibPed->GetPedMeanRms(idet,ism,jrow,kcol);
957                       Int_t   pedrms1    = (Int_t) pedmeanrms%100;
958                       Float_t pedrms     = (Float_t)pedrms1/10.;
959                       Float_t pedmean    = 
960                           (Float_t) (pedmeanrms - pedrms1)/1000.0;
961                       if (adc > 0.)
962                       {
963                           adc += (pedmean + 3.0*pedrms);
964                           AddDigit(trno,trpid,detno,ism,jrow,kcol,adc);
965                       }
966
967                     }
968                 } // column loop
969             } // row    loop
970           treeD->Fill();
971           ResetDigit();
972         } // supermodule loop
973     } // detector loop
974   fPMDLoader->WriteDigits("OVERWRITE");
975   fPMDLoader->UnloadDigits();
976   ResetCellADC();
977 }
978 //____________________________________________________________________________
979 void AliPMDDigitizer::TrackAssignment2CPVCell()
980 {
981   // This block assigns the cell id when there are
982   // multiple tracks in a cell according to the
983   // energy deposition
984   // This method added by Ajay
985   Bool_t jsort = false;
986
987   Int_t i, j, k;
988
989   Int_t   *status1;
990   Int_t   *status2;
991   Int_t   *trnarray;  
992   Float_t *fracEdp;
993   Float_t *trEdp;
994   
995   Int_t   ****cpvTrack;
996   Float_t ****cpvEdep;
997
998   cpvTrack = new Int_t ***[fgkTotUM];
999   cpvEdep  = new Float_t ***[fgkTotUM];
1000   for (i=0; i<fgkTotUM; i++)
1001     {
1002       cpvTrack[i] = new Int_t **[fgkRow];
1003       cpvEdep[i]  = new Float_t **[fgkRow];
1004     }
1005
1006   for (i = 0; i < fgkTotUM; i++)
1007     {
1008       for (j = 0; j < fgkRow; j++)
1009         {
1010           cpvTrack[i][j] = new Int_t *[fgkCol];
1011           cpvEdep[i][j]  = new Float_t *[fgkCol];
1012         }
1013     }
1014   for (i = 0; i < fgkTotUM; i++)
1015     {
1016       for (j = 0; j < fgkRow; j++)
1017         {
1018           for (k = 0; k < fgkCol; k++)
1019             {
1020               Int_t nn = fCPVCounter[i][j][k];
1021               if(nn > 0)
1022                 {
1023                   cpvTrack[i][j][k] = new Int_t[nn];
1024                   cpvEdep[i][j][k] = new Float_t[nn];
1025                 }
1026               else
1027                 {
1028                   nn = 1;
1029                   cpvTrack[i][j][k] = new Int_t[nn];
1030                   cpvEdep[i][j][k] = new Float_t[nn];
1031                 }                     
1032               fCPVCounter[i][j][k] = 0;
1033             }
1034         }
1035     }
1036
1037
1038   Int_t nentries = fCPVCell.GetEntries();
1039  
1040   Int_t   mtrackno, ism, ixp, iyp;
1041   Float_t edep;
1042   for (i = 0; i < nentries; i++)
1043     {
1044       AliPMDcell* cpvcell = (AliPMDcell*)fCPVCell.UncheckedAt(i);
1045       
1046       mtrackno = cpvcell->GetTrackNumber();
1047       ism      = cpvcell->GetSMNumber();
1048       ixp      = cpvcell->GetX();
1049       iyp      = cpvcell->GetY();
1050       edep     = cpvcell->GetEdep();
1051       Int_t nn = fCPVCounter[ism][ixp][iyp];
1052       cpvTrack[ism][ixp][iyp][nn] = (Int_t) mtrackno;
1053       cpvEdep[ism][ixp][iyp][nn] = edep;
1054       fCPVCounter[ism][ixp][iyp]++;
1055     }
1056   
1057   Int_t iz, il;
1058   Int_t im, ix, iy;
1059   Int_t nn;
1060   for (im=0; im<fgkTotUM; im++)
1061     {
1062       for (ix=0; ix<fgkRow; ix++)
1063         {
1064           for (iy=0; iy<fgkCol; iy++)
1065             {
1066               nn = fCPVCounter[im][ix][iy];
1067               if (nn > 1)
1068                 {
1069                   // This block handles if a cell is fired
1070                   // many times by many tracks
1071                   status1  = new Int_t[nn];
1072                   status2  = new Int_t[nn];
1073                   trnarray = new Int_t[nn];
1074                   for (iz = 0; iz < nn; iz++)
1075                     {
1076                       status1[iz] = cpvTrack[im][ix][iy][iz];
1077                     }
1078                   TMath::Sort(nn,status1,status2,jsort);
1079                   Int_t trackOld = -99999;
1080                   Int_t track, trCount = 0;
1081                   for (iz = 0; iz < nn; iz++)
1082                     {
1083                       track = status1[status2[iz]];
1084                       if (trackOld != track)
1085                         {
1086                           trnarray[trCount] = track;
1087                           trCount++;
1088                         }                             
1089                       trackOld = track;
1090                     }
1091                   delete [] status1;
1092                   delete [] status2;
1093                   Float_t totEdp = 0.;
1094                   trEdp = new Float_t[trCount];
1095                   fracEdp = new Float_t[trCount];
1096                   for (il = 0; il < trCount; il++)
1097                     {
1098                       trEdp[il] = 0.;
1099                       track = trnarray[il];
1100                       for (iz = 0; iz < nn; iz++)
1101                         {
1102                           if (track == cpvTrack[im][ix][iy][iz])
1103                             {
1104                               trEdp[il] += cpvEdep[im][ix][iy][iz];
1105                             }
1106                         }
1107                       totEdp += trEdp[il];
1108                     }
1109                   Int_t ilOld = 0;
1110                   Float_t fracOld = 0.;
1111                   
1112                   for (il = 0; il < trCount; il++)
1113                     {
1114                       fracEdp[il] = trEdp[il]/totEdp;
1115                       if (fracOld < fracEdp[il])
1116                         {
1117                           fracOld = fracEdp[il];
1118                           ilOld = il;
1119                         }
1120                     }
1121                   fCPVTrackNo[im][ix][iy] = trnarray[ilOld];
1122                   delete [] fracEdp;
1123                   delete [] trEdp;
1124                   delete [] trnarray;
1125                 }
1126               else if (nn == 1)
1127                 {
1128                   // This only handles if a cell is fired
1129                   // by only one track
1130                   
1131                   fCPVTrackNo[im][ix][iy] = cpvTrack[im][ix][iy][0];
1132                   
1133                 }
1134               else if (nn ==0)
1135                 {
1136                   // This is if no cell is fired
1137                   fCPVTrackNo[im][ix][iy] = -999;
1138                 }
1139             } // end of iy
1140         } // end of ix
1141     } // end of im
1142   
1143   // Delete all the pointers
1144   
1145  for (i = 0; i < fgkTotUM; i++)
1146     {
1147       for (j = 0; j < fgkRow; j++)
1148         {
1149           for (k = 0; k < fgkCol; k++)
1150             {
1151               delete []cpvTrack[i][j][k];
1152               delete []cpvEdep[i][j][k];
1153             }
1154         }
1155     }
1156  
1157   for (i = 0; i < fgkTotUM; i++)
1158     {
1159       for (j = 0; j < fgkRow; j++)
1160         {
1161           delete [] cpvTrack[i][j];
1162           delete [] cpvEdep[i][j];
1163         }
1164     }
1165   
1166   for (i = 0; i < fgkTotUM; i++)
1167     {
1168       delete [] cpvTrack[i];
1169       delete [] cpvEdep[i];
1170     }
1171   delete [] cpvTrack;
1172   delete [] cpvEdep;
1173   
1174   // 
1175   // End of the cell id assignment
1176   //
1177 }
1178 //____________________________________________________________________________
1179
1180 void AliPMDDigitizer::MergeSDigits(Int_t filenumber, Int_t troffset)
1181 {
1182   // merging sdigits
1183   fRunLoader = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(filenumber));
1184   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
1185   fPMDLoader->LoadSDigits("read");
1186   TTree* treeS = fPMDLoader->TreeS();
1187   AliPMDsdigit  *pmdsdigit;
1188   TBranch *branch = treeS->GetBranch("PMDSDigit");
1189   if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
1190   branch->SetAddress(&fSDigits);
1191
1192   Int_t   itrackno, itrackpid, idet, ism;
1193   Int_t   ixp, iyp;
1194   Float_t edep;
1195   Int_t nmodules = (Int_t) treeS->GetEntries();
1196   AliDebug(1,Form("Number of Modules in the treeS = %d",nmodules));
1197   AliDebug(1,Form("Track Offset = %d",troffset));
1198   for (Int_t imodule = 0; imodule < nmodules; imodule++)
1199     {
1200       treeS->GetEntry(imodule);
1201       Int_t nentries = fSDigits->GetLast();
1202       AliDebug(2,Form("Number of Entries per Module = %d",nentries));
1203       for (Int_t ient = 0; ient < nentries+1; ient++)
1204         {
1205           pmdsdigit = (AliPMDsdigit*)fSDigits->UncheckedAt(ient);
1206           itrackno  = pmdsdigit->GetTrackNumber();
1207           itrackpid = pmdsdigit->GetTrackPid();
1208           idet      = pmdsdigit->GetDetector();
1209           ism       = pmdsdigit->GetSMNumber();
1210           ixp       = pmdsdigit->GetRow();
1211           iyp       = pmdsdigit->GetColumn();
1212           edep      = pmdsdigit->GetCellEdep();
1213           if (idet == 0)
1214             {
1215               if (fPRE[ism][ixp][iyp] < edep)
1216                 {
1217                   fPRETrackNo[ism][ixp][iyp] = troffset + itrackno;
1218                   fPRETrackPid[ism][ixp][iyp] = itrackpid;
1219                 }
1220               fPRE[ism][ixp][iyp] += edep;
1221             }
1222           else if (idet == 1)
1223             {
1224               if (fCPV[ism][ixp][iyp] < edep)
1225                 {
1226                   fCPVTrackNo[ism][ixp][iyp] = troffset + itrackno;
1227                   fCPVTrackPid[ism][ixp][iyp] = itrackpid;
1228                 }
1229               fCPV[ism][ixp][iyp] += edep;
1230             }
1231         }
1232     }
1233
1234 }
1235 // ----------------------------------------------------------------------
1236 void AliPMDDigitizer::TrackAssignment2Cell()
1237 {
1238   // 
1239   // This block assigns the cell id when there are
1240   // multiple tracks in a cell according to the
1241   // energy deposition
1242   //
1243   Bool_t jsort = false;
1244
1245   Int_t i, j, k;
1246
1247   Int_t   *status1;
1248   Int_t   *status2;
1249   Int_t   *trnarray;
1250   Float_t *fracEdp;
1251   Float_t *trEdp;
1252   
1253   Int_t   ****pmdTrack;
1254   Float_t ****pmdEdep;
1255
1256   pmdTrack = new Int_t ***[fgkTotUM];
1257   pmdEdep  = new Float_t ***[fgkTotUM];
1258   for (i=0; i<fgkTotUM; i++)
1259     {
1260       pmdTrack[i] = new Int_t **[fgkRow];
1261       pmdEdep[i]  = new Float_t **[fgkRow];
1262     }
1263
1264   for (i = 0; i < fgkTotUM; i++)
1265     {
1266       for (j = 0; j < fgkRow; j++)
1267         {
1268           pmdTrack[i][j] = new Int_t *[fgkCol];
1269           pmdEdep[i][j]  = new Float_t *[fgkCol];
1270         }
1271     }
1272   
1273   for (i = 0; i < fgkTotUM; i++)
1274     {
1275       for (j = 0; j < fgkRow; j++)
1276         {
1277           for (k = 0; k < fgkCol; k++)
1278             {
1279               Int_t nn = fPRECounter[i][j][k];
1280               if(nn > 0)
1281                 {
1282                   pmdTrack[i][j][k] = new Int_t[nn];
1283                   pmdEdep[i][j][k] = new Float_t[nn];
1284                 }
1285               else
1286                 {
1287                   nn = 1;
1288                   pmdTrack[i][j][k] = new Int_t[nn];
1289                   pmdEdep[i][j][k] = new Float_t[nn];
1290                 }
1291               fPRECounter[i][j][k] = 0;
1292             }
1293         }
1294     }
1295
1296
1297   Int_t nentries = fCell.GetEntries();
1298
1299   Int_t   mtrackno, ism, ixp, iyp;
1300   Float_t edep;
1301
1302   for (i = 0; i < nentries; i++)
1303     {
1304       AliPMDcell* cell = (AliPMDcell*)fCell.UncheckedAt(i);
1305       
1306       mtrackno = cell->GetTrackNumber();
1307       ism      = cell->GetSMNumber();
1308       ixp      = cell->GetX();
1309       iyp      = cell->GetY();
1310       edep     = cell->GetEdep();
1311       Int_t nn = fPRECounter[ism][ixp][iyp];
1312       pmdTrack[ism][ixp][iyp][nn] = (Int_t) mtrackno;
1313       pmdEdep[ism][ixp][iyp][nn] = edep;
1314       fPRECounter[ism][ixp][iyp]++;
1315     }
1316   
1317   Int_t iz, il;
1318   Int_t im, ix, iy;
1319   Int_t nn;
1320   
1321   for (im=0; im<fgkTotUM; im++)
1322     {
1323       for (ix=0; ix<fgkRow; ix++)
1324         {
1325           for (iy=0; iy<fgkCol; iy++)
1326             {
1327               nn = fPRECounter[im][ix][iy];
1328               if (nn > 1)
1329                 {
1330                   // This block handles if a cell is fired
1331                   // many times by many tracks
1332                   status1  = new Int_t[nn];
1333                   status2  = new Int_t[nn];
1334                   trnarray = new Int_t[nn];
1335                   for (iz = 0; iz < nn; iz++)
1336                     {
1337                       status1[iz] = pmdTrack[im][ix][iy][iz];
1338                     }
1339                   TMath::Sort(nn,status1,status2,jsort);
1340                   Int_t trackOld = -99999;
1341                   Int_t track, trCount = 0;
1342                   for (iz = 0; iz < nn; iz++)
1343                     {
1344                       track = status1[status2[iz]];
1345                       if (trackOld != track)
1346                         {
1347                           trnarray[trCount] = track;
1348                           trCount++;
1349                         }
1350                       trackOld = track;
1351                     }
1352                   delete [] status1;
1353                   delete [] status2;
1354                   Float_t totEdp = 0.;
1355                   trEdp = new Float_t[trCount];
1356                   fracEdp = new Float_t[trCount];
1357                   for (il = 0; il < trCount; il++)
1358                     {
1359                       trEdp[il] = 0.;
1360                       track = trnarray[il];
1361                       for (iz = 0; iz < nn; iz++)
1362                         {
1363                           if (track == pmdTrack[im][ix][iy][iz])
1364                             {
1365                               trEdp[il] += pmdEdep[im][ix][iy][iz];
1366                             }
1367                         }
1368                       totEdp += trEdp[il];
1369                     }
1370                   Int_t ilOld = 0;
1371                   Float_t fracOld = 0.;
1372                   
1373                   for (il = 0; il < trCount; il++)
1374                     {
1375                       fracEdp[il] = trEdp[il]/totEdp;
1376                       if (fracOld < fracEdp[il])
1377                         {
1378                           fracOld = fracEdp[il];
1379                           ilOld = il;
1380                         }
1381                     }
1382                   fPRETrackNo[im][ix][iy] = trnarray[ilOld];
1383                   delete [] fracEdp;
1384                   delete [] trEdp;
1385                   delete [] trnarray;
1386                 }
1387               else if (nn == 1)
1388                 {
1389                   // This only handles if a cell is fired
1390                   // by only one track
1391                   
1392                   fPRETrackNo[im][ix][iy] = pmdTrack[im][ix][iy][0];
1393                   
1394                 }
1395               else if (nn ==0)
1396                 {
1397                   // This is if no cell is fired
1398                   fPRETrackNo[im][ix][iy] = -999;
1399                 }
1400             } // end of iy
1401         } // end of ix
1402     } // end of im
1403   
1404   // Delete all the pointers
1405   
1406   for (i = 0; i < fgkTotUM; i++)
1407     {
1408       for (j = 0; j < fgkRow; j++)
1409         {
1410           for (k = 0; k < fgkCol; k++)
1411             {
1412               delete [] pmdTrack[i][j][k];
1413               delete [] pmdEdep[i][j][k];
1414             }
1415         }
1416     }
1417   
1418   for (i = 0; i < fgkTotUM; i++)
1419     {
1420       for (j = 0; j < fgkRow; j++)
1421         {
1422           delete [] pmdTrack[i][j];
1423           delete [] pmdEdep[i][j];
1424         }
1425     }
1426   
1427   for (i = 0; i < fgkTotUM; i++)
1428     {
1429       delete [] pmdTrack[i];
1430       delete [] pmdEdep[i];
1431     }
1432   delete [] pmdTrack;
1433   delete [] pmdEdep;
1434   // 
1435   // End of the cell id assignment
1436   //
1437 }
1438 //____________________________________________________________________________
1439 void AliPMDDigitizer::MeV2ADC(Float_t mev, Float_t & adc) const
1440 {
1441   // This converts the simulated edep to ADC according to the
1442   // Test Beam Data
1443   //PS Test in May 2009, Voltage @ 1350 V
1444   // KeV - ADC conversion for 12bit ADC
1445   // MPV data used for the fit and taken here
1446
1447   const Float_t kConstant   = -0.1602;
1448   const Float_t kErConstant =  0.9914;
1449   const Float_t kSlope      = 77.47;
1450   const Float_t kErSlope    =  3.16;
1451   
1452   Float_t cons   = gRandom->Gaus(kConstant,kErConstant);
1453   Float_t slop   = gRandom->Gaus(kSlope,kErSlope);
1454   
1455   Float_t adc12bit = slop*mev*0.001 + cons;
1456   
1457   if (adc12bit < 0.) adc12bit = 0.;
1458                            
1459   if(adc12bit < 1600.0)
1460     {
1461       adc = (Float_t) adc12bit;
1462     }
1463   else if (adc12bit >= 1600.0)
1464     {
1465       adc = 1600.0;
1466     }
1467 }
1468 //____________________________________________________________________________
1469 void AliPMDDigitizer::AddSDigit(Int_t trnumber, Int_t trpid, Int_t det,
1470                                 Int_t smnumber, Int_t irow, Int_t icol,
1471                                 Float_t adc)
1472 {
1473   // Add SDigit
1474   //
1475   if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
1476   TClonesArray &lsdigits = *fSDigits;
1477   new(lsdigits[fNsdigit++])  AliPMDsdigit(trnumber,trpid,det,smnumber,irow,icol,adc);
1478 }
1479 //____________________________________________________________________________
1480
1481 void AliPMDDigitizer::AddDigit(Int_t trnumber, Int_t trpid, Int_t det,
1482                                Int_t smnumber, Int_t irow, Int_t icol,
1483                                Float_t adc)
1484 {
1485   // Add Digit
1486   //
1487   if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
1488   TClonesArray &ldigits = *fDigits;
1489   new(ldigits[fNdigit++]) AliPMDdigit(trnumber,trpid, det,smnumber,irow,icol,adc);
1490 }
1491 //____________________________________________________________________________
1492
1493 void AliPMDDigitizer::SetZPosition(Float_t zpos)
1494 {
1495   fZPos = zpos;
1496 }
1497 //____________________________________________________________________________
1498 Float_t AliPMDDigitizer::GetZPosition() const
1499 {
1500   return fZPos;
1501 }
1502 //____________________________________________________________________________
1503
1504 void AliPMDDigitizer::ResetCell()
1505 {
1506   // clears the cell array and also the counter
1507   //  for each cell
1508   //
1509   fCPVCell.Delete();
1510   fCell.Delete();
1511   for (Int_t i = 0; i < fgkTotUM; i++)
1512     {
1513       for (Int_t j = 0; j < fgkRow; j++)
1514         {
1515           for (Int_t k = 0; k < fgkCol; k++)
1516             {
1517               fCPVCounter[i][j][k] = 0; 
1518               fPRECounter[i][j][k] = 0;
1519             }
1520         }
1521     }
1522 }
1523 //____________________________________________________________________________
1524 void AliPMDDigitizer::ResetSDigit()
1525 {
1526   // Clears SDigits
1527   fNsdigit = 0;
1528   if (fSDigits) fSDigits->Delete();
1529 }
1530 //____________________________________________________________________________
1531 void AliPMDDigitizer::ResetDigit()
1532 {
1533   // Clears Digits
1534   fNdigit = 0;
1535   if (fDigits) fDigits->Delete();
1536 }
1537 //____________________________________________________________________________
1538
1539 void AliPMDDigitizer::ResetCellADC()
1540 {
1541   // Clears individual cells edep and track number
1542   for (Int_t i = 0; i < fgkTotUM; i++)
1543     {
1544       for (Int_t j = 0; j < fgkRow; j++)
1545         {
1546           for (Int_t k = 0; k < fgkCol; k++)
1547             {
1548               fCPV[i][j][k]         = 0.;
1549               fPRE[i][j][k]         = 0.;
1550               fCPVTrackNo[i][j][k]  = 0;
1551               fPRETrackNo[i][j][k]  = 0;
1552               fCPVTrackPid[i][j][k] = -1;
1553               fPRETrackPid[i][j][k] = -1;
1554             }
1555         }
1556     }
1557 }
1558 //____________________________________________________________________________
1559
1560 void AliPMDDigitizer::UnLoad(Option_t *option)
1561 {
1562   // Unloads all the root files
1563   //
1564   const char *cS = strstr(option,"S");
1565   const char *cD = strstr(option,"D");
1566
1567   fRunLoader->UnloadgAlice();
1568   fRunLoader->UnloadHeader();
1569   fRunLoader->UnloadKinematics();
1570
1571   if (cS)
1572     {
1573       fPMDLoader->UnloadHits();
1574     }
1575   if (cD)
1576     {
1577       fPMDLoader->UnloadHits();
1578       fPMDLoader->UnloadSDigits();
1579     }
1580 }
1581
1582 //----------------------------------------------------------------------
1583 Float_t AliPMDDigitizer::Gain(Int_t det, Int_t smn, Int_t row, Int_t col) const
1584 {
1585   // returns of the gain of the cell
1586   // Added this method by ZA
1587
1588   //cout<<" I am here in gain "<<fCalibData<< "smn,row, col "<<smn
1589   //<<" "<<row<<" "<<col<<endl;
1590
1591   if(!fCalibGain) {
1592     AliError("No calibration data loaded from CDB!!!");
1593     return 1;
1594   }
1595
1596   Float_t GainFact;
1597   GainFact = fCalibGain->GetGainFact(det,smn,row,col);
1598   return GainFact;
1599 }
1600 //----------------------------------------------------------------------
1601 AliPMDCalibData* AliPMDDigitizer::GetCalibGain() const
1602 {
1603   // The run number will be centralized in AliCDBManager,
1604   // you don't need to set it here!
1605   // Added this method by ZA
1606   // Cleaned up by Alberto
1607   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("PMD/Calib/Gain");
1608   
1609   if(!entry) AliFatal("Calibration object retrieval failed!");
1610   
1611   AliPMDCalibData *calibdata=0;
1612   if (entry) calibdata = (AliPMDCalibData*) entry->GetObject();
1613   
1614   if (!calibdata)  AliFatal("No calibration data from calibration database !");
1615   
1616   return calibdata;
1617 }
1618 //----------------------------------------------------------------------
1619 AliPMDPedestal* AliPMDDigitizer::GetCalibPed() const
1620 {
1621   // The run number will be centralized in AliCDBManager,
1622   // you don't need to set it here!
1623
1624   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("PMD/Calib/Ped");
1625   
1626   if(!entry) AliFatal("Pedestal object retrieval failed!");
1627   
1628   AliPMDPedestal *pedestal=0;
1629   if (entry) pedestal = (AliPMDPedestal*) entry->GetObject();
1630   
1631   if (!pedestal)  AliFatal("No pedestal data from calibration database !");
1632   
1633   return pedestal;
1634 }