]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PMD/AliPMDDigitizer.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Try to get correct calibration object - in meant time the TRD track written to the...
[u/mrichter/AliRoot.git] / PMD / AliPMDDigitizer.cxx
... / ...
CommitLineData
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
61ClassImp(AliPMDDigitizer)
62
63AliPMDDigitizer::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//____________________________________________________________________________
97AliPMDDigitizer::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//____________________________________________________________________________
118AliPMDDigitizer & AliPMDDigitizer::operator=(const AliPMDDigitizer& /*digitizer*/)
119{
120 // Assignment operator
121 AliError("Assignement operator not allowed ");
122
123 return *this;
124}
125//____________________________________________________________________________
126AliPMDDigitizer::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//____________________________________________________________________________
162AliPMDDigitizer::~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//____________________________________________________________________________
182void 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//____________________________________________________________________________
246void 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
480void 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
754void 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//____________________________________________________________________________
840void 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
947void 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// ----------------------------------------------------------------------
1000void 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//____________________________________________________________________________
1202void 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//____________________________________________________________________________
1239void 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
1250void 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
1261void AliPMDDigitizer::SetZPosition(Float_t zpos)
1262{
1263 fZPos = zpos;
1264}
1265//____________________________________________________________________________
1266Float_t AliPMDDigitizer::GetZPosition() const
1267{
1268 return fZPos;
1269}
1270//____________________________________________________________________________
1271
1272void 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//____________________________________________________________________________
1290void AliPMDDigitizer::ResetSDigit()
1291{
1292 // Clears SDigits
1293 fNsdigit = 0;
1294 if (fSDigits) fSDigits->Delete();
1295}
1296//____________________________________________________________________________
1297void AliPMDDigitizer::ResetDigit()
1298{
1299 // Clears Digits
1300 fNdigit = 0;
1301 if (fDigits) fDigits->Delete();
1302}
1303//____________________________________________________________________________
1304
1305void 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
1324void 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//----------------------------------------------------------------------
1347Float_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//----------------------------------------------------------------------
1365AliPMDCalibData* 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//----------------------------------------------------------------------
1383AliPMDPedestal* 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}