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