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