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