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