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