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