1) First Year Geometry option added, name "EMCAL_FIRSTYEAR". Not useful yet until...
[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>
24c5571f 24#include <TBRIK.h>
25#include <TNode.h>
b2b13ee2 26#include <TTree.h>
24c5571f 27#include <TGeometry.h>
b2b13ee2 28#include <TObjArray.h>
29#include <TClonesArray.h>
30#include <TFile.h>
31#include <TNtuple.h>
32#include <TParticle.h>
e939a978 33#include <TRandom.h>
b2b13ee2 34
ecee2a1a 35#include "AliLog.h"
b2b13ee2 36#include "AliRun.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"
09a06455 46#include "AliCDBManager.h"
47#include "AliCDBStorage.h"
48#include "AliCDBEntry.h"
49#include "AliMC.h"
b2b13ee2 50
09a06455 51#include "AliPMD.h"
52#include "AliPMDhit.h"
b2b13ee2 53#include "AliPMDcell.h"
54#include "AliPMDsdigit.h"
55#include "AliPMDdigit.h"
09a06455 56#include "AliPMDCalibData.h"
35535af7 57#include "AliPMDPedestal.h"
b2b13ee2 58#include "AliPMDDigitizer.h"
09a06455 59
b2b13ee2 60
61ClassImp(AliPMDDigitizer)
a7d110b8 62
01960205 63AliPMDDigitizer::AliPMDDigitizer() :
64 fRunLoader(0),
65 fPMDHit(0),
66 fPMD(0),
67 fPMDLoader(0),
35535af7 68 fCalibGain(GetCalibGain()),
69 fCalibPed(GetCalibPed()),
ebd83c56 70 fSDigits(0),
71 fDigits(0),
4d59e381 72 fCPVCell(0),
ebd83c56 73 fCell(0),
01960205 74 fNsdigit(0),
75 fNdigit(0),
76 fDetNo(0),
2332574a 77 fZPos(361.5) // in units of cm, default position of PMD
b2b13ee2 78{
36031625 79 // Default Constructor
80 //
36031625 81 for (Int_t i = 0; i < fgkTotUM; i++)
b2b13ee2 82 {
36031625 83 for (Int_t j = 0; j < fgkRow; j++)
b2b13ee2 84 {
36031625 85 for (Int_t k = 0; k < fgkCol; k++)
b2b13ee2 86 {
920e13db 87 fCPV[i][j][k] = 0.;
88 fPRE[i][j][k] = 0.;
89 fCPVCounter[i][j][k] = 0;
90 fPRECounter[i][j][k] = 0;
91 fCPVTrackNo[i][j][k] = -1;
92 fPRETrackNo[i][j][k] = -1;
93 fCPVTrackPid[i][j][k] = -1;
94 fPRETrackPid[i][j][k] = -1;
b2b13ee2 95 }
96 }
97 }
2332574a 98
b2b13ee2 99}
7088dbe9 100//____________________________________________________________________________
2332574a 101AliPMDDigitizer::AliPMDDigitizer(const AliPMDDigitizer& digitizer):
102 AliDigitizer(digitizer),
103 fRunLoader(0),
104 fPMDHit(0),
105 fPMD(0),
106 fPMDLoader(0),
35535af7 107 fCalibGain(GetCalibGain()),
108 fCalibPed(GetCalibPed()),
2332574a 109 fSDigits(0),
110 fDigits(0),
4d59e381 111 fCPVCell(0),
2332574a 112 fCell(0),
113 fNsdigit(0),
114 fNdigit(0),
115 fDetNo(0),
116 fZPos(361.5) // in units of cm, default position of PMD
a48edddd 117{
118 // copy constructor
119 AliError("Copy constructor not allowed ");
120
121}
122//____________________________________________________________________________
a48edddd 123AliPMDDigitizer & AliPMDDigitizer::operator=(const AliPMDDigitizer& /*digitizer*/)
124{
125 // Assignment operator
126 AliError("Assignement operator not allowed ");
127
128 return *this;
129}
130//____________________________________________________________________________
2332574a 131AliPMDDigitizer::AliPMDDigitizer(AliRunDigitizer* manager):
132 AliDigitizer(manager),
ebd83c56 133 fRunLoader(0),
134 fPMDHit(0),
135 fPMD(0),
136 fPMDLoader(0),
35535af7 137 fCalibGain(GetCalibGain()),
138 fCalibPed(GetCalibPed()),
ebd83c56 139 fSDigits(new TClonesArray("AliPMDsdigit", 1000)),
140 fDigits(new TClonesArray("AliPMDdigit", 1000)),
4d59e381 141 fCPVCell(0),
ebd83c56 142 fCell(0),
ebd83c56 143 fNsdigit(0),
144 fNdigit(0),
145 fDetNo(0),
146 fZPos(361.5)// in units of cm, This is the default position of PMD
7088dbe9 147{
148 // ctor which should be used
149
09a06455 150
ebd83c56 151 for (Int_t i = 0; i < fgkTotUM; i++)
152 {
153 for (Int_t j = 0; j < fgkRow; j++)
154 {
155 for (Int_t k = 0; k < fgkCol; k++)
156 {
920e13db 157 fCPV[i][j][k] = 0.;
158 fPRE[i][j][k] = 0.;
159 fCPVCounter[i][j][k] = 0;
160 fPRECounter[i][j][k] = 0;
161 fCPVTrackNo[i][j][k] = -1;
162 fPRETrackNo[i][j][k] = -1;
163 fCPVTrackPid[i][j][k] = -1;
164 fPRETrackPid[i][j][k] = -1;
ebd83c56 165 }
166 }
167 }
7088dbe9 168}
2332574a 169
7088dbe9 170//____________________________________________________________________________
b2b13ee2 171AliPMDDigitizer::~AliPMDDigitizer()
172{
36031625 173 // Default Destructor
174 //
01960205 175 if (fSDigits) {
176 fSDigits->Delete();
177 delete fSDigits;
178 fSDigits=0;
179 }
180 if (fDigits) {
181 fDigits->Delete();
182 delete fDigits;
183 fDigits=0;
184 }
4d59e381 185 fCPVCell.Delete();
ebd83c56 186 fCell.Delete();
b2b13ee2 187}
188//
189// Member functions
190//
7088dbe9 191//____________________________________________________________________________
85a5290f 192void AliPMDDigitizer::OpengAliceFile(const char *file, Option_t *option)
b2b13ee2 193{
36031625 194 // Loads galice.root file and corresponding header, kinematics
195 // hits and sdigits or digits depending on the option
196 //
f540341d 197
e191bb57 198 TString evfoldname = AliConfig::GetDefaultEventFolderName();
f540341d 199 fRunLoader = AliRunLoader::GetRunLoader(evfoldname);
200 if (!fRunLoader)
5e2cc5bc 201 fRunLoader = AliRunLoader::Open(file,AliConfig::GetDefaultEventFolderName(), "UPDATE");
b2b13ee2 202
203 if (!fRunLoader)
204 {
ecee2a1a 205 AliError(Form("Can not open session for file %s.",file));
b2b13ee2 206 }
b2b13ee2 207
be9188a5 208 const char *cHS = strstr(option,"HS");
209 const char *cHD = strstr(option,"HD");
210 const char *cSD = strstr(option,"SD");
b2b13ee2 211
be9188a5 212 if(cHS || cHD)
ecee2a1a 213 {
be9188a5 214 if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
215 if (!fRunLoader->TreeE()) fRunLoader->LoadHeader();
216 if (!fRunLoader->TreeK()) fRunLoader->LoadKinematics();
217
218 gAlice = fRunLoader->GetAliRun();
219
220 if (gAlice)
221 {
222 AliDebug(1,"Alirun object found");
223 }
224 else
225 {
226 AliError("Could not found Alirun object");
227 }
f540341d 228
be9188a5 229 fPMD = (AliPMD*)gAlice->GetDetector("PMD");
230 }
231
36031625 232 fPMDLoader = fRunLoader->GetLoader("PMDLoader");
233 if (fPMDLoader == 0x0)
b2b13ee2 234 {
ecee2a1a 235 AliError("Can not find PMDLoader");
b2b13ee2 236 }
237
b2b13ee2 238
239 if (cHS)
240 {
36031625 241 fPMDLoader->LoadHits("READ");
242 fPMDLoader->LoadSDigits("recreate");
b2b13ee2 243 }
244 else if (cHD)
245 {
36031625 246 fPMDLoader->LoadHits("READ");
247 fPMDLoader->LoadDigits("recreate");
b2b13ee2 248 }
249 else if (cSD)
250 {
36031625 251 fPMDLoader->LoadSDigits("READ");
252 fPMDLoader->LoadDigits("recreate");
b2b13ee2 253 }
b2b13ee2 254}
7088dbe9 255//____________________________________________________________________________
b2b13ee2 256void AliPMDDigitizer::Hits2SDigits(Int_t ievt)
257{
36031625 258 // This reads the PMD Hits tree and assigns the right track number
259 // to a cell and stores in the summable digits tree
260 //
b2b13ee2 261
36031625 262 const Int_t kPi0 = 111;
263 const Int_t kGamma = 22;
b2b13ee2 264 Int_t npmd;
265 Int_t trackno;
b2b13ee2 266 Int_t smnumber;
267 Int_t trackpid;
268 Int_t mtrackno;
269 Int_t mtrackpid;
270
271 Float_t xPos, yPos, zPos;
8599803d 272 Int_t xpad = -1, ypad = -1;
b2b13ee2 273 Float_t edep;
274 Float_t vx = -999.0, vy = -999.0, vz = -999.0;
275
a7d110b8 276
ebd83c56 277 if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
b2b13ee2 278 ResetSDigit();
279
ecee2a1a 280 AliDebug(1,Form("Event Number = %d",ievt));
b2b13ee2 281 Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
ecee2a1a 282 AliDebug(1,Form("Number of Particles = %d",nparticles));
8c24e5f2 283
920e13db 284 //
285
b2b13ee2 286 fRunLoader->GetEvent(ievt);
b2b13ee2 287 // ------------------------------------------------------- //
288 // Pointer to specific detector hits.
289 // Get pointers to Alice detectors and Hits containers
290
ebd83c56 291 TTree* treeH = fPMDLoader->TreeH();
b2b13ee2 292
ebd83c56 293 Int_t ntracks = (Int_t) treeH->GetEntries();
ecee2a1a 294 AliDebug(1,Form("Number of Tracks in the TreeH = %d", ntracks));
ebd83c56 295 TTree* treeS = fPMDLoader->TreeS();
296 if (treeS == 0x0)
b2b13ee2 297 {
36031625 298 fPMDLoader->MakeTree("S");
ebd83c56 299 treeS = fPMDLoader->TreeS();
b2b13ee2 300 }
301 Int_t bufsize = 16000;
5e2cc5bc 302 treeS->Branch("PMDSDigit", &fSDigits, bufsize);
b2b13ee2 303
ebd83c56 304 TClonesArray* hits = 0;
305 if (fPMD) hits = fPMD->Hits();
b2b13ee2 306
307 // Start loop on tracks in the hits containers
a7d110b8 308
5e2cc5bc 309 for (Int_t track=0; track<ntracks;track++)
b2b13ee2 310 {
cad3294f 311 gAlice->GetMCApp()->ResetHits();
ebd83c56 312 treeH->GetEvent(track);
5e2cc5bc 313 if (fPMD)
b2b13ee2 314 {
ebd83c56 315 npmd = hits->GetEntriesFast();
4d59e381 316 for (Int_t ipmd = 0; ipmd < npmd; ipmd++)
b2b13ee2 317 {
ebd83c56 318 fPMDHit = (AliPMDhit*) hits->UncheckedAt(ipmd);
36031625 319 trackno = fPMDHit->GetTrack();
b2b13ee2 320 // get kinematics of the particles
36031625 321
ebd83c56 322 TParticle* mparticle = gAlice->GetMCApp()->Particle(trackno);
323 trackpid = mparticle->GetPdgCode();
4d59e381 324 Int_t ks = mparticle->GetStatusCode();
b2b13ee2 325 Int_t imo;
36031625 326 Int_t tracknoOld=0, trackpidOld=0, statusOld = 0;
4d59e381 327
8599803d 328 if (mparticle->GetFirstMother() == -1)
329 {
36031625 330 tracknoOld = trackno;
331 trackpidOld = trackpid;
332 statusOld = -1;
8599803d 333 }
8599803d 334 Int_t igstatus = 0;
4d59e381 335 //------------------modified by Mriganka ----------------------
336 if(ks==1||(imo = mparticle->GetFirstMother())<0 ){
337 vx = mparticle->Vx();
338 vy = mparticle->Vy();
339 vz = mparticle->Vz();
340
341 if(trackpid==kGamma||trackpid==11||trackpid==-11||
342 trackpid==kPi0)igstatus=1;
343 }
24c5571f 344
345
4d59e381 346 while(((imo = mparticle->GetFirstMother()) >= 0)&&
347 (ks = mparticle->GetStatusCode() <1) )
5f7ce82b 348 {
b208c6a3 349 mparticle = gAlice->GetMCApp()->Particle(imo);
4d59e381 350 trackpid = mparticle->GetPdgCode();
351 ks = mparticle->GetStatusCode();
b2b13ee2 352 vx = mparticle->Vx();
353 vy = mparticle->Vy();
354 vz = mparticle->Vz();
24c5571f 355
356 trackno=imo;
357
5f7ce82b 358 }
24c5571f 359
5f7ce82b 360 if(trackpid==kGamma||trackpid==11||trackpid==-11||
361 trackpid==kPi0)igstatus=1;
362 mtrackpid=trackpid;
363 mtrackno=trackno;
364 trackpid=trackpidOld;
365 trackno=tracknoOld;
d934d257 366
24c5571f 367 //-----------------end of modification----------------
36031625 368 xPos = fPMDHit->X();
369 yPos = fPMDHit->Y();
370 zPos = fPMDHit->Z();
7088dbe9 371
c35544fe 372 edep = fPMDHit->GetEnergy();
373 Int_t vol1 = fPMDHit->GetVolume(1); // Column
374 Int_t vol2 = fPMDHit->GetVolume(2); // Row
24c5571f 375 Int_t vol7 = fPMDHit->GetVolume(4); // Serial Module No
f117e3aa 376
c35544fe 377
8599803d 378 // -----------------------------------------//
f117e3aa 379 // In new geometry after adding electronics //
8599803d 380 // For Super Module 1 & 2 //
8599803d 381 // nrow = 48, ncol = 96 //
f117e3aa 382 // For Super Module 3 & 4 //
383 // nrow = 96, ncol = 48 //
8599803d 384 // -----------------------------------------//
f117e3aa 385
24c5571f 386 if (vol7 < 24)
387 {
388 smnumber = vol7;
389 }
390 else
391 {
392 smnumber = vol7 - 24;
393 }
394 Int_t vol8 = smnumber/6 + 1; // fake supermodule
8599803d 395
f117e3aa 396 if (vol8 == 1 || vol8 == 2)
b2b13ee2 397 {
36031625 398 xpad = vol2;
399 ypad = vol1;
b2b13ee2 400 }
f117e3aa 401 else if (vol8 == 3 || vol8 == 4)
402 {
403 xpad = vol1;
404 ypad = vol2;
405 }
b2b13ee2 406
ecee2a1a 407 AliDebug(2,Form("Zposition = %f Edeposition = %f",zPos,edep));
5f7ce82b 408 //Float_t zposition = TMath::Abs(zPos);
f117e3aa 409 if (zPos < fZPos)
b2b13ee2 410 {
411 // CPV
412 fDetNo = 1;
413 }
f117e3aa 414 else if (zPos > fZPos)
b2b13ee2 415 {
416 // PMD
417 fDetNo = 0;
418 }
24c5571f 419 //Int_t smn = smnumber - 1;
420 Int_t smn = smnumber;
8599803d 421 Int_t ixx = xpad - 1;
422 Int_t iyy = ypad - 1;
b2b13ee2 423 if (fDetNo == 0)
424 {
36031625 425 fPRE[smn][ixx][iyy] += edep;
426 fPRECounter[smn][ixx][iyy]++;
b2b13ee2 427
ebd83c56 428 AliPMDcell* cell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
ebd83c56 429 fCell.Add(cell);
b2b13ee2 430 }
431 else if(fDetNo == 1)
432 {
8599803d 433 fCPV[smn][ixx][iyy] += edep;
4d59e381 434 fCPVCounter[smn][ixx][iyy]++;
435 AliPMDcell* cpvcell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
436 fCPVCell.Add(cpvcell);
b2b13ee2 437 }
438 }
439 }
440 } // Track Loop ended
4d59e381 441 TrackAssignment2CPVCell();
b2b13ee2 442 TrackAssignment2Cell();
b2b13ee2 443 ResetCell();
444
445 Float_t deltaE = 0.;
446 Int_t detno = 0;
b2b13ee2 447 Int_t trno = -1;
920e13db 448 Int_t trpid = -99;
b2b13ee2 449
450 for (Int_t idet = 0; idet < 2; idet++)
451 {
36031625 452 for (Int_t ism = 0; ism < fgkTotUM; ism++)
b2b13ee2 453 {
36031625 454 for (Int_t jrow = 0; jrow < fgkRow; jrow++)
b2b13ee2 455 {
36031625 456 for (Int_t kcol = 0; kcol < fgkCol; kcol++)
b2b13ee2 457 {
b2b13ee2 458 if (idet == 0)
459 {
36031625 460 deltaE = fPRE[ism][jrow][kcol];
461 trno = fPRETrackNo[ism][jrow][kcol];
b2b13ee2 462 detno = 0;
463 }
464 else if (idet == 1)
465 {
466 deltaE = fCPV[ism][jrow][kcol];
467 trno = fCPVTrackNo[ism][jrow][kcol];
468 detno = 1;
469 }
470 if (deltaE > 0.)
471 {
8c24e5f2 472 // Natasha
473 TParticle *mparticle = gAlice->GetMCApp()->Particle(trno);
474 trpid = mparticle->GetPdgCode();
920e13db 475 AddSDigit(trno,trpid,detno,ism,jrow,kcol,deltaE);
b2b13ee2 476 }
477 }
478 }
ebd83c56 479 treeS->Fill();
b2b13ee2 480 ResetSDigit();
481 }
482 }
36031625 483 fPMDLoader->WriteSDigits("OVERWRITE");
b2b13ee2 484 ResetCellADC();
b2b13ee2 485}
7088dbe9 486//____________________________________________________________________________
b2b13ee2 487
488void AliPMDDigitizer::Hits2Digits(Int_t ievt)
489{
36031625 490 // This reads the PMD Hits tree and assigns the right track number
491 // to a cell and stores in the digits tree
492 //
493 const Int_t kPi0 = 111;
494 const Int_t kGamma = 22;
b2b13ee2 495 Int_t npmd;
496 Int_t trackno;
b2b13ee2 497 Int_t smnumber;
498 Int_t trackpid;
499 Int_t mtrackno;
500 Int_t mtrackpid;
501
502 Float_t xPos, yPos, zPos;
8599803d 503 Int_t xpad = -1, ypad = -1;
b2b13ee2 504 Float_t edep;
505 Float_t vx = -999.0, vy = -999.0, vz = -999.0;
506
ebd83c56 507 if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
b2b13ee2 508 ResetDigit();
509
5e2cc5bc 510 AliDebug(1,Form("Event Number = %d",ievt));
b2b13ee2 511 Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
ecee2a1a 512 AliDebug(1,Form("Number of Particles = %d", nparticles));
920e13db 513
b2b13ee2 514 fRunLoader->GetEvent(ievt);
b2b13ee2 515 // ------------------------------------------------------- //
516 // Pointer to specific detector hits.
517 // Get pointers to Alice detectors and Hits containers
518
b208c6a3 519 fPMD = (AliPMD*)gAlice->GetDetector("PMD");
36031625 520 fPMDLoader = fRunLoader->GetLoader("PMDLoader");
b2b13ee2 521
36031625 522 if (fPMDLoader == 0x0)
b2b13ee2 523 {
ecee2a1a 524 AliError("Can not find PMD or PMDLoader");
b2b13ee2 525 }
ebd83c56 526 TTree* treeH = fPMDLoader->TreeH();
527 Int_t ntracks = (Int_t) treeH->GetEntries();
ecee2a1a 528 AliDebug(1,Form("Number of Tracks in the TreeH = %d", ntracks));
36031625 529 fPMDLoader->LoadDigits("recreate");
ebd83c56 530 TTree* treeD = fPMDLoader->TreeD();
531 if (treeD == 0x0)
b2b13ee2 532 {
36031625 533 fPMDLoader->MakeTree("D");
ebd83c56 534 treeD = fPMDLoader->TreeD();
b2b13ee2 535 }
536 Int_t bufsize = 16000;
5e2cc5bc 537 treeD->Branch("PMDDigit", &fDigits, bufsize);
b2b13ee2 538
ebd83c56 539 TClonesArray* hits = 0;
540 if (fPMD) hits = fPMD->Hits();
b2b13ee2 541
542 // Start loop on tracks in the hits containers
543
5e2cc5bc 544 for (Int_t track=0; track<ntracks;track++)
b2b13ee2 545 {
cad3294f 546 gAlice->GetMCApp()->ResetHits();
ebd83c56 547 treeH->GetEvent(track);
b2b13ee2 548
5e2cc5bc 549 if (fPMD)
b2b13ee2 550 {
ebd83c56 551 npmd = hits->GetEntriesFast();
4d59e381 552 for (Int_t ipmd = 0; ipmd < npmd; ipmd++)
b2b13ee2 553 {
ebd83c56 554 fPMDHit = (AliPMDhit*) hits->UncheckedAt(ipmd);
36031625 555 trackno = fPMDHit->GetTrack();
b2b13ee2 556
557 // get kinematics of the particles
558
ebd83c56 559 TParticle* mparticle = gAlice->GetMCApp()->Particle(trackno);
560 trackpid = mparticle->GetPdgCode();
4d59e381 561 Int_t ks = mparticle->GetStatusCode();
b2b13ee2 562 Int_t imo;
36031625 563 Int_t tracknoOld=0, trackpidOld=0, statusOld = 0;
8599803d 564 if (mparticle->GetFirstMother() == -1)
565 {
36031625 566 tracknoOld = trackno;
567 trackpidOld = trackpid;
568 statusOld = -1;
8599803d 569 }
570
571 Int_t igstatus = 0;
4d59e381 572 //-----------------------modified by Mriganka ------------------
573 if(ks==1||(imo = mparticle->GetFirstMother())<0 ){
574 vx = mparticle->Vx();
575 vy = mparticle->Vy();
576 vz = mparticle->Vz();
577
578 if(trackpid==kGamma||trackpid==11||trackpid==-11||trackpid==kPi0)
579 igstatus=1;
580 }
581
582
583 while(((imo = mparticle->GetFirstMother()) >= 0)&&
584 (ks = mparticle->GetStatusCode() <1) )
b2b13ee2 585 {
b208c6a3 586 mparticle = gAlice->GetMCApp()->Particle(imo);
4d59e381 587 trackpid = mparticle->GetPdgCode();
588 ks = mparticle->GetStatusCode();
b2b13ee2 589 vx = mparticle->Vx();
590 vy = mparticle->Vy();
591 vz = mparticle->Vz();
24c5571f 592
593 trackno=imo;
594
d934d257 595 }
4d59e381 596
597 if(trackpid==kGamma||trackpid==11||trackpid==-11||trackpid==kPi0)
598 igstatus=1;
599 mtrackpid=trackpid;
600 mtrackno=trackno;
601 trackpid=trackpidOld;
602 trackno=tracknoOld;
b2b13ee2 603
4d59e381 604 //-----------------end of modification----------------
36031625 605 xPos = fPMDHit->X();
606 yPos = fPMDHit->Y();
607 zPos = fPMDHit->Z();
f117e3aa 608 edep = fPMDHit->GetEnergy();
c35544fe 609 Int_t vol1 = fPMDHit->GetVolume(1); // Column
610 Int_t vol2 = fPMDHit->GetVolume(2); // Row
24c5571f 611 Int_t vol7 = fPMDHit->GetVolume(4); // Serial Module No
8599803d 612
613 // -----------------------------------------//
f117e3aa 614 // In new geometry after adding electronics //
8599803d 615 // For Super Module 1 & 2 //
8599803d 616 // nrow = 48, ncol = 96 //
f117e3aa 617 // For Super Module 3 & 4 //
618 // nrow = 96, ncol = 48 //
8599803d 619 // -----------------------------------------//
b2b13ee2 620
24c5571f 621 if (vol7 < 24)
622 {
623 smnumber = vol7;
624 }
625 else
626 {
627 smnumber = vol7 - 24;
628 }
629 Int_t vol8 = smnumber/6 + 1; // fake supermodule
8599803d 630
f117e3aa 631 if (vol8 == 1 || vol8 == 2)
8599803d 632 {
36031625 633 xpad = vol2;
634 ypad = vol1;
b2b13ee2 635 }
f117e3aa 636 else if (vol8 == 3 || vol8 == 4)
637 {
638 xpad = vol1;
639 ypad = vol2;
640 }
b2b13ee2 641
ecee2a1a 642 AliDebug(2,Form("ZPosition = %f Edeposition = %d",zPos,edep));
f117e3aa 643 //Float_t zposition = TMath::Abs(zPos);
8599803d 644
f117e3aa 645 if (zPos < fZPos)
b2b13ee2 646 {
647 // CPV
648 fDetNo = 1;
649 }
f117e3aa 650 else if (zPos > fZPos)
b2b13ee2 651 {
652 // PMD
653 fDetNo = 0;
654 }
8599803d 655
24c5571f 656 //Int_t smn = smnumber - 1;
657 Int_t smn = smnumber;
8599803d 658 Int_t ixx = xpad - 1;
659 Int_t iyy = ypad - 1;
660 if (fDetNo == 0)
b2b13ee2 661 {
36031625 662 fPRE[smn][ixx][iyy] += edep;
663 fPRECounter[smn][ixx][iyy]++;
8599803d 664
ebd83c56 665 AliPMDcell* cell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
8599803d 666
ebd83c56 667 fCell.Add(cell);
b2b13ee2 668 }
8599803d 669 else if(fDetNo == 1)
b2b13ee2 670 {
8599803d 671 fCPV[smn][ixx][iyy] += edep;
4d59e381 672 fCPVCounter[smn][ixx][iyy]++;
673 AliPMDcell* cpvcell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
674 fCPVCell.Add(cpvcell);
b2b13ee2 675 }
676 }
677 }
678 } // Track Loop ended
4d59e381 679 TrackAssignment2CPVCell();
b2b13ee2 680 TrackAssignment2Cell();
681 ResetCell();
682
50555ba1 683 Float_t gain1;
7088dbe9 684 Float_t adc;
b2b13ee2 685 Float_t deltaE = 0.;
686 Int_t detno = 0;
b2b13ee2 687 Int_t trno = 1;
920e13db 688 Int_t trpid = -99;
689
b2b13ee2 690 for (Int_t idet = 0; idet < 2; idet++)
5e2cc5bc 691 {
36031625 692 for (Int_t ism = 0; ism < fgkTotUM; ism++)
5e2cc5bc 693 {
36031625 694 for (Int_t jrow = 0; jrow < fgkRow; jrow++)
5e2cc5bc 695 {
36031625 696 for (Int_t kcol = 0; kcol < fgkCol; kcol++)
5e2cc5bc 697 {
b2b13ee2 698 if (idet == 0)
5e2cc5bc 699 {
700 deltaE = fPRE[ism][jrow][kcol];
36031625 701 trno = fPRETrackNo[ism][jrow][kcol];
920e13db 702 detno = 0;
5e2cc5bc 703 }
b2b13ee2 704 else if (idet == 1)
5e2cc5bc 705 {
706 deltaE = fCPV[ism][jrow][kcol];
8599803d 707 trno = fCPVTrackNo[ism][jrow][kcol];
920e13db 708 detno = 1;
5e2cc5bc 709 }
b2b13ee2 710 if (deltaE > 0.)
5e2cc5bc 711 {
7088dbe9 712 MeV2ADC(deltaE,adc);
5e2cc5bc 713
714 // To decalibrate the adc values
715 //
716 gain1 = Gain(idet,ism,jrow,kcol);
717 if (gain1 != 0.)
718 {
719 Int_t adcDecalib = (Int_t)(adc/gain1);
720 adc = (Float_t) adcDecalib;
721 }
722 else if(gain1 == 0.)
723 {
724 adc = 0.;
725 }
35535af7 726
727 // Pedestal Decalibration
728 Int_t pedmeanrms =
729 fCalibPed->GetPedMeanRms(idet,ism,jrow,kcol);
37edc588 730 Int_t pedrms1 = (Int_t) pedmeanrms%100;
35535af7 731 Float_t pedrms = (Float_t)pedrms1/10.;
732 Float_t pedmean =
733 (Float_t) (pedmeanrms - pedrms1)/1000.0;
35535af7 734 if (adc > 0.)
735 {
736 adc += (pedmean + 3.0*pedrms);
8c24e5f2 737 TParticle *mparticle
738 = gAlice->GetMCApp()->Particle(trno);
739 trpid = mparticle->GetPdgCode();
740
920e13db 741 AddDigit(trno,trpid,detno,ism,jrow,kcol,adc);
35535af7 742 }
5e2cc5bc 743 }
744 } // column loop
745 } // row loop
03965b13 746 treeD->Fill();
747 ResetDigit();
5e2cc5bc 748 } // supermodule loop
749 } // detector loop
03965b13 750
36031625 751 fPMDLoader->WriteDigits("OVERWRITE");
b2b13ee2 752 ResetCellADC();
920e13db 753
b2b13ee2 754}
7088dbe9 755//____________________________________________________________________________
b2b13ee2 756
757
758void AliPMDDigitizer::SDigits2Digits(Int_t ievt)
759{
36031625 760 // This reads the PMD sdigits tree and converts energy deposition
761 // in a cell to ADC and stores in the digits tree
762 //
ecee2a1a 763
b2b13ee2 764 fRunLoader->GetEvent(ievt);
765
ebd83c56 766 TTree* treeS = fPMDLoader->TreeS();
b2b13ee2 767 AliPMDsdigit *pmdsdigit;
ebd83c56 768 TBranch *branch = treeS->GetBranch("PMDSDigit");
ecee2a1a 769 if(!branch)
770 {
771 AliError("PMD Sdigit branch does not exist");
772 return;
773 }
ebd83c56 774 if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
b2b13ee2 775 branch->SetAddress(&fSDigits);
776
ebd83c56 777 TTree* treeD = fPMDLoader->TreeD();
778 if (treeD == 0x0)
b2b13ee2 779 {
36031625 780 fPMDLoader->MakeTree("D");
ebd83c56 781 treeD = fPMDLoader->TreeD();
b2b13ee2 782 }
783 Int_t bufsize = 16000;
ebd83c56 784 if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
5e2cc5bc 785 treeD->Branch("PMDDigit", &fDigits, bufsize);
b2b13ee2 786
920e13db 787 Int_t trno, trpid, det, smn;
7088dbe9 788 Int_t irow, icol;
b2b13ee2 789 Float_t edep, adc;
790
ebd83c56 791 Int_t nmodules = (Int_t) treeS->GetEntries();
ecee2a1a 792 AliDebug(1,Form("Number of modules = %d",nmodules));
b2b13ee2 793
794 for (Int_t imodule = 0; imodule < nmodules; imodule++)
795 {
5e2cc5bc 796 treeS->GetEntry(imodule);
b2b13ee2 797 Int_t nentries = fSDigits->GetLast();
ecee2a1a 798 AliDebug(2,Form("Number of entries per module = %d",nentries+1));
b2b13ee2 799 for (Int_t ient = 0; ient < nentries+1; ient++)
800 {
801 pmdsdigit = (AliPMDsdigit*)fSDigits->UncheckedAt(ient);
802 trno = pmdsdigit->GetTrackNumber();
920e13db 803 trpid = pmdsdigit->GetTrackPid();
b2b13ee2 804 det = pmdsdigit->GetDetector();
805 smn = pmdsdigit->GetSMNumber();
7088dbe9 806 irow = pmdsdigit->GetRow();
807 icol = pmdsdigit->GetColumn();
b2b13ee2 808 edep = pmdsdigit->GetCellEdep();
809
810 MeV2ADC(edep,adc);
5e2cc5bc 811
812 // To decalibrte the adc values
813 //
814 Float_t gain1 = Gain(det,smn,irow,icol);
815 if (gain1 != 0.)
816 {
04a35d8d 817 Int_t adcDecalib = (Int_t)(adc/gain1);
818 adc = (Float_t) adcDecalib;
5e2cc5bc 819 }
820 else if(gain1 == 0.)
821 {
822 adc = 0.;
823 }
35535af7 824 // Pedestal Decalibration
825 Int_t pedmeanrms = fCalibPed->GetPedMeanRms(det,smn,irow,icol);
37edc588 826 Int_t pedrms1 = (Int_t) pedmeanrms%100;
35535af7 827 Float_t pedrms = (Float_t)pedrms1/10.;
828 Float_t pedmean = (Float_t) (pedmeanrms - pedrms1)/1000.0;
35535af7 829 if(adc > 0.)
830 {
831 adc += (pedmean + 3.0*pedrms);
920e13db 832 AddDigit(trno,trpid,det,smn,irow,icol,adc);
35535af7 833 }
5e2cc5bc 834
b2b13ee2 835 }
ebd83c56 836 treeD->Fill();
b2b13ee2 837 ResetDigit();
838 }
36031625 839 fPMDLoader->WriteDigits("OVERWRITE");
ecee2a1a 840
b2b13ee2 841}
7088dbe9 842//____________________________________________________________________________
843void AliPMDDigitizer::Exec(Option_t *option)
844{
845 // Does the event merging and digitization
7088dbe9 846 const char *cdeb = strstr(option,"deb");
847 if(cdeb)
848 {
ecee2a1a 849 AliDebug(100," *** PMD Exec is called ***");
7088dbe9 850 }
ecee2a1a 851
7088dbe9 852 Int_t ninputs = fManager->GetNinputs();
ecee2a1a 853 AliDebug(1,Form("Number of files to be processed = %d",ninputs));
7088dbe9 854 ResetCellADC();
855
856 for (Int_t i = 0; i < ninputs; i++)
857 {
858 Int_t troffset = fManager->GetMask(i);
859 MergeSDigits(i, troffset);
860 }
b2b13ee2 861
7088dbe9 862 fRunLoader = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
863 fPMD = (AliPMD*)gAlice->GetDetector("PMD");
864 fPMDLoader = fRunLoader->GetLoader("PMDLoader");
865 if (fPMDLoader == 0x0)
866 {
ecee2a1a 867 AliError("Can not find PMD or PMDLoader");
7088dbe9 868 }
06f2c79d 869 fPMDLoader->LoadDigits("update");
ebd83c56 870 TTree* treeD = fPMDLoader->TreeD();
871 if (treeD == 0x0)
7088dbe9 872 {
873 fPMDLoader->MakeTree("D");
ebd83c56 874 treeD = fPMDLoader->TreeD();
7088dbe9 875 }
876 Int_t bufsize = 16000;
ebd83c56 877 if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
5e2cc5bc 878 treeD->Branch("PMDDigit", &fDigits, bufsize);
7088dbe9 879
880 Float_t adc;
881 Float_t deltaE = 0.;
920e13db 882 Int_t detno = 0;
883 Int_t trno = 1;
884 Int_t trpid = -99;
7088dbe9 885
886 for (Int_t idet = 0; idet < 2; idet++)
887 {
888 for (Int_t ism = 0; ism < fgkTotUM; ism++)
889 {
890 for (Int_t jrow = 0; jrow < fgkRow; jrow++)
891 {
892 for (Int_t kcol = 0; kcol < fgkCol; kcol++)
893 {
894 if (idet == 0)
895 {
896 deltaE = fPRE[ism][jrow][kcol];
897 trno = fPRETrackNo[ism][jrow][kcol];
920e13db 898 trpid = fPRETrackPid[ism][jrow][kcol];
899 detno = 0;
7088dbe9 900 }
901 else if (idet == 1)
902 {
903 deltaE = fCPV[ism][jrow][kcol];
904 trno = fCPVTrackNo[ism][jrow][kcol];
920e13db 905 trpid = fCPVTrackPid[ism][jrow][kcol];
906 detno = 1;
7088dbe9 907 }
908 if (deltaE > 0.)
909 {
910 MeV2ADC(deltaE,adc);
35535af7 911
5f566ada 912 //
35535af7 913 // Gain decalibration
5f566ada 914 //
915 Float_t gain1 = Gain(idet,ism,jrow,kcol);
35535af7 916
5f566ada 917 if (gain1 != 0.)
918 {
919 Int_t adcDecalib = (Int_t)(adc/gain1);
920 adc = (Float_t) adcDecalib;
921 }
922 else if(gain1 == 0.)
923 {
924 adc = 0.;
925 }
35535af7 926 // Pedestal Decalibration
927 Int_t pedmeanrms =
928 fCalibPed->GetPedMeanRms(idet,ism,jrow,kcol);
37edc588 929 Int_t pedrms1 = (Int_t) pedmeanrms%100;
35535af7 930 Float_t pedrms = (Float_t)pedrms1/10.;
931 Float_t pedmean =
932 (Float_t) (pedmeanrms - pedrms1)/1000.0;
35535af7 933 if (adc > 0.)
934 {
935 adc += (pedmean + 3.0*pedrms);
920e13db 936 AddDigit(trno,trpid,detno,ism,jrow,kcol,adc);
35535af7 937 }
5f566ada 938
7088dbe9 939 }
940 } // column loop
941 } // row loop
ebd83c56 942 treeD->Fill();
7088dbe9 943 ResetDigit();
944 } // supermodule loop
945 } // detector loop
5e2cc5bc 946 fPMDLoader->WriteDigits("OVERWRITE");
06f2c79d 947 fPMDLoader->UnloadDigits();
948 ResetCellADC();
7088dbe9 949}
950//____________________________________________________________________________
4d59e381 951void AliPMDDigitizer::TrackAssignment2CPVCell()
952{
953 // This block assigns the cell id when there are
954 // multiple tracks in a cell according to the
955 // energy deposition
956 // This method added by Ajay
957 Bool_t jsort = false;
958
959 Int_t i, j, k;
960
961 Int_t *status1;
962 Int_t *status2;
963 Int_t *trnarray;
964 Float_t *fracEdp;
965 Float_t *trEdp;
966
967 Int_t ****cpvTrack;
968 Float_t ****cpvEdep;
969
970 cpvTrack = new Int_t ***[fgkTotUM];
971 cpvEdep = new Float_t ***[fgkTotUM];
972 for (i=0; i<fgkTotUM; i++)
973 {
974 cpvTrack[i] = new Int_t **[fgkRow];
975 cpvEdep[i] = new Float_t **[fgkRow];
976 }
977
978 for (i = 0; i < fgkTotUM; i++)
979 {
980 for (j = 0; j < fgkRow; j++)
981 {
982 cpvTrack[i][j] = new Int_t *[fgkCol];
983 cpvEdep[i][j] = new Float_t *[fgkCol];
984 }
985 }
986 for (i = 0; i < fgkTotUM; i++)
987 {
988 for (j = 0; j < fgkRow; j++)
989 {
990 for (k = 0; k < fgkCol; k++)
991 {
992 Int_t nn = fCPVCounter[i][j][k];
993 if(nn > 0)
994 {
995 cpvTrack[i][j][k] = new Int_t[nn];
996 cpvEdep[i][j][k] = new Float_t[nn];
997 }
998 else
999 {
1000 nn = 1;
1001 cpvTrack[i][j][k] = new Int_t[nn];
1002 cpvEdep[i][j][k] = new Float_t[nn];
1003 }
1004 fCPVCounter[i][j][k] = 0;
1005 }
1006 }
1007 }
1008
1009
1010 Int_t nentries = fCPVCell.GetEntries();
1011
1012 Int_t mtrackno, ism, ixp, iyp;
1013 Float_t edep;
1014 for (i = 0; i < nentries; i++)
1015 {
1016 AliPMDcell* cpvcell = (AliPMDcell*)fCPVCell.UncheckedAt(i);
1017
1018 mtrackno = cpvcell->GetTrackNumber();
1019 ism = cpvcell->GetSMNumber();
1020 ixp = cpvcell->GetX();
1021 iyp = cpvcell->GetY();
1022 edep = cpvcell->GetEdep();
1023 Int_t nn = fCPVCounter[ism][ixp][iyp];
1024 cpvTrack[ism][ixp][iyp][nn] = (Int_t) mtrackno;
1025 cpvEdep[ism][ixp][iyp][nn] = edep;
1026 fCPVCounter[ism][ixp][iyp]++;
1027 }
1028
1029 Int_t iz, il;
1030 Int_t im, ix, iy;
1031 Int_t nn;
1032 for (im=0; im<fgkTotUM; im++)
1033 {
1034 for (ix=0; ix<fgkRow; ix++)
1035 {
1036 for (iy=0; iy<fgkCol; iy++)
1037 {
1038 nn = fCPVCounter[im][ix][iy];
1039 if (nn > 1)
1040 {
1041 // This block handles if a cell is fired
1042 // many times by many tracks
1043 status1 = new Int_t[nn];
1044 status2 = new Int_t[nn];
1045 trnarray = new Int_t[nn];
1046 for (iz = 0; iz < nn; iz++)
1047 {
1048 status1[iz] = cpvTrack[im][ix][iy][iz];
1049 }
1050 TMath::Sort(nn,status1,status2,jsort);
1051 Int_t trackOld = -99999;
1052 Int_t track, trCount = 0;
1053 for (iz = 0; iz < nn; iz++)
1054 {
1055 track = status1[status2[iz]];
1056 if (trackOld != track)
1057 {
1058 trnarray[trCount] = track;
1059 trCount++;
1060 }
1061 trackOld = track;
1062 }
1063 delete [] status1;
1064 delete [] status2;
1065 Float_t totEdp = 0.;
1066 trEdp = new Float_t[trCount];
1067 fracEdp = new Float_t[trCount];
1068 for (il = 0; il < trCount; il++)
1069 {
1070 trEdp[il] = 0.;
1071 track = trnarray[il];
1072 for (iz = 0; iz < nn; iz++)
1073 {
1074 if (track == cpvTrack[im][ix][iy][iz])
1075 {
1076 trEdp[il] += cpvEdep[im][ix][iy][iz];
1077 }
1078 }
1079 totEdp += trEdp[il];
1080 }
1081 Int_t ilOld = 0;
1082 Float_t fracOld = 0.;
1083
1084 for (il = 0; il < trCount; il++)
1085 {
1086 fracEdp[il] = trEdp[il]/totEdp;
1087 if (fracOld < fracEdp[il])
1088 {
1089 fracOld = fracEdp[il];
1090 ilOld = il;
1091 }
1092 }
1093 fCPVTrackNo[im][ix][iy] = trnarray[ilOld];
1094 delete [] fracEdp;
1095 delete [] trEdp;
1096 delete [] trnarray;
1097 }
1098 else if (nn == 1)
1099 {
1100 // This only handles if a cell is fired
1101 // by only one track
1102
1103 fCPVTrackNo[im][ix][iy] = cpvTrack[im][ix][iy][0];
1104
1105 }
1106 else if (nn ==0)
1107 {
1108 // This is if no cell is fired
1109 fCPVTrackNo[im][ix][iy] = -999;
1110 }
1111 } // end of iy
1112 } // end of ix
1113 } // end of im
1114
1115 // Delete all the pointers
1116
1117 for (i = 0; i < fgkTotUM; i++)
1118 {
1119 for (j = 0; j < fgkRow; j++)
1120 {
1121 for (k = 0; k < fgkCol; k++)
1122 {
1123 delete []cpvTrack[i][j][k];
1124 delete []cpvEdep[i][j][k];
1125 }
1126 }
1127 }
1128
1129 for (i = 0; i < fgkTotUM; i++)
1130 {
1131 for (j = 0; j < fgkRow; j++)
1132 {
1133 delete [] cpvTrack[i][j];
1134 delete [] cpvEdep[i][j];
1135 }
1136 }
1137
1138 for (i = 0; i < fgkTotUM; i++)
1139 {
1140 delete [] cpvTrack[i];
1141 delete [] cpvEdep[i];
1142 }
1143 delete [] cpvTrack;
1144 delete [] cpvEdep;
1145
1146 //
1147 // End of the cell id assignment
1148 //
1149}
1150//____________________________________________________________________________
7088dbe9 1151
1152void AliPMDDigitizer::MergeSDigits(Int_t filenumber, Int_t troffset)
1153{
1154 // merging sdigits
1155 fRunLoader = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(filenumber));
1156 fPMDLoader = fRunLoader->GetLoader("PMDLoader");
1157 fPMDLoader->LoadSDigits("read");
ebd83c56 1158 TTree* treeS = fPMDLoader->TreeS();
7088dbe9 1159 AliPMDsdigit *pmdsdigit;
ebd83c56 1160 TBranch *branch = treeS->GetBranch("PMDSDigit");
1161 if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
7088dbe9 1162 branch->SetAddress(&fSDigits);
1163
920e13db 1164 Int_t itrackno, itrackpid, idet, ism;
7088dbe9 1165 Int_t ixp, iyp;
1166 Float_t edep;
ebd83c56 1167 Int_t nmodules = (Int_t) treeS->GetEntries();
ecee2a1a 1168 AliDebug(1,Form("Number of Modules in the treeS = %d",nmodules));
1169 AliDebug(1,Form("Track Offset = %d",troffset));
7088dbe9 1170 for (Int_t imodule = 0; imodule < nmodules; imodule++)
1171 {
5e2cc5bc 1172 treeS->GetEntry(imodule);
7088dbe9 1173 Int_t nentries = fSDigits->GetLast();
ecee2a1a 1174 AliDebug(2,Form("Number of Entries per Module = %d",nentries));
7088dbe9 1175 for (Int_t ient = 0; ient < nentries+1; ient++)
1176 {
1177 pmdsdigit = (AliPMDsdigit*)fSDigits->UncheckedAt(ient);
1178 itrackno = pmdsdigit->GetTrackNumber();
920e13db 1179 itrackpid = pmdsdigit->GetTrackPid();
7088dbe9 1180 idet = pmdsdigit->GetDetector();
1181 ism = pmdsdigit->GetSMNumber();
1182 ixp = pmdsdigit->GetRow();
1183 iyp = pmdsdigit->GetColumn();
1184 edep = pmdsdigit->GetCellEdep();
7088dbe9 1185 if (idet == 0)
1186 {
1187 if (fPRE[ism][ixp][iyp] < edep)
1188 {
1189 fPRETrackNo[ism][ixp][iyp] = troffset + itrackno;
920e13db 1190 fPRETrackPid[ism][ixp][iyp] = itrackpid;
7088dbe9 1191 }
1192 fPRE[ism][ixp][iyp] += edep;
1193 }
1194 else if (idet == 1)
1195 {
1196 if (fCPV[ism][ixp][iyp] < edep)
1197 {
1198 fCPVTrackNo[ism][ixp][iyp] = troffset + itrackno;
920e13db 1199 fCPVTrackPid[ism][ixp][iyp] = itrackpid;
7088dbe9 1200 }
1201 fCPV[ism][ixp][iyp] += edep;
1202 }
1203 }
1204 }
1205
1206}
1207// ----------------------------------------------------------------------
b2b13ee2 1208void AliPMDDigitizer::TrackAssignment2Cell()
1209{
b2b13ee2 1210 //
75bb524e 1211 // This block assigns the cell id when there are
b2b13ee2 1212 // multiple tracks in a cell according to the
1213 // energy deposition
1214 //
75bb524e 1215 Bool_t jsort = false;
8599803d 1216
b2b13ee2 1217 Int_t i, j, k;
1218
4d59e381 1219 Int_t *status1;
1220 Int_t *status2;
1221 Int_t *trnarray;
36031625 1222 Float_t *fracEdp;
1223 Float_t *trEdp;
4d59e381 1224
36031625 1225 Int_t ****pmdTrack;
1226 Float_t ****pmdEdep;
b2b13ee2 1227
36031625 1228 pmdTrack = new Int_t ***[fgkTotUM];
1229 pmdEdep = new Float_t ***[fgkTotUM];
1230 for (i=0; i<fgkTotUM; i++)
b2b13ee2 1231 {
36031625 1232 pmdTrack[i] = new Int_t **[fgkRow];
1233 pmdEdep[i] = new Float_t **[fgkRow];
b2b13ee2 1234 }
1235
36031625 1236 for (i = 0; i < fgkTotUM; i++)
b2b13ee2 1237 {
36031625 1238 for (j = 0; j < fgkRow; j++)
b2b13ee2 1239 {
36031625 1240 pmdTrack[i][j] = new Int_t *[fgkCol];
1241 pmdEdep[i][j] = new Float_t *[fgkCol];
b2b13ee2 1242 }
1243 }
1244
36031625 1245 for (i = 0; i < fgkTotUM; i++)
b2b13ee2 1246 {
36031625 1247 for (j = 0; j < fgkRow; j++)
b2b13ee2 1248 {
36031625 1249 for (k = 0; k < fgkCol; k++)
b2b13ee2 1250 {
36031625 1251 Int_t nn = fPRECounter[i][j][k];
b2b13ee2 1252 if(nn > 0)
1253 {
36031625 1254 pmdTrack[i][j][k] = new Int_t[nn];
1255 pmdEdep[i][j][k] = new Float_t[nn];
b2b13ee2 1256 }
1257 else
1258 {
1259 nn = 1;
36031625 1260 pmdTrack[i][j][k] = new Int_t[nn];
1261 pmdEdep[i][j][k] = new Float_t[nn];
5e2cc5bc 1262 }
36031625 1263 fPRECounter[i][j][k] = 0;
b2b13ee2 1264 }
1265 }
1266 }
1267
1268
ebd83c56 1269 Int_t nentries = fCell.GetEntries();
b2b13ee2 1270
1271 Int_t mtrackno, ism, ixp, iyp;
1272 Float_t edep;
1273
1274 for (i = 0; i < nentries; i++)
1275 {
ebd83c56 1276 AliPMDcell* cell = (AliPMDcell*)fCell.UncheckedAt(i);
b2b13ee2 1277
ebd83c56 1278 mtrackno = cell->GetTrackNumber();
1279 ism = cell->GetSMNumber();
1280 ixp = cell->GetX();
1281 iyp = cell->GetY();
1282 edep = cell->GetEdep();
36031625 1283 Int_t nn = fPRECounter[ism][ixp][iyp];
36031625 1284 pmdTrack[ism][ixp][iyp][nn] = (Int_t) mtrackno;
1285 pmdEdep[ism][ixp][iyp][nn] = edep;
1286 fPRECounter[ism][ixp][iyp]++;
b2b13ee2 1287 }
1288
1289 Int_t iz, il;
1290 Int_t im, ix, iy;
1291 Int_t nn;
1292
36031625 1293 for (im=0; im<fgkTotUM; im++)
b2b13ee2 1294 {
36031625 1295 for (ix=0; ix<fgkRow; ix++)
b2b13ee2 1296 {
36031625 1297 for (iy=0; iy<fgkCol; iy++)
b2b13ee2 1298 {
36031625 1299 nn = fPRECounter[im][ix][iy];
b2b13ee2 1300 if (nn > 1)
1301 {
1302 // This block handles if a cell is fired
1303 // many times by many tracks
75bb524e 1304 status1 = new Int_t[nn];
1305 status2 = new Int_t[nn];
1306 trnarray = new Int_t[nn];
b2b13ee2 1307 for (iz = 0; iz < nn; iz++)
1308 {
36031625 1309 status1[iz] = pmdTrack[im][ix][iy][iz];
b2b13ee2 1310 }
75bb524e 1311 TMath::Sort(nn,status1,status2,jsort);
36031625 1312 Int_t trackOld = -99999;
1313 Int_t track, trCount = 0;
b2b13ee2 1314 for (iz = 0; iz < nn; iz++)
1315 {
75bb524e 1316 track = status1[status2[iz]];
36031625 1317 if (trackOld != track)
b2b13ee2 1318 {
36031625 1319 trnarray[trCount] = track;
1320 trCount++;
5e2cc5bc 1321 }
36031625 1322 trackOld = track;
b2b13ee2 1323 }
6eb80a7c 1324 delete [] status1;
1325 delete [] status2;
36031625 1326 Float_t totEdp = 0.;
1327 trEdp = new Float_t[trCount];
1328 fracEdp = new Float_t[trCount];
1329 for (il = 0; il < trCount; il++)
b2b13ee2 1330 {
36031625 1331 trEdp[il] = 0.;
75bb524e 1332 track = trnarray[il];
b2b13ee2 1333 for (iz = 0; iz < nn; iz++)
1334 {
36031625 1335 if (track == pmdTrack[im][ix][iy][iz])
b2b13ee2 1336 {
36031625 1337 trEdp[il] += pmdEdep[im][ix][iy][iz];
b2b13ee2 1338 }
1339 }
36031625 1340 totEdp += trEdp[il];
b2b13ee2 1341 }
36031625 1342 Int_t ilOld = 0;
1343 Float_t fracOld = 0.;
b2b13ee2 1344
36031625 1345 for (il = 0; il < trCount; il++)
b2b13ee2 1346 {
36031625 1347 fracEdp[il] = trEdp[il]/totEdp;
1348 if (fracOld < fracEdp[il])
b2b13ee2 1349 {
36031625 1350 fracOld = fracEdp[il];
1351 ilOld = il;
b2b13ee2 1352 }
1353 }
36031625 1354 fPRETrackNo[im][ix][iy] = trnarray[ilOld];
6eb80a7c 1355 delete [] fracEdp;
1356 delete [] trEdp;
1357 delete [] trnarray;
b2b13ee2 1358 }
1359 else if (nn == 1)
1360 {
1361 // This only handles if a cell is fired
1362 // by only one track
1363
36031625 1364 fPRETrackNo[im][ix][iy] = pmdTrack[im][ix][iy][0];
b2b13ee2 1365
1366 }
1367 else if (nn ==0)
1368 {
1369 // This is if no cell is fired
36031625 1370 fPRETrackNo[im][ix][iy] = -999;
b2b13ee2 1371 }
1372 } // end of iy
1373 } // end of ix
1374 } // end of im
1375
1376 // Delete all the pointers
1377
36031625 1378 for (i = 0; i < fgkTotUM; i++)
b2b13ee2 1379 {
36031625 1380 for (j = 0; j < fgkRow; j++)
b2b13ee2 1381 {
36031625 1382 for (k = 0; k < fgkCol; k++)
b2b13ee2 1383 {
36031625 1384 delete [] pmdTrack[i][j][k];
1385 delete [] pmdEdep[i][j][k];
b2b13ee2 1386 }
1387 }
1388 }
1389
36031625 1390 for (i = 0; i < fgkTotUM; i++)
b2b13ee2 1391 {
36031625 1392 for (j = 0; j < fgkRow; j++)
b2b13ee2 1393 {
36031625 1394 delete [] pmdTrack[i][j];
1395 delete [] pmdEdep[i][j];
b2b13ee2 1396 }
1397 }
1398
36031625 1399 for (i = 0; i < fgkTotUM; i++)
b2b13ee2 1400 {
36031625 1401 delete [] pmdTrack[i];
1402 delete [] pmdEdep[i];
b2b13ee2 1403 }
6eb80a7c 1404 delete [] pmdTrack;
1405 delete [] pmdEdep;
b2b13ee2 1406 //
1407 // End of the cell id assignment
1408 //
1409}
7088dbe9 1410//____________________________________________________________________________
a7d110b8 1411void AliPMDDigitizer::MeV2ADC(Float_t mev, Float_t & adc) const
b2b13ee2 1412{
36031625 1413 // This converts the simulated edep to ADC according to the
1414 // Test Beam Data
7ba814f3 1415 //PS Test in May 2009
4d59e381 1416 // KeV - ADC conversion for 12bit ADC
7ba814f3 1417 // Data provided by Ajay and agreed by Dr. Viyogi
1418 // MPV data used for the fit and taken here
1419
1420 const Float_t kConstant = -9.22;
1421 const Float_t kErConstant = 0.602225;
1422 const Float_t kSlope = 79.1119;
1423 const Float_t kErSlope = 0.844047;
4d59e381 1424
b8e69f1f 1425 Float_t cons = gRandom->Gaus(kConstant,kErConstant);
1426 Float_t slop = gRandom->Gaus(kSlope,kErSlope);
4d59e381 1427
1428 Float_t adc12bit = slop*mev*0.001 + cons;
1429
7ba814f3 1430 if (adc12bit < 0.) adc12bit = 0.;
04a35d8d 1431
4d59e381 1432 if(adc12bit < 1600.0)
b8e69f1f 1433 {
1434 adc = (Float_t) adc12bit;
1435 }
4d59e381 1436 else if (adc12bit >= 1600.0)
b8e69f1f 1437 {
4d59e381 1438 adc = 1600.0;
b8e69f1f 1439 }
b2b13ee2 1440}
7088dbe9 1441//____________________________________________________________________________
920e13db 1442void AliPMDDigitizer::AddSDigit(Int_t trnumber, Int_t trpid, Int_t det,
1443 Int_t smnumber, Int_t irow, Int_t icol,
1444 Float_t adc)
b2b13ee2 1445{
36031625 1446 // Add SDigit
1447 //
ebd83c56 1448 if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
b2b13ee2 1449 TClonesArray &lsdigits = *fSDigits;
920e13db 1450 new(lsdigits[fNsdigit++]) AliPMDsdigit(trnumber,trpid,det,smnumber,irow,icol,adc);
b2b13ee2 1451}
7088dbe9 1452//____________________________________________________________________________
b2b13ee2 1453
920e13db 1454void AliPMDDigitizer::AddDigit(Int_t trnumber, Int_t trpid, Int_t det,
1455 Int_t smnumber, Int_t irow, Int_t icol,
1456 Float_t adc)
b2b13ee2 1457{
36031625 1458 // Add Digit
1459 //
ebd83c56 1460 if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
b2b13ee2 1461 TClonesArray &ldigits = *fDigits;
920e13db 1462 new(ldigits[fNdigit++]) AliPMDdigit(trnumber,trpid, det,smnumber,irow,icol,adc);
b2b13ee2 1463}
7088dbe9 1464//____________________________________________________________________________
b2b13ee2 1465
b2b13ee2 1466void AliPMDDigitizer::SetZPosition(Float_t zpos)
1467{
1468 fZPos = zpos;
1469}
7088dbe9 1470//____________________________________________________________________________
b2b13ee2 1471Float_t AliPMDDigitizer::GetZPosition() const
1472{
1473 return fZPos;
1474}
7088dbe9 1475//____________________________________________________________________________
b2b13ee2 1476
1477void AliPMDDigitizer::ResetCell()
1478{
36031625 1479 // clears the cell array and also the counter
1480 // for each cell
1481 //
4d59e381 1482 fCPVCell.Delete();
ebd83c56 1483 fCell.Delete();
36031625 1484 for (Int_t i = 0; i < fgkTotUM; i++)
b2b13ee2 1485 {
36031625 1486 for (Int_t j = 0; j < fgkRow; j++)
b2b13ee2 1487 {
36031625 1488 for (Int_t k = 0; k < fgkCol; k++)
b2b13ee2 1489 {
4d59e381 1490 fCPVCounter[i][j][k] = 0;
5e2cc5bc 1491 fPRECounter[i][j][k] = 0;
b2b13ee2 1492 }
1493 }
1494 }
1495}
7088dbe9 1496//____________________________________________________________________________
b2b13ee2 1497void AliPMDDigitizer::ResetSDigit()
1498{
36031625 1499 // Clears SDigits
b2b13ee2 1500 fNsdigit = 0;
ebd83c56 1501 if (fSDigits) fSDigits->Delete();
b2b13ee2 1502}
7088dbe9 1503//____________________________________________________________________________
b2b13ee2 1504void AliPMDDigitizer::ResetDigit()
1505{
36031625 1506 // Clears Digits
b2b13ee2 1507 fNdigit = 0;
ebd83c56 1508 if (fDigits) fDigits->Delete();
b2b13ee2 1509}
7088dbe9 1510//____________________________________________________________________________
b2b13ee2 1511
1512void AliPMDDigitizer::ResetCellADC()
1513{
7088dbe9 1514 // Clears individual cells edep and track number
36031625 1515 for (Int_t i = 0; i < fgkTotUM; i++)
b2b13ee2 1516 {
36031625 1517 for (Int_t j = 0; j < fgkRow; j++)
b2b13ee2 1518 {
36031625 1519 for (Int_t k = 0; k < fgkCol; k++)
b2b13ee2 1520 {
920e13db 1521 fCPV[i][j][k] = 0.;
1522 fPRE[i][j][k] = 0.;
1523 fCPVTrackNo[i][j][k] = 0;
1524 fPRETrackNo[i][j][k] = 0;
1525 fCPVTrackPid[i][j][k] = -1;
1526 fPRETrackPid[i][j][k] = -1;
b2b13ee2 1527 }
1528 }
1529 }
1530}
7088dbe9 1531//____________________________________________________________________________
b2b13ee2 1532
1533void AliPMDDigitizer::UnLoad(Option_t *option)
1534{
36031625 1535 // Unloads all the root files
1536 //
b2b13ee2 1537 const char *cS = strstr(option,"S");
1538 const char *cD = strstr(option,"D");
1539
1540 fRunLoader->UnloadgAlice();
1541 fRunLoader->UnloadHeader();
1542 fRunLoader->UnloadKinematics();
1543
1544 if (cS)
1545 {
36031625 1546 fPMDLoader->UnloadHits();
b2b13ee2 1547 }
1548 if (cD)
1549 {
36031625 1550 fPMDLoader->UnloadHits();
1551 fPMDLoader->UnloadSDigits();
b2b13ee2 1552 }
1553}
09a06455 1554
1555//----------------------------------------------------------------------
1556Float_t AliPMDDigitizer::Gain(Int_t det, Int_t smn, Int_t row, Int_t col) const
1557{
1558 // returns of the gain of the cell
1559 // Added this method by ZA
1560
1561 //cout<<" I am here in gain "<<fCalibData<< "smn,row, col "<<smn
1562 //<<" "<<row<<" "<<col<<endl;
1563
35535af7 1564 if(!fCalibGain) {
09a06455 1565 AliError("No calibration data loaded from CDB!!!");
50555ba1 1566 return 1;
09a06455 1567 }
1568
1569 Float_t GainFact;
35535af7 1570 GainFact = fCalibGain->GetGainFact(det,smn,row,col);
09a06455 1571 return GainFact;
1572}
1573//----------------------------------------------------------------------
35535af7 1574AliPMDCalibData* AliPMDDigitizer::GetCalibGain() const
09a06455 1575{
1576 // The run number will be centralized in AliCDBManager,
1577 // you don't need to set it here!
1578 // Added this method by ZA
0dd3d6f9 1579 // Cleaned up by Alberto
35535af7 1580 AliCDBEntry *entry = AliCDBManager::Instance()->Get("PMD/Calib/Gain");
09a06455 1581
0dd3d6f9 1582 if(!entry) AliFatal("Calibration object retrieval failed!");
09a06455 1583
1584 AliPMDCalibData *calibdata=0;
1585 if (entry) calibdata = (AliPMDCalibData*) entry->GetObject();
1586
0dd3d6f9 1587 if (!calibdata) AliFatal("No calibration data from calibration database !");
09a06455 1588
1589 return calibdata;
1590}
35535af7 1591//----------------------------------------------------------------------
1592AliPMDPedestal* AliPMDDigitizer::GetCalibPed() const
1593{
1594 // The run number will be centralized in AliCDBManager,
1595 // you don't need to set it here!
1596
1597 AliCDBEntry *entry = AliCDBManager::Instance()->Get("PMD/Calib/Ped");
1598
1599 if(!entry) AliFatal("Pedestal object retrieval failed!");
1600
1601 AliPMDPedestal *pedestal=0;
1602 if (entry) pedestal = (AliPMDPedestal*) entry->GetObject();
1603
1604 if (!pedestal) AliFatal("No pedestal data from calibration database !");
1605
1606 return pedestal;
1607}