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