Coding conventions obeyed
[u/mrichter/AliRoot.git] / PMD / AliPMDDigitizer.cxx
CommitLineData
b9746a9e 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 **************************************************************************/
b2b13ee2 15//-----------------------------------------------------//
16// //
8599803d 17// Source File : PMDDigitizer.cxx, Version 00 //
b2b13ee2 18// //
19// Date : September 20 2002 //
20// //
21//-----------------------------------------------------//
22
36031625 23#include <Riostream.h>
b2b13ee2 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 "AliRun.h"
35#include "AliPMD.h"
36#include "AliHit.h"
37#include "AliDetector.h"
38#include "AliRunLoader.h"
39#include "AliLoader.h"
40#include "AliConfig.h"
41#include "AliMagF.h"
42#include "AliRunDigitizer.h"
43#include "AliHeader.h"
44
45#include "AliPMDcell.h"
46#include "AliPMDsdigit.h"
47#include "AliPMDdigit.h"
48#include "AliPMDDigitizer.h"
49#include "AliPMDClustering.h"
5d12ce38 50#include "AliMC.h"
b2b13ee2 51
52ClassImp(AliPMDDigitizer)
53//
54// Constructor
55//
56AliPMDDigitizer::AliPMDDigitizer()
57{
36031625 58 // Default Constructor
59 //
b2b13ee2 60 if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
61 fNsdigit = 0;
62 if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
63 fNdigit = 0;
64
36031625 65 for (Int_t i = 0; i < fgkTotUM; i++)
b2b13ee2 66 {
36031625 67 for (Int_t j = 0; j < fgkRow; j++)
b2b13ee2 68 {
36031625 69 for (Int_t k = 0; k < fgkCol; k++)
b2b13ee2 70 {
71 fCPV[i][j][k] = 0.;
36031625 72 fPRE[i][j][k] = 0.;
b2b13ee2 73 }
74 }
75 }
76
77 if (!fCell) fCell = new TObjArray();
78
79 fZPos = 361.5; // in units of cm, This is the default position of PMD
b2b13ee2 80}
81AliPMDDigitizer::~AliPMDDigitizer()
82{
36031625 83 // Default Destructor
84 //
b2b13ee2 85 delete fSDigits;
86 delete fDigits;
87 delete fCell;
88}
89//
90// Member functions
91//
92void AliPMDDigitizer::OpengAliceFile(Char_t *file, Option_t *option)
93{
36031625 94 // Loads galice.root file and corresponding header, kinematics
95 // hits and sdigits or digits depending on the option
96 //
b2b13ee2 97 fRunLoader = AliRunLoader::Open(file,AliConfig::fgkDefaultEventFolderName,
98 "UPDATE");
99
100 if (!fRunLoader)
101 {
102 Error("Open","Can not open session for file %s.",file);
103 }
104
105 fRunLoader->LoadgAlice();
106 fRunLoader->LoadHeader();
107 fRunLoader->LoadKinematics();
108
36031625 109 fAlice = fRunLoader->GetAliRun();
b2b13ee2 110
36031625 111 if (fAlice)
b2b13ee2 112 {
113 printf("<AliPMDdigitizer::Open> ");
114 printf("AliRun object found on file.\n");
115 }
116 else
117 {
118 printf("<AliPMDdigitizer::Open> ");
119 printf("Could not find AliRun object.\n");
120 }
121
36031625 122 fPMD = (AliPMD*)fAlice->GetDetector("PMD");
123 fPMDLoader = fRunLoader->GetLoader("PMDLoader");
124 if (fPMDLoader == 0x0)
b2b13ee2 125 {
126 cerr<<"Hits2Digits : Can not find PMD or PMDLoader\n";
127 }
128
129 const char *cHS = strstr(option,"HS");
130 const char *cHD = strstr(option,"HD");
131 const char *cSD = strstr(option,"SD");
132
133 if (cHS)
134 {
36031625 135 fPMDLoader->LoadHits("READ");
136 fPMDLoader->LoadSDigits("recreate");
b2b13ee2 137 }
138 else if (cHD)
139 {
36031625 140 fPMDLoader->LoadHits("READ");
141 fPMDLoader->LoadDigits("recreate");
b2b13ee2 142 }
143 else if (cSD)
144 {
36031625 145 fPMDLoader->LoadSDigits("READ");
146 fPMDLoader->LoadDigits("recreate");
b2b13ee2 147 }
148
149}
150void AliPMDDigitizer::Hits2SDigits(Int_t ievt)
151{
36031625 152 // This reads the PMD Hits tree and assigns the right track number
153 // to a cell and stores in the summable digits tree
154 //
155 // cout << " -------- Beginning of Hits2SDigits ----------- " << endl;
b2b13ee2 156
36031625 157 const Int_t kPi0 = 111;
158 const Int_t kGamma = 22;
b2b13ee2 159 Int_t npmd;
160 Int_t trackno;
b2b13ee2 161 Int_t smnumber;
162 Int_t trackpid;
163 Int_t mtrackno;
164 Int_t mtrackpid;
165
166 Float_t xPos, yPos, zPos;
8599803d 167 Int_t xpad = -1, ypad = -1;
b2b13ee2 168 Float_t edep;
169 Float_t vx = -999.0, vy = -999.0, vz = -999.0;
170
171
172 ResetSDigit();
173
174 printf("Event Number = %d \n",ievt);
175 Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
176 printf("Number of Particles = %d \n", nparticles);
177 fRunLoader->GetEvent(ievt);
36031625 178 fPArray = fAlice->GetMCApp()->Particles();
b2b13ee2 179 // ------------------------------------------------------- //
180 // Pointer to specific detector hits.
181 // Get pointers to Alice detectors and Hits containers
182
36031625 183 fTreeH = fPMDLoader->TreeH();
b2b13ee2 184
36031625 185 Int_t ntracks = (Int_t) fTreeH->GetEntries();
b2b13ee2 186 printf("Number of Tracks in the TreeH = %d \n", ntracks);
187
36031625 188 fTreeS = fPMDLoader->TreeS();
189 if (fTreeS == 0x0)
b2b13ee2 190 {
36031625 191 fPMDLoader->MakeTree("S");
192 fTreeS = fPMDLoader->TreeS();
b2b13ee2 193 }
194 Int_t bufsize = 16000;
36031625 195 fTreeS->Branch("PMDSDigit", &fSDigits, bufsize);
b2b13ee2 196
36031625 197 if (fPMD) fHits = fPMD->Hits();
b2b13ee2 198
199 // Start loop on tracks in the hits containers
200
201
202 for (Int_t track=0; track<ntracks;track++)
203 {
36031625 204 fAlice->ResetHits();
205 fTreeH->GetEvent(track);
b2b13ee2 206
36031625 207 if (fPMD)
b2b13ee2 208 {
36031625 209 npmd = fHits->GetEntriesFast();
b2b13ee2 210 for (int ipmd = 0; ipmd < npmd; ipmd++)
211 {
36031625 212 fPMDHit = (AliPMDhit*) fHits->UncheckedAt(ipmd);
213 trackno = fPMDHit->GetTrack();
b2b13ee2 214
215 // get kinematics of the particles
36031625 216
217 fParticle = fAlice->GetMCApp()->Particle(trackno);
218 trackpid = fParticle->GetPdgCode();
b2b13ee2 219
220 Int_t igatr = -999;
221 Int_t ichtr = -999;
8599803d 222 Int_t igapid = -999;
b2b13ee2 223 Int_t imo;
224 Int_t igen = 0;
36031625 225 Int_t idmo = -999;
8599803d 226
36031625 227 TParticle* mparticle = fParticle;
228 Int_t tracknoOld=0, trackpidOld=0, statusOld = 0;
8599803d 229 if (mparticle->GetFirstMother() == -1)
230 {
36031625 231 tracknoOld = trackno;
232 trackpidOld = trackpid;
233 statusOld = -1;
8599803d 234 }
235
236 Int_t igstatus = 0;
b2b13ee2 237 while((imo = mparticle->GetFirstMother()) >= 0)
238 {
239 igen++;
36031625 240
241 mparticle = fAlice->GetMCApp()->Particle(imo);
242 idmo = mparticle->GetPdgCode();
b2b13ee2 243
244 vx = mparticle->Vx();
245 vy = mparticle->Vy();
246 vz = mparticle->Vz();
247
36031625 248 //printf("==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, idmo, vx, vy, vz);
249 //fprintf(fpw1,"==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, idmo, vx, vy, vz);
250 if ((idmo == kGamma || idmo == -11 || idmo == 11) && vx == 0. && vy == 0. && vz == 0.)
b2b13ee2 251 {
252 igatr = imo;
36031625 253 igapid = idmo;
8599803d 254 igstatus = 1;
255 }
256 if(igstatus == 0)
257 {
36031625 258 if (idmo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
8599803d 259 {
260 igatr = imo;
36031625 261 igapid = idmo;
8599803d 262 }
b2b13ee2 263 }
264 ichtr = imo;
265 }
266
36031625 267 if (idmo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
b2b13ee2 268 {
269 mtrackno = igatr;
8599803d 270 mtrackpid = igapid;
b2b13ee2 271 }
272 else
273 {
274 mtrackno = ichtr;
36031625 275 mtrackpid = idmo;
b2b13ee2 276 }
36031625 277 if (statusOld == -1)
8599803d 278 {
36031625 279 mtrackno = tracknoOld;
280 mtrackpid = trackpidOld;
8599803d 281 }
282
36031625 283 xPos = fPMDHit->X();
284 yPos = fPMDHit->Y();
285 zPos = fPMDHit->Z();
286 edep = fPMDHit->fEnergy;
287 Int_t vol1 = fPMDHit->fVolume[1]; // Column
288 Int_t vol2 = fPMDHit->fVolume[2]; // Row
289 Int_t vol3 = fPMDHit->fVolume[3]; // UnitModule
290 Int_t vol6 = fPMDHit->fVolume[6]; // SuperModule
8599803d 291 // -----------------------------------------//
292 // For Super Module 1 & 2 //
293 // nrow = 96, ncol = 48 //
294 // For Super Module 3 & 4 //
295 // nrow = 48, ncol = 96 //
296 // -----------------------------------------//
b2b13ee2 297
36031625 298 smnumber = (vol6-1)*6 + vol3;
8599803d 299
36031625 300 if (vol6 == 1 || vol6 == 2)
8599803d 301 {
36031625 302 xpad = vol1;
303 ypad = vol2;
8599803d 304 }
36031625 305 else if (vol6 == 3 || vol6 == 4)
b2b13ee2 306 {
36031625 307 xpad = vol2;
308 ypad = vol1;
b2b13ee2 309 }
b2b13ee2 310
75bb524e 311 //cout << "zpos = " << zPos << " edep = " << edep << endl;
b2b13ee2 312
8599803d 313 Float_t zposition = TMath::Abs(zPos);
314 if (zposition < fZPos)
b2b13ee2 315 {
316 // CPV
317 fDetNo = 1;
318 }
8599803d 319 else if (zposition > fZPos)
b2b13ee2 320 {
321 // PMD
322 fDetNo = 0;
323 }
8599803d 324 Int_t smn = smnumber - 1;
325 Int_t ixx = xpad - 1;
326 Int_t iyy = ypad - 1;
b2b13ee2 327 if (fDetNo == 0)
328 {
36031625 329 fPRE[smn][ixx][iyy] += edep;
330 fPRECounter[smn][ixx][iyy]++;
b2b13ee2 331
36031625 332 fPMDcell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
b2b13ee2 333
36031625 334 fCell->Add(fPMDcell);
b2b13ee2 335 }
336 else if(fDetNo == 1)
337 {
8599803d 338 fCPV[smn][ixx][iyy] += edep;
339 fCPVTrackNo[smn][ixx][iyy] = mtrackno;
b2b13ee2 340 }
341 }
342 }
343 } // Track Loop ended
8599803d 344
b2b13ee2 345 TrackAssignment2Cell();
b2b13ee2 346 ResetCell();
347
348 Float_t deltaE = 0.;
349 Int_t detno = 0;
b2b13ee2 350 Int_t trno = -1;
351 Int_t cellno = 0;
352
353 for (Int_t idet = 0; idet < 2; idet++)
354 {
36031625 355 for (Int_t ism = 0; ism < fgkTotUM; ism++)
b2b13ee2 356 {
36031625 357 for (Int_t jrow = 0; jrow < fgkRow; jrow++)
b2b13ee2 358 {
36031625 359 for (Int_t kcol = 0; kcol < fgkCol; kcol++)
b2b13ee2 360 {
36031625 361 cellno = jrow*fgkCol + kcol;
b2b13ee2 362 if (idet == 0)
363 {
36031625 364 deltaE = fPRE[ism][jrow][kcol];
365 trno = fPRETrackNo[ism][jrow][kcol];
b2b13ee2 366 detno = 0;
367 }
368 else if (idet == 1)
369 {
370 deltaE = fCPV[ism][jrow][kcol];
371 trno = fCPVTrackNo[ism][jrow][kcol];
372 detno = 1;
373 }
374 if (deltaE > 0.)
375 {
b2b13ee2 376 AddSDigit(trno,detno,ism,cellno,deltaE);
377 }
378 }
379 }
36031625 380 fTreeS->Fill();
b2b13ee2 381 ResetSDigit();
382 }
383 }
36031625 384 fPMDLoader->WriteSDigits("OVERWRITE");
b2b13ee2 385 ResetCellADC();
386
387 // cout << " -------- End of Hits2SDigit ----------- " << endl;
388}
389
390void AliPMDDigitizer::Hits2Digits(Int_t ievt)
391{
36031625 392 // This reads the PMD Hits tree and assigns the right track number
393 // to a cell and stores in the digits tree
394 //
395 const Int_t kPi0 = 111;
396 const Int_t kGamma = 22;
b2b13ee2 397 Int_t npmd;
398 Int_t trackno;
b2b13ee2 399 Int_t smnumber;
400 Int_t trackpid;
401 Int_t mtrackno;
402 Int_t mtrackpid;
403
404 Float_t xPos, yPos, zPos;
8599803d 405 Int_t xpad = -1, ypad = -1;
b2b13ee2 406 Float_t edep;
407 Float_t vx = -999.0, vy = -999.0, vz = -999.0;
408
409
410 ResetDigit();
411
412 printf("Event Number = %d \n",ievt);
413
414 Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
415 printf("Number of Particles = %d \n", nparticles);
416 fRunLoader->GetEvent(ievt);
36031625 417 fPArray = fAlice->GetMCApp()->Particles();
b2b13ee2 418 // ------------------------------------------------------- //
419 // Pointer to specific detector hits.
420 // Get pointers to Alice detectors and Hits containers
421
36031625 422 fPMD = (AliPMD*)fAlice->GetDetector("PMD");
423 fPMDLoader = fRunLoader->GetLoader("PMDLoader");
b2b13ee2 424
36031625 425 if (fPMDLoader == 0x0)
b2b13ee2 426 {
427 cerr<<"Hits2Digits method : Can not find PMD or PMDLoader\n";
428 }
36031625 429 fTreeH = fPMDLoader->TreeH();
430 Int_t ntracks = (Int_t) fTreeH->GetEntries();
b2b13ee2 431 printf("Number of Tracks in the TreeH = %d \n", ntracks);
36031625 432 fPMDLoader->LoadDigits("recreate");
433 fTreeD = fPMDLoader->TreeD();
434 if (fTreeD == 0x0)
b2b13ee2 435 {
36031625 436 fPMDLoader->MakeTree("D");
437 fTreeD = fPMDLoader->TreeD();
b2b13ee2 438 }
439 Int_t bufsize = 16000;
36031625 440 fTreeD->Branch("PMDDigit", &fDigits, bufsize);
b2b13ee2 441
36031625 442 if (fPMD) fHits = fPMD->Hits();
b2b13ee2 443
444 // Start loop on tracks in the hits containers
445
446 for (Int_t track=0; track<ntracks;track++)
447 {
36031625 448 fAlice->ResetHits();
449 fTreeH->GetEvent(track);
b2b13ee2 450
36031625 451 if (fPMD)
b2b13ee2 452 {
36031625 453 npmd = fHits->GetEntriesFast();
b2b13ee2 454 for (int ipmd = 0; ipmd < npmd; ipmd++)
455 {
36031625 456 fPMDHit = (AliPMDhit*) fHits->UncheckedAt(ipmd);
457 trackno = fPMDHit->GetTrack();
b2b13ee2 458
459 // get kinematics of the particles
460
36031625 461 fParticle = fAlice->GetMCApp()->Particle(trackno);
462 trackpid = fParticle->GetPdgCode();
b2b13ee2 463
464 Int_t igatr = -999;
465 Int_t ichtr = -999;
8599803d 466 Int_t igapid = -999;
b2b13ee2 467 Int_t imo;
468 Int_t igen = 0;
36031625 469 Int_t idmo = -999;
8599803d 470
36031625 471 TParticle* mparticle = fParticle;
472 Int_t tracknoOld=0, trackpidOld=0, statusOld = 0;
8599803d 473 if (mparticle->GetFirstMother() == -1)
474 {
36031625 475 tracknoOld = trackno;
476 trackpidOld = trackpid;
477 statusOld = -1;
8599803d 478 }
479
480 Int_t igstatus = 0;
b2b13ee2 481 while((imo = mparticle->GetFirstMother()) >= 0)
482 {
483 igen++;
36031625 484
485 mparticle = fAlice->GetMCApp()->Particle(imo);
486 idmo = mparticle->GetPdgCode();
b2b13ee2 487
488 vx = mparticle->Vx();
489 vy = mparticle->Vy();
490 vz = mparticle->Vz();
491
36031625 492 //printf("==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, idmo, vx, vy, vz);
493 //fprintf(fpw1,"==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, idmo, vx, vy, vz);
494 if ((idmo == kGamma || idmo == -11 || idmo == 11) && vx == 0. && vy == 0. && vz == 0.)
b2b13ee2 495 {
496 igatr = imo;
36031625 497 igapid = idmo;
8599803d 498 igstatus = 1;
499 }
500 if(igstatus == 0)
501 {
36031625 502 if (idmo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
8599803d 503 {
504 igatr = imo;
36031625 505 igapid = idmo;
8599803d 506 }
b2b13ee2 507 }
508 ichtr = imo;
509 }
510
36031625 511 if (idmo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
b2b13ee2 512 {
513 mtrackno = igatr;
8599803d 514 mtrackpid = igapid;
b2b13ee2 515 }
516 else
517 {
518 mtrackno = ichtr;
36031625 519 mtrackpid = idmo;
b2b13ee2 520 }
36031625 521 if (statusOld == -1)
8599803d 522 {
36031625 523 mtrackno = tracknoOld;
524 mtrackpid = trackpidOld;
8599803d 525 }
b2b13ee2 526
36031625 527 xPos = fPMDHit->X();
528 yPos = fPMDHit->Y();
529 zPos = fPMDHit->Z();
530 edep = fPMDHit->fEnergy;
531 Int_t vol1 = fPMDHit->fVolume[1]; // Column
532 Int_t vol2 = fPMDHit->fVolume[2]; // Row
533 Int_t vol3 = fPMDHit->fVolume[3]; // UnitModule
534 Int_t vol6 = fPMDHit->fVolume[6]; // SuperModule
8599803d 535
536 // -----------------------------------------//
537 // For Super Module 1 & 2 //
538 // nrow = 96, ncol = 48 //
539 // For Super Module 3 & 4 //
540 // nrow = 48, ncol = 96 //
541 // -----------------------------------------//
b2b13ee2 542
36031625 543 smnumber = (vol6-1)*6 + vol3;
8599803d 544
36031625 545 if (vol6 == 1 || vol6 == 2)
b2b13ee2 546 {
36031625 547 xpad = vol1;
548 ypad = vol2;
8599803d 549 }
36031625 550 else if (vol6 == 3 || vol6 == 4)
8599803d 551 {
36031625 552 xpad = vol2;
553 ypad = vol1;
b2b13ee2 554 }
b2b13ee2 555
556 //cout << "-zpos = " << -zPos << endl;
8599803d 557
558 Float_t zposition = TMath::Abs(zPos);
559
560 if (zposition < fZPos)
b2b13ee2 561 {
562 // CPV
563 fDetNo = 1;
564 }
8599803d 565 else if (zposition > fZPos)
b2b13ee2 566 {
567 // PMD
568 fDetNo = 0;
569 }
8599803d 570
571 Int_t smn = smnumber - 1;
572 Int_t ixx = xpad - 1;
573 Int_t iyy = ypad - 1;
574 if (fDetNo == 0)
b2b13ee2 575 {
36031625 576 fPRE[smn][ixx][iyy] += edep;
577 fPRECounter[smn][ixx][iyy]++;
8599803d 578
36031625 579 fPMDcell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
8599803d 580
36031625 581 fCell->Add(fPMDcell);
b2b13ee2 582 }
8599803d 583 else if(fDetNo == 1)
b2b13ee2 584 {
8599803d 585 fCPV[smn][ixx][iyy] += edep;
586 fCPVTrackNo[smn][ixx][iyy] = mtrackno;
b2b13ee2 587 }
588 }
589 }
590 } // Track Loop ended
591
592 TrackAssignment2Cell();
593 ResetCell();
594
595 Float_t deltaE = 0.;
596 Int_t detno = 0;
b2b13ee2 597 Int_t trno = 1;
598 Int_t cellno;
599
600 for (Int_t idet = 0; idet < 2; idet++)
601 {
36031625 602 for (Int_t ism = 0; ism < fgkTotUM; ism++)
b2b13ee2 603 {
36031625 604 for (Int_t jrow = 0; jrow < fgkRow; jrow++)
b2b13ee2 605 {
36031625 606 for (Int_t kcol = 0; kcol < fgkCol; kcol++)
b2b13ee2 607 {
36031625 608 cellno = jrow*fgkCol + kcol;
b2b13ee2 609 if (idet == 0)
610 {
36031625 611 deltaE = fPRE[ism][jrow][kcol];
612 trno = fPRETrackNo[ism][jrow][kcol];
b2b13ee2 613 detno = 0;
614 }
615 else if (idet == 1)
616 {
617 deltaE = fCPV[ism][jrow][kcol];
8599803d 618 trno = fCPVTrackNo[ism][jrow][kcol];
b2b13ee2 619 detno = 1;
620 }
621 if (deltaE > 0.)
622 {
b2b13ee2 623 AddDigit(trno,detno,ism,cellno,deltaE);
624 }
625 } // column loop
626 } // row loop
627 } // supermodule loop
36031625 628 fTreeD->Fill();
b2b13ee2 629 ResetDigit();
630 } // detector loop
631
36031625 632 fPMDLoader->WriteDigits("OVERWRITE");
b2b13ee2 633 ResetCellADC();
634
635 // cout << " -------- End of Hits2Digit ----------- " << endl;
636}
637
638
639void AliPMDDigitizer::SDigits2Digits(Int_t ievt)
640{
36031625 641 // This reads the PMD sdigits tree and converts energy deposition
642 // in a cell to ADC and stores in the digits tree
643 //
b2b13ee2 644 // cout << " -------- Beginning of SDigits2Digit ----------- " << endl;
645 fRunLoader->GetEvent(ievt);
646
36031625 647 fTreeS = fPMDLoader->TreeS();
b2b13ee2 648 AliPMDsdigit *pmdsdigit;
36031625 649 TBranch *branch = fTreeS->GetBranch("PMDSDigit");
b2b13ee2 650 branch->SetAddress(&fSDigits);
651
36031625 652 fTreeD = fPMDLoader->TreeD();
653 if (fTreeD == 0x0)
b2b13ee2 654 {
36031625 655 fPMDLoader->MakeTree("D");
656 fTreeD = fPMDLoader->TreeD();
b2b13ee2 657 }
658 Int_t bufsize = 16000;
36031625 659 fTreeD->Branch("PMDDigit", &fDigits, bufsize);
b2b13ee2 660
661 Int_t trno, det, smn;
662 Int_t cellno;
663 Float_t edep, adc;
664
36031625 665 Int_t nmodules = (Int_t) fTreeS->GetEntries();
b2b13ee2 666
667 for (Int_t imodule = 0; imodule < nmodules; imodule++)
668 {
36031625 669 fTreeS->GetEntry(imodule);
b2b13ee2 670 Int_t nentries = fSDigits->GetLast();
671 //cout << " nentries = " << nentries << endl;
672 for (Int_t ient = 0; ient < nentries+1; ient++)
673 {
674 pmdsdigit = (AliPMDsdigit*)fSDigits->UncheckedAt(ient);
675 trno = pmdsdigit->GetTrackNumber();
676 det = pmdsdigit->GetDetector();
677 smn = pmdsdigit->GetSMNumber();
678 cellno = pmdsdigit->GetCellNumber();
679 edep = pmdsdigit->GetCellEdep();
680
681 MeV2ADC(edep,adc);
b2b13ee2 682 AddDigit(trno,det,smn,cellno,adc);
683 }
36031625 684 fTreeD->Fill();
b2b13ee2 685 ResetDigit();
686 }
36031625 687 fPMDLoader->WriteDigits("OVERWRITE");
b2b13ee2 688 // cout << " -------- End of SDigits2Digit ----------- " << endl;
689}
690
691void AliPMDDigitizer::TrackAssignment2Cell()
692{
b2b13ee2 693 //
75bb524e 694 // This block assigns the cell id when there are
b2b13ee2 695 // multiple tracks in a cell according to the
696 // energy deposition
697 //
75bb524e 698 Bool_t jsort = false;
8599803d 699
b2b13ee2 700 Int_t i, j, k;
701
36031625 702 Float_t *fracEdp;
703 Float_t *trEdp;
75bb524e 704 Int_t *status1;
705 Int_t *status2;
706 Int_t *trnarray;
36031625 707 Int_t ****pmdTrack;
708 Float_t ****pmdEdep;
b2b13ee2 709
36031625 710 pmdTrack = new Int_t ***[fgkTotUM];
711 pmdEdep = new Float_t ***[fgkTotUM];
712 for (i=0; i<fgkTotUM; i++)
b2b13ee2 713 {
36031625 714 pmdTrack[i] = new Int_t **[fgkRow];
715 pmdEdep[i] = new Float_t **[fgkRow];
b2b13ee2 716 }
717
36031625 718 for (i = 0; i < fgkTotUM; i++)
b2b13ee2 719 {
36031625 720 for (j = 0; j < fgkRow; j++)
b2b13ee2 721 {
36031625 722 pmdTrack[i][j] = new Int_t *[fgkCol];
723 pmdEdep[i][j] = new Float_t *[fgkCol];
b2b13ee2 724 }
725 }
726
36031625 727 for (i = 0; i < fgkTotUM; i++)
b2b13ee2 728 {
36031625 729 for (j = 0; j < fgkRow; j++)
b2b13ee2 730 {
36031625 731 for (k = 0; k < fgkCol; k++)
b2b13ee2 732 {
36031625 733 Int_t nn = fPRECounter[i][j][k];
b2b13ee2 734 if(nn > 0)
735 {
36031625 736 pmdTrack[i][j][k] = new Int_t[nn];
737 pmdEdep[i][j][k] = new Float_t[nn];
b2b13ee2 738 }
739 else
740 {
741 nn = 1;
36031625 742 pmdTrack[i][j][k] = new Int_t[nn];
743 pmdEdep[i][j][k] = new Float_t[nn];
b2b13ee2 744 }
36031625 745 fPRECounter[i][j][k] = 0;
b2b13ee2 746 }
747 }
748 }
749
750
751 Int_t nentries = fCell->GetEntries();
752
753 Int_t mtrackno, ism, ixp, iyp;
754 Float_t edep;
755
756 for (i = 0; i < nentries; i++)
757 {
36031625 758 fPMDcell = (AliPMDcell*)fCell->UncheckedAt(i);
b2b13ee2 759
36031625 760 mtrackno = fPMDcell->GetTrackNumber();
761 ism = fPMDcell->GetSMNumber();
762 ixp = fPMDcell->GetX();
763 iyp = fPMDcell->GetY();
764 edep = fPMDcell->GetEdep();
765 Int_t nn = fPRECounter[ism][ixp][iyp];
b2b13ee2 766 // cout << " nn = " << nn << endl;
36031625 767 pmdTrack[ism][ixp][iyp][nn] = (Int_t) mtrackno;
768 pmdEdep[ism][ixp][iyp][nn] = edep;
769 fPRECounter[ism][ixp][iyp]++;
b2b13ee2 770 }
771
772 Int_t iz, il;
773 Int_t im, ix, iy;
774 Int_t nn;
775
36031625 776 for (im=0; im<fgkTotUM; im++)
b2b13ee2 777 {
36031625 778 for (ix=0; ix<fgkRow; ix++)
b2b13ee2 779 {
36031625 780 for (iy=0; iy<fgkCol; iy++)
b2b13ee2 781 {
36031625 782 nn = fPRECounter[im][ix][iy];
b2b13ee2 783 if (nn > 1)
784 {
785 // This block handles if a cell is fired
786 // many times by many tracks
75bb524e 787 status1 = new Int_t[nn];
788 status2 = new Int_t[nn];
789 trnarray = new Int_t[nn];
b2b13ee2 790 for (iz = 0; iz < nn; iz++)
791 {
36031625 792 status1[iz] = pmdTrack[im][ix][iy][iz];
b2b13ee2 793 }
75bb524e 794 TMath::Sort(nn,status1,status2,jsort);
36031625 795 Int_t trackOld = -99999;
796 Int_t track, trCount = 0;
b2b13ee2 797 for (iz = 0; iz < nn; iz++)
798 {
75bb524e 799 track = status1[status2[iz]];
36031625 800 if (trackOld != track)
b2b13ee2 801 {
36031625 802 trnarray[trCount] = track;
803 trCount++;
b2b13ee2 804 }
36031625 805 trackOld = track;
b2b13ee2 806 }
75bb524e 807 delete status1;
808 delete status2;
36031625 809 Float_t totEdp = 0.;
810 trEdp = new Float_t[trCount];
811 fracEdp = new Float_t[trCount];
812 for (il = 0; il < trCount; il++)
b2b13ee2 813 {
36031625 814 trEdp[il] = 0.;
75bb524e 815 track = trnarray[il];
b2b13ee2 816 for (iz = 0; iz < nn; iz++)
817 {
36031625 818 if (track == pmdTrack[im][ix][iy][iz])
b2b13ee2 819 {
36031625 820 trEdp[il] += pmdEdep[im][ix][iy][iz];
b2b13ee2 821 }
822 }
36031625 823 totEdp += trEdp[il];
b2b13ee2 824 }
36031625 825 Int_t ilOld = 0;
826 Float_t fracOld = 0.;
b2b13ee2 827
36031625 828 for (il = 0; il < trCount; il++)
b2b13ee2 829 {
36031625 830 fracEdp[il] = trEdp[il]/totEdp;
831 if (fracOld < fracEdp[il])
b2b13ee2 832 {
36031625 833 fracOld = fracEdp[il];
834 ilOld = il;
b2b13ee2 835 }
836 }
36031625 837 fPRETrackNo[im][ix][iy] = trnarray[ilOld];
838 delete fracEdp;
839 delete trEdp;
75bb524e 840 delete trnarray;
b2b13ee2 841 }
842 else if (nn == 1)
843 {
844 // This only handles if a cell is fired
845 // by only one track
846
36031625 847 fPRETrackNo[im][ix][iy] = pmdTrack[im][ix][iy][0];
b2b13ee2 848
849 }
850 else if (nn ==0)
851 {
852 // This is if no cell is fired
36031625 853 fPRETrackNo[im][ix][iy] = -999;
b2b13ee2 854 }
855 } // end of iy
856 } // end of ix
857 } // end of im
858
859 // Delete all the pointers
860
36031625 861 for (i = 0; i < fgkTotUM; i++)
b2b13ee2 862 {
36031625 863 for (j = 0; j < fgkRow; j++)
b2b13ee2 864 {
36031625 865 for (k = 0; k < fgkCol; k++)
b2b13ee2 866 {
36031625 867 delete [] pmdTrack[i][j][k];
868 delete [] pmdEdep[i][j][k];
b2b13ee2 869 }
870 }
871 }
872
36031625 873 for (i = 0; i < fgkTotUM; i++)
b2b13ee2 874 {
36031625 875 for (j = 0; j < fgkRow; j++)
b2b13ee2 876 {
36031625 877 delete [] pmdTrack[i][j];
878 delete [] pmdEdep[i][j];
b2b13ee2 879 }
880 }
881
36031625 882 for (i = 0; i < fgkTotUM; i++)
b2b13ee2 883 {
36031625 884 delete [] pmdTrack[i];
885 delete [] pmdEdep[i];
b2b13ee2 886 }
36031625 887 delete pmdTrack;
888 delete pmdEdep;
b2b13ee2 889 //
890 // End of the cell id assignment
891 //
892}
893
894
895void AliPMDDigitizer::MeV2ADC(Float_t mev, Float_t & adc)
896{
36031625 897 // This converts the simulated edep to ADC according to the
898 // Test Beam Data
b2b13ee2 899 // To be done
36031625 900 //
b2b13ee2 901 adc = mev*1.;
902}
903void AliPMDDigitizer::AddSDigit(Int_t trnumber, Int_t det, Int_t smnumber,
904 Int_t cellnumber, Float_t adc)
905{
36031625 906 // Add SDigit
907 //
b2b13ee2 908 TClonesArray &lsdigits = *fSDigits;
909 AliPMDsdigit *newcell;
910 newcell = new AliPMDsdigit(trnumber,det,smnumber,cellnumber,adc);
911 new(lsdigits[fNsdigit++]) AliPMDsdigit(newcell);
912 delete newcell;
913}
914
915void AliPMDDigitizer::AddDigit(Int_t trnumber, Int_t det, Int_t smnumber,
916 Int_t cellnumber, Float_t adc)
917{
36031625 918 // Add Digit
919 //
b2b13ee2 920 TClonesArray &ldigits = *fDigits;
921 AliPMDdigit *newcell;
922 newcell = new AliPMDdigit(trnumber,det,smnumber,cellnumber,adc);
923 new(ldigits[fNdigit++]) AliPMDdigit(newcell);
924 delete newcell;
925}
926
b2b13ee2 927void AliPMDDigitizer::SetZPosition(Float_t zpos)
928{
929 fZPos = zpos;
930}
931Float_t AliPMDDigitizer::GetZPosition() const
932{
933 return fZPos;
934}
935
936void AliPMDDigitizer::ResetCell()
937{
36031625 938 // clears the cell array and also the counter
939 // for each cell
940 //
b2b13ee2 941 fCell->Clear();
36031625 942 for (Int_t i = 0; i < fgkTotUM; i++)
b2b13ee2 943 {
36031625 944 for (Int_t j = 0; j < fgkRow; j++)
b2b13ee2 945 {
36031625 946 for (Int_t k = 0; k < fgkCol; k++)
b2b13ee2 947 {
36031625 948 fPRECounter[i][j][k] = 0;
b2b13ee2 949 }
950 }
951 }
952}
953void AliPMDDigitizer::ResetSDigit()
954{
36031625 955 // Clears SDigits
b2b13ee2 956 fNsdigit = 0;
957 if (fSDigits) fSDigits->Clear();
958}
959void AliPMDDigitizer::ResetDigit()
960{
36031625 961 // Clears Digits
b2b13ee2 962 fNdigit = 0;
963 if (fDigits) fDigits->Clear();
964}
965
966void AliPMDDigitizer::ResetCellADC()
967{
36031625 968 // Clears individual cells edep
969 for (Int_t i = 0; i < fgkTotUM; i++)
b2b13ee2 970 {
36031625 971 for (Int_t j = 0; j < fgkRow; j++)
b2b13ee2 972 {
36031625 973 for (Int_t k = 0; k < fgkCol; k++)
b2b13ee2 974 {
975 fCPV[i][j][k] = 0.;
36031625 976 fPRE[i][j][k] = 0.;
b2b13ee2 977 }
978 }
979 }
980}
981
982void AliPMDDigitizer::UnLoad(Option_t *option)
983{
36031625 984 // Unloads all the root files
985 //
b2b13ee2 986 const char *cS = strstr(option,"S");
987 const char *cD = strstr(option,"D");
988
989 fRunLoader->UnloadgAlice();
990 fRunLoader->UnloadHeader();
991 fRunLoader->UnloadKinematics();
992
993 if (cS)
994 {
36031625 995 fPMDLoader->UnloadHits();
b2b13ee2 996 }
997 if (cD)
998 {
36031625 999 fPMDLoader->UnloadHits();
1000 fPMDLoader->UnloadSDigits();
b2b13ee2 1001 }
1002}