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