]> 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
Removed CDB in /home/basanta/ALISOFT/AliRoot
[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 "AliPMDDigitizer.h"
58
59
60ClassImp(AliPMDDigitizer)
61
62AliPMDDigitizer::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//____________________________________________________________________________
96AliPMDDigitizer::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//____________________________________________________________________________
116AliPMDDigitizer & AliPMDDigitizer::operator=(const AliPMDDigitizer& /*digitizer*/)
117{
118 // Assignment operator
119 AliError("Assignement operator not allowed ");
120
121 return *this;
122}
123//____________________________________________________________________________
124AliPMDDigitizer::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//____________________________________________________________________________
159AliPMDDigitizer::~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//____________________________________________________________________________
179void 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//____________________________________________________________________________
243void 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
477void 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
738void 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//____________________________________________________________________________
813void 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
906void 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// ----------------------------------------------------------------------
959void 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//____________________________________________________________________________
1161void 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//____________________________________________________________________________
1198void 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
1209void 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
1220void AliPMDDigitizer::SetZPosition(Float_t zpos)
1221{
1222 fZPos = zpos;
1223}
1224//____________________________________________________________________________
1225Float_t AliPMDDigitizer::GetZPosition() const
1226{
1227 return fZPos;
1228}
1229//____________________________________________________________________________
1230
1231void 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//____________________________________________________________________________
1249void AliPMDDigitizer::ResetSDigit()
1250{
1251 // Clears SDigits
1252 fNsdigit = 0;
1253 if (fSDigits) fSDigits->Delete();
1254}
1255//____________________________________________________________________________
1256void AliPMDDigitizer::ResetDigit()
1257{
1258 // Clears Digits
1259 fNdigit = 0;
1260 if (fDigits) fDigits->Delete();
1261}
1262//____________________________________________________________________________
1263
1264void 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
1283void 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//----------------------------------------------------------------------
1306Float_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//----------------------------------------------------------------------
1324AliPMDCalibData* 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}