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