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