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