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