Removing the tasks from the digitization (Ruben)
[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"
f21fc003 43#include "AliDigitizationInput.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//____________________________________________________________________________
f21fc003 131AliPMDDigitizer::AliPMDDigitizer(AliDigitizationInput* digInput):
132 AliDigitizer(digInput),
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;
414d1234 264 Int_t npmd = 0;
265 Int_t trackno = 0;
266 Int_t smnumber = 0;
267 Int_t trackpid = 0;
268 Int_t mtrackno = 0;
269 Int_t mtrackpid = 0;
270
271 Float_t xPos = 0., yPos = 0., zPos = 0.;
8599803d 272 Int_t xpad = -1, ypad = -1;
414d1234 273 Float_t edep = 0.;
b2b13ee2 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;
ffdc86b5 335
336 Int_t trnotemp = trackno; // Modified on 25th Nov 2009
4d59e381 337 if(ks==1||(imo = mparticle->GetFirstMother())<0 ){
338 vx = mparticle->Vx();
339 vy = mparticle->Vy();
340 vz = mparticle->Vz();
341
342 if(trackpid==kGamma||trackpid==11||trackpid==-11||
343 trackpid==kPi0)igstatus=1;
344 }
24c5571f 345
346
4d59e381 347 while(((imo = mparticle->GetFirstMother()) >= 0)&&
348 (ks = mparticle->GetStatusCode() <1) )
5f7ce82b 349 {
b208c6a3 350 mparticle = gAlice->GetMCApp()->Particle(imo);
4d59e381 351 trackpid = mparticle->GetPdgCode();
352 ks = mparticle->GetStatusCode();
b2b13ee2 353 vx = mparticle->Vx();
354 vy = mparticle->Vy();
355 vz = mparticle->Vz();
24c5571f 356
ffdc86b5 357 // Modified on 25th Nov 2009
358
359 trnotemp = trackno;
360 if(trackpid == 111)
361 {
362 trackno = trnotemp;
363 }
364 if(trackpid != 111)
365 {
366 trackno=imo;
5f7ce82b 367 }
ffdc86b5 368 // end of modification on 25th Nov 2009
369 }
24c5571f 370
5f7ce82b 371 if(trackpid==kGamma||trackpid==11||trackpid==-11||
372 trackpid==kPi0)igstatus=1;
373 mtrackpid=trackpid;
374 mtrackno=trackno;
375 trackpid=trackpidOld;
376 trackno=tracknoOld;
d934d257 377
24c5571f 378 //-----------------end of modification----------------
ae909aed 379 Float_t ptime = fPMDHit->GetTime()*1e6; // time in microsec
380 if (ptime < 0. || ptime > 1.2) continue;
381
36031625 382 xPos = fPMDHit->X();
383 yPos = fPMDHit->Y();
384 zPos = fPMDHit->Z();
ae909aed 385
c35544fe 386 edep = fPMDHit->GetEnergy();
387 Int_t vol1 = fPMDHit->GetVolume(1); // Column
388 Int_t vol2 = fPMDHit->GetVolume(2); // Row
24c5571f 389 Int_t vol7 = fPMDHit->GetVolume(4); // Serial Module No
f117e3aa 390
c35544fe 391
8599803d 392 // -----------------------------------------//
f117e3aa 393 // In new geometry after adding electronics //
8599803d 394 // For Super Module 1 & 2 //
8599803d 395 // nrow = 48, ncol = 96 //
f117e3aa 396 // For Super Module 3 & 4 //
397 // nrow = 96, ncol = 48 //
8599803d 398 // -----------------------------------------//
f117e3aa 399
24c5571f 400 if (vol7 < 24)
401 {
402 smnumber = vol7;
403 }
404 else
405 {
406 smnumber = vol7 - 24;
407 }
408 Int_t vol8 = smnumber/6 + 1; // fake supermodule
8599803d 409
f117e3aa 410 if (vol8 == 1 || vol8 == 2)
b2b13ee2 411 {
36031625 412 xpad = vol2;
413 ypad = vol1;
b2b13ee2 414 }
f117e3aa 415 else if (vol8 == 3 || vol8 == 4)
416 {
417 xpad = vol1;
418 ypad = vol2;
419 }
b2b13ee2 420
ecee2a1a 421 AliDebug(2,Form("Zposition = %f Edeposition = %f",zPos,edep));
c505f70b 422
423 if (vol7 < 24)
b2b13ee2 424 {
c505f70b 425 // PRE
426 fDetNo = 0;
b2b13ee2 427 }
c505f70b 428 else
b2b13ee2 429 {
c505f70b 430 // CPV
431 fDetNo = 1;
b2b13ee2 432 }
ffdc86b5 433
24c5571f 434 Int_t smn = smnumber;
8599803d 435 Int_t ixx = xpad - 1;
436 Int_t iyy = ypad - 1;
b2b13ee2 437 if (fDetNo == 0)
438 {
36031625 439 fPRE[smn][ixx][iyy] += edep;
440 fPRECounter[smn][ixx][iyy]++;
b2b13ee2 441
ebd83c56 442 AliPMDcell* cell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
ebd83c56 443 fCell.Add(cell);
b2b13ee2 444 }
445 else if(fDetNo == 1)
446 {
8599803d 447 fCPV[smn][ixx][iyy] += edep;
4d59e381 448 fCPVCounter[smn][ixx][iyy]++;
449 AliPMDcell* cpvcell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
450 fCPVCell.Add(cpvcell);
b2b13ee2 451 }
452 }
453 }
454 } // Track Loop ended
4d59e381 455 TrackAssignment2CPVCell();
b2b13ee2 456 TrackAssignment2Cell();
b2b13ee2 457 ResetCell();
458
459 Float_t deltaE = 0.;
460 Int_t detno = 0;
b2b13ee2 461 Int_t trno = -1;
920e13db 462 Int_t trpid = -99;
b2b13ee2 463
464 for (Int_t idet = 0; idet < 2; idet++)
465 {
36031625 466 for (Int_t ism = 0; ism < fgkTotUM; ism++)
b2b13ee2 467 {
36031625 468 for (Int_t jrow = 0; jrow < fgkRow; jrow++)
b2b13ee2 469 {
36031625 470 for (Int_t kcol = 0; kcol < fgkCol; kcol++)
b2b13ee2 471 {
b2b13ee2 472 if (idet == 0)
473 {
36031625 474 deltaE = fPRE[ism][jrow][kcol];
475 trno = fPRETrackNo[ism][jrow][kcol];
b2b13ee2 476 detno = 0;
477 }
478 else if (idet == 1)
479 {
480 deltaE = fCPV[ism][jrow][kcol];
481 trno = fCPVTrackNo[ism][jrow][kcol];
482 detno = 1;
483 }
484 if (deltaE > 0.)
485 {
8c24e5f2 486 // Natasha
487 TParticle *mparticle = gAlice->GetMCApp()->Particle(trno);
488 trpid = mparticle->GetPdgCode();
920e13db 489 AddSDigit(trno,trpid,detno,ism,jrow,kcol,deltaE);
b2b13ee2 490 }
491 }
492 }
ebd83c56 493 treeS->Fill();
b2b13ee2 494 ResetSDigit();
495 }
496 }
36031625 497 fPMDLoader->WriteSDigits("OVERWRITE");
b2b13ee2 498 ResetCellADC();
b2b13ee2 499}
7088dbe9 500//____________________________________________________________________________
b2b13ee2 501
502void AliPMDDigitizer::Hits2Digits(Int_t ievt)
503{
36031625 504 // This reads the PMD Hits tree and assigns the right track number
505 // to a cell and stores in the digits tree
506 //
507 const Int_t kPi0 = 111;
508 const Int_t kGamma = 22;
414d1234 509 Int_t npmd = 0;
510 Int_t trackno = 0;
511 Int_t smnumber = 0;
512 Int_t trackpid = 0;
513 Int_t mtrackno = 0;
514 Int_t mtrackpid = 0;
515
516 Float_t xPos = 0., yPos = 0., zPos = 0.;
8599803d 517 Int_t xpad = -1, ypad = -1;
414d1234 518 Float_t edep = 0.;
b2b13ee2 519 Float_t vx = -999.0, vy = -999.0, vz = -999.0;
520
ebd83c56 521 if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
b2b13ee2 522 ResetDigit();
523
5e2cc5bc 524 AliDebug(1,Form("Event Number = %d",ievt));
b2b13ee2 525 Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
ecee2a1a 526 AliDebug(1,Form("Number of Particles = %d", nparticles));
920e13db 527
b2b13ee2 528 fRunLoader->GetEvent(ievt);
b2b13ee2 529 // ------------------------------------------------------- //
530 // Pointer to specific detector hits.
531 // Get pointers to Alice detectors and Hits containers
532
b208c6a3 533 fPMD = (AliPMD*)gAlice->GetDetector("PMD");
36031625 534 fPMDLoader = fRunLoader->GetLoader("PMDLoader");
b2b13ee2 535
36031625 536 if (fPMDLoader == 0x0)
b2b13ee2 537 {
ecee2a1a 538 AliError("Can not find PMD or PMDLoader");
b2b13ee2 539 }
ebd83c56 540 TTree* treeH = fPMDLoader->TreeH();
541 Int_t ntracks = (Int_t) treeH->GetEntries();
ecee2a1a 542 AliDebug(1,Form("Number of Tracks in the TreeH = %d", ntracks));
36031625 543 fPMDLoader->LoadDigits("recreate");
ebd83c56 544 TTree* treeD = fPMDLoader->TreeD();
545 if (treeD == 0x0)
b2b13ee2 546 {
36031625 547 fPMDLoader->MakeTree("D");
ebd83c56 548 treeD = fPMDLoader->TreeD();
b2b13ee2 549 }
550 Int_t bufsize = 16000;
5e2cc5bc 551 treeD->Branch("PMDDigit", &fDigits, bufsize);
b2b13ee2 552
ebd83c56 553 TClonesArray* hits = 0;
554 if (fPMD) hits = fPMD->Hits();
b2b13ee2 555
556 // Start loop on tracks in the hits containers
557
5e2cc5bc 558 for (Int_t track=0; track<ntracks;track++)
b2b13ee2 559 {
cad3294f 560 gAlice->GetMCApp()->ResetHits();
ebd83c56 561 treeH->GetEvent(track);
b2b13ee2 562
5e2cc5bc 563 if (fPMD)
b2b13ee2 564 {
ebd83c56 565 npmd = hits->GetEntriesFast();
4d59e381 566 for (Int_t ipmd = 0; ipmd < npmd; ipmd++)
b2b13ee2 567 {
ebd83c56 568 fPMDHit = (AliPMDhit*) hits->UncheckedAt(ipmd);
36031625 569 trackno = fPMDHit->GetTrack();
b2b13ee2 570
571 // get kinematics of the particles
572
ebd83c56 573 TParticle* mparticle = gAlice->GetMCApp()->Particle(trackno);
574 trackpid = mparticle->GetPdgCode();
4d59e381 575 Int_t ks = mparticle->GetStatusCode();
b2b13ee2 576 Int_t imo;
36031625 577 Int_t tracknoOld=0, trackpidOld=0, statusOld = 0;
8599803d 578 if (mparticle->GetFirstMother() == -1)
579 {
36031625 580 tracknoOld = trackno;
581 trackpidOld = trackpid;
582 statusOld = -1;
8599803d 583 }
584
585 Int_t igstatus = 0;
ffdc86b5 586
587 Int_t trnotemp = trackno; // modified on 25th Nov 2009
4d59e381 588 if(ks==1||(imo = mparticle->GetFirstMother())<0 ){
589 vx = mparticle->Vx();
590 vy = mparticle->Vy();
591 vz = mparticle->Vz();
592
593 if(trackpid==kGamma||trackpid==11||trackpid==-11||trackpid==kPi0)
594 igstatus=1;
595 }
596
597
598 while(((imo = mparticle->GetFirstMother()) >= 0)&&
599 (ks = mparticle->GetStatusCode() <1) )
b2b13ee2 600 {
b208c6a3 601 mparticle = gAlice->GetMCApp()->Particle(imo);
4d59e381 602 trackpid = mparticle->GetPdgCode();
603 ks = mparticle->GetStatusCode();
b2b13ee2 604 vx = mparticle->Vx();
605 vy = mparticle->Vy();
606 vz = mparticle->Vz();
24c5571f 607
ffdc86b5 608 // Modified on 25th Nov 2009
609
610 trnotemp = trackno;
611 if(trackpid == 111)
612 {
613 trackno = trnotemp;
614 }
615 if(trackpid != 111)
616 {
617 trackno=imo;
d934d257 618 }
ffdc86b5 619 }
4d59e381 620
621 if(trackpid==kGamma||trackpid==11||trackpid==-11||trackpid==kPi0)
622 igstatus=1;
623 mtrackpid=trackpid;
624 mtrackno=trackno;
625 trackpid=trackpidOld;
626 trackno=tracknoOld;
ae909aed 627
628 Float_t ptime = fPMDHit->GetTime()*1e6;
629 if (ptime < 0. || ptime > 1.2) continue;
b2b13ee2 630
36031625 631 xPos = fPMDHit->X();
632 yPos = fPMDHit->Y();
633 zPos = fPMDHit->Z();
f117e3aa 634 edep = fPMDHit->GetEnergy();
c35544fe 635 Int_t vol1 = fPMDHit->GetVolume(1); // Column
636 Int_t vol2 = fPMDHit->GetVolume(2); // Row
24c5571f 637 Int_t vol7 = fPMDHit->GetVolume(4); // Serial Module No
8599803d 638
639 // -----------------------------------------//
f117e3aa 640 // In new geometry after adding electronics //
8599803d 641 // For Super Module 1 & 2 //
8599803d 642 // nrow = 48, ncol = 96 //
f117e3aa 643 // For Super Module 3 & 4 //
644 // nrow = 96, ncol = 48 //
8599803d 645 // -----------------------------------------//
b2b13ee2 646
24c5571f 647 if (vol7 < 24)
648 {
649 smnumber = vol7;
650 }
651 else
652 {
653 smnumber = vol7 - 24;
654 }
655 Int_t vol8 = smnumber/6 + 1; // fake supermodule
8599803d 656
f117e3aa 657 if (vol8 == 1 || vol8 == 2)
8599803d 658 {
36031625 659 xpad = vol2;
660 ypad = vol1;
b2b13ee2 661 }
f117e3aa 662 else if (vol8 == 3 || vol8 == 4)
663 {
664 xpad = vol1;
665 ypad = vol2;
666 }
b2b13ee2 667
66d5d07e 668 AliDebug(2,Form("ZPosition = %f Edeposition = %f",zPos,edep));
8599803d 669
c505f70b 670 if (vol7 < 24)
671 {
672 // PRE
673 fDetNo = 0;
674 }
675 else
b2b13ee2 676 {
b2b13ee2 677 fDetNo = 1;
678 }
414d1234 679
24c5571f 680 Int_t smn = smnumber;
8599803d 681 Int_t ixx = xpad - 1;
682 Int_t iyy = ypad - 1;
683 if (fDetNo == 0)
b2b13ee2 684 {
36031625 685 fPRE[smn][ixx][iyy] += edep;
686 fPRECounter[smn][ixx][iyy]++;
8599803d 687
ebd83c56 688 AliPMDcell* cell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
8599803d 689
ebd83c56 690 fCell.Add(cell);
b2b13ee2 691 }
8599803d 692 else if(fDetNo == 1)
b2b13ee2 693 {
8599803d 694 fCPV[smn][ixx][iyy] += edep;
4d59e381 695 fCPVCounter[smn][ixx][iyy]++;
696 AliPMDcell* cpvcell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
697 fCPVCell.Add(cpvcell);
b2b13ee2 698 }
699 }
700 }
701 } // Track Loop ended
4d59e381 702 TrackAssignment2CPVCell();
b2b13ee2 703 TrackAssignment2Cell();
704 ResetCell();
705
414d1234 706 Float_t gain1 = 1.;
707 Float_t adc = 0. ;
b2b13ee2 708 Float_t deltaE = 0.;
414d1234 709 Int_t detno = 0;
710 Int_t trno = 1;
711 Int_t trpid = -99;
920e13db 712
b2b13ee2 713 for (Int_t idet = 0; idet < 2; idet++)
5e2cc5bc 714 {
36031625 715 for (Int_t ism = 0; ism < fgkTotUM; ism++)
5e2cc5bc 716 {
36031625 717 for (Int_t jrow = 0; jrow < fgkRow; jrow++)
5e2cc5bc 718 {
36031625 719 for (Int_t kcol = 0; kcol < fgkCol; kcol++)
5e2cc5bc 720 {
b2b13ee2 721 if (idet == 0)
5e2cc5bc 722 {
723 deltaE = fPRE[ism][jrow][kcol];
36031625 724 trno = fPRETrackNo[ism][jrow][kcol];
920e13db 725 detno = 0;
5e2cc5bc 726 }
b2b13ee2 727 else if (idet == 1)
5e2cc5bc 728 {
729 deltaE = fCPV[ism][jrow][kcol];
8599803d 730 trno = fCPVTrackNo[ism][jrow][kcol];
920e13db 731 detno = 1;
5e2cc5bc 732 }
b2b13ee2 733 if (deltaE > 0.)
5e2cc5bc 734 {
7088dbe9 735 MeV2ADC(deltaE,adc);
5e2cc5bc 736
737 // To decalibrate the adc values
738 //
739 gain1 = Gain(idet,ism,jrow,kcol);
740 if (gain1 != 0.)
741 {
742 Int_t adcDecalib = (Int_t)(adc/gain1);
743 adc = (Float_t) adcDecalib;
744 }
745 else if(gain1 == 0.)
746 {
747 adc = 0.;
748 }
35535af7 749
750 // Pedestal Decalibration
751 Int_t pedmeanrms =
752 fCalibPed->GetPedMeanRms(idet,ism,jrow,kcol);
37edc588 753 Int_t pedrms1 = (Int_t) pedmeanrms%100;
35535af7 754 Float_t pedrms = (Float_t)pedrms1/10.;
755 Float_t pedmean =
756 (Float_t) (pedmeanrms - pedrms1)/1000.0;
35535af7 757 if (adc > 0.)
758 {
759 adc += (pedmean + 3.0*pedrms);
8c24e5f2 760 TParticle *mparticle
761 = gAlice->GetMCApp()->Particle(trno);
762 trpid = mparticle->GetPdgCode();
763
920e13db 764 AddDigit(trno,trpid,detno,ism,jrow,kcol,adc);
35535af7 765 }
5e2cc5bc 766 }
767 } // column loop
768 } // row loop
03965b13 769 treeD->Fill();
770 ResetDigit();
5e2cc5bc 771 } // supermodule loop
772 } // detector loop
03965b13 773
36031625 774 fPMDLoader->WriteDigits("OVERWRITE");
b2b13ee2 775 ResetCellADC();
920e13db 776
b2b13ee2 777}
7088dbe9 778//____________________________________________________________________________
b2b13ee2 779
780
781void AliPMDDigitizer::SDigits2Digits(Int_t ievt)
782{
36031625 783 // This reads the PMD sdigits tree and converts energy deposition
784 // in a cell to ADC and stores in the digits tree
785 //
ecee2a1a 786
b2b13ee2 787 fRunLoader->GetEvent(ievt);
788
ebd83c56 789 TTree* treeS = fPMDLoader->TreeS();
b2b13ee2 790 AliPMDsdigit *pmdsdigit;
ebd83c56 791 TBranch *branch = treeS->GetBranch("PMDSDigit");
ecee2a1a 792 if(!branch)
793 {
794 AliError("PMD Sdigit branch does not exist");
795 return;
796 }
ebd83c56 797 if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
b2b13ee2 798 branch->SetAddress(&fSDigits);
799
ebd83c56 800 TTree* treeD = fPMDLoader->TreeD();
801 if (treeD == 0x0)
b2b13ee2 802 {
36031625 803 fPMDLoader->MakeTree("D");
ebd83c56 804 treeD = fPMDLoader->TreeD();
b2b13ee2 805 }
806 Int_t bufsize = 16000;
ebd83c56 807 if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
5e2cc5bc 808 treeD->Branch("PMDDigit", &fDigits, bufsize);
b2b13ee2 809
414d1234 810 Int_t trno = 1, trpid = 0, det = 0, smn = 0;
811 Int_t irow = 0, icol = 0;
812 Float_t edep = 0., adc = 0.;
b2b13ee2 813
ebd83c56 814 Int_t nmodules = (Int_t) treeS->GetEntries();
ecee2a1a 815 AliDebug(1,Form("Number of modules = %d",nmodules));
b2b13ee2 816
817 for (Int_t imodule = 0; imodule < nmodules; imodule++)
818 {
5e2cc5bc 819 treeS->GetEntry(imodule);
b2b13ee2 820 Int_t nentries = fSDigits->GetLast();
ecee2a1a 821 AliDebug(2,Form("Number of entries per module = %d",nentries+1));
b2b13ee2 822 for (Int_t ient = 0; ient < nentries+1; ient++)
823 {
824 pmdsdigit = (AliPMDsdigit*)fSDigits->UncheckedAt(ient);
825 trno = pmdsdigit->GetTrackNumber();
920e13db 826 trpid = pmdsdigit->GetTrackPid();
b2b13ee2 827 det = pmdsdigit->GetDetector();
828 smn = pmdsdigit->GetSMNumber();
7088dbe9 829 irow = pmdsdigit->GetRow();
830 icol = pmdsdigit->GetColumn();
b2b13ee2 831 edep = pmdsdigit->GetCellEdep();
832
833 MeV2ADC(edep,adc);
5e2cc5bc 834
835 // To decalibrte the adc values
836 //
837 Float_t gain1 = Gain(det,smn,irow,icol);
838 if (gain1 != 0.)
839 {
04a35d8d 840 Int_t adcDecalib = (Int_t)(adc/gain1);
841 adc = (Float_t) adcDecalib;
5e2cc5bc 842 }
843 else if(gain1 == 0.)
844 {
845 adc = 0.;
846 }
35535af7 847 // Pedestal Decalibration
848 Int_t pedmeanrms = fCalibPed->GetPedMeanRms(det,smn,irow,icol);
37edc588 849 Int_t pedrms1 = (Int_t) pedmeanrms%100;
35535af7 850 Float_t pedrms = (Float_t)pedrms1/10.;
851 Float_t pedmean = (Float_t) (pedmeanrms - pedrms1)/1000.0;
35535af7 852 if(adc > 0.)
853 {
854 adc += (pedmean + 3.0*pedrms);
920e13db 855 AddDigit(trno,trpid,det,smn,irow,icol,adc);
35535af7 856 }
5e2cc5bc 857
b2b13ee2 858 }
ebd83c56 859 treeD->Fill();
b2b13ee2 860 ResetDigit();
861 }
36031625 862 fPMDLoader->WriteDigits("OVERWRITE");
ecee2a1a 863
b2b13ee2 864}
7088dbe9 865//____________________________________________________________________________
f21fc003 866void AliPMDDigitizer::Digitize(Option_t *option)
7088dbe9 867{
868 // Does the event merging and digitization
7088dbe9 869 const char *cdeb = strstr(option,"deb");
870 if(cdeb)
871 {
ecee2a1a 872 AliDebug(100," *** PMD Exec is called ***");
7088dbe9 873 }
ecee2a1a 874
f21fc003 875 Int_t ninputs = fDigInput->GetNinputs();
ecee2a1a 876 AliDebug(1,Form("Number of files to be processed = %d",ninputs));
7088dbe9 877 ResetCellADC();
878
879 for (Int_t i = 0; i < ninputs; i++)
880 {
f21fc003 881 Int_t troffset = fDigInput->GetMask(i);
7088dbe9 882 MergeSDigits(i, troffset);
883 }
b2b13ee2 884
f21fc003 885 fRunLoader = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
7088dbe9 886 fPMD = (AliPMD*)gAlice->GetDetector("PMD");
887 fPMDLoader = fRunLoader->GetLoader("PMDLoader");
888 if (fPMDLoader == 0x0)
889 {
ecee2a1a 890 AliError("Can not find PMD or PMDLoader");
7088dbe9 891 }
06f2c79d 892 fPMDLoader->LoadDigits("update");
ebd83c56 893 TTree* treeD = fPMDLoader->TreeD();
894 if (treeD == 0x0)
7088dbe9 895 {
896 fPMDLoader->MakeTree("D");
ebd83c56 897 treeD = fPMDLoader->TreeD();
7088dbe9 898 }
899 Int_t bufsize = 16000;
ebd83c56 900 if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
5e2cc5bc 901 treeD->Branch("PMDDigit", &fDigits, bufsize);
7088dbe9 902
414d1234 903 Float_t adc = 0.;
7088dbe9 904 Float_t deltaE = 0.;
920e13db 905 Int_t detno = 0;
906 Int_t trno = 1;
907 Int_t trpid = -99;
7088dbe9 908
909 for (Int_t idet = 0; idet < 2; idet++)
910 {
911 for (Int_t ism = 0; ism < fgkTotUM; ism++)
912 {
913 for (Int_t jrow = 0; jrow < fgkRow; jrow++)
914 {
915 for (Int_t kcol = 0; kcol < fgkCol; kcol++)
916 {
917 if (idet == 0)
918 {
919 deltaE = fPRE[ism][jrow][kcol];
920 trno = fPRETrackNo[ism][jrow][kcol];
920e13db 921 trpid = fPRETrackPid[ism][jrow][kcol];
922 detno = 0;
7088dbe9 923 }
924 else if (idet == 1)
925 {
926 deltaE = fCPV[ism][jrow][kcol];
927 trno = fCPVTrackNo[ism][jrow][kcol];
920e13db 928 trpid = fCPVTrackPid[ism][jrow][kcol];
929 detno = 1;
7088dbe9 930 }
931 if (deltaE > 0.)
932 {
933 MeV2ADC(deltaE,adc);
35535af7 934
5f566ada 935 //
35535af7 936 // Gain decalibration
5f566ada 937 //
938 Float_t gain1 = Gain(idet,ism,jrow,kcol);
35535af7 939
5f566ada 940 if (gain1 != 0.)
941 {
942 Int_t adcDecalib = (Int_t)(adc/gain1);
943 adc = (Float_t) adcDecalib;
944 }
945 else if(gain1 == 0.)
946 {
947 adc = 0.;
948 }
35535af7 949 // Pedestal Decalibration
950 Int_t pedmeanrms =
951 fCalibPed->GetPedMeanRms(idet,ism,jrow,kcol);
37edc588 952 Int_t pedrms1 = (Int_t) pedmeanrms%100;
35535af7 953 Float_t pedrms = (Float_t)pedrms1/10.;
954 Float_t pedmean =
955 (Float_t) (pedmeanrms - pedrms1)/1000.0;
35535af7 956 if (adc > 0.)
957 {
958 adc += (pedmean + 3.0*pedrms);
920e13db 959 AddDigit(trno,trpid,detno,ism,jrow,kcol,adc);
35535af7 960 }
5f566ada 961
7088dbe9 962 }
963 } // column loop
964 } // row loop
ebd83c56 965 treeD->Fill();
7088dbe9 966 ResetDigit();
967 } // supermodule loop
968 } // detector loop
5e2cc5bc 969 fPMDLoader->WriteDigits("OVERWRITE");
06f2c79d 970 fPMDLoader->UnloadDigits();
971 ResetCellADC();
7088dbe9 972}
973//____________________________________________________________________________
4d59e381 974void AliPMDDigitizer::TrackAssignment2CPVCell()
975{
976 // This block assigns the cell id when there are
977 // multiple tracks in a cell according to the
978 // energy deposition
979 // This method added by Ajay
980 Bool_t jsort = false;
981
414d1234 982 Int_t i = 0, j = 0, k = 0;
4d59e381 983
984 Int_t *status1;
985 Int_t *status2;
986 Int_t *trnarray;
987 Float_t *fracEdp;
988 Float_t *trEdp;
989
990 Int_t ****cpvTrack;
991 Float_t ****cpvEdep;
992
993 cpvTrack = new Int_t ***[fgkTotUM];
994 cpvEdep = new Float_t ***[fgkTotUM];
995 for (i=0; i<fgkTotUM; i++)
996 {
997 cpvTrack[i] = new Int_t **[fgkRow];
998 cpvEdep[i] = new Float_t **[fgkRow];
999 }
1000
1001 for (i = 0; i < fgkTotUM; i++)
1002 {
1003 for (j = 0; j < fgkRow; j++)
1004 {
1005 cpvTrack[i][j] = new Int_t *[fgkCol];
1006 cpvEdep[i][j] = new Float_t *[fgkCol];
1007 }
1008 }
1009 for (i = 0; i < fgkTotUM; i++)
1010 {
1011 for (j = 0; j < fgkRow; j++)
1012 {
1013 for (k = 0; k < fgkCol; k++)
1014 {
1015 Int_t nn = fCPVCounter[i][j][k];
1016 if(nn > 0)
1017 {
1018 cpvTrack[i][j][k] = new Int_t[nn];
1019 cpvEdep[i][j][k] = new Float_t[nn];
1020 }
1021 else
1022 {
1023 nn = 1;
1024 cpvTrack[i][j][k] = new Int_t[nn];
1025 cpvEdep[i][j][k] = new Float_t[nn];
1026 }
1027 fCPVCounter[i][j][k] = 0;
1028 }
1029 }
1030 }
1031
1032
1033 Int_t nentries = fCPVCell.GetEntries();
1034
414d1234 1035 Int_t mtrackno = 0, ism = 0, ixp = 0, iyp = 0;
1036 Float_t edep = 0.;
4d59e381 1037 for (i = 0; i < nentries; i++)
1038 {
1039 AliPMDcell* cpvcell = (AliPMDcell*)fCPVCell.UncheckedAt(i);
1040
1041 mtrackno = cpvcell->GetTrackNumber();
1042 ism = cpvcell->GetSMNumber();
1043 ixp = cpvcell->GetX();
1044 iyp = cpvcell->GetY();
1045 edep = cpvcell->GetEdep();
1046 Int_t nn = fCPVCounter[ism][ixp][iyp];
1047 cpvTrack[ism][ixp][iyp][nn] = (Int_t) mtrackno;
1048 cpvEdep[ism][ixp][iyp][nn] = edep;
1049 fCPVCounter[ism][ixp][iyp]++;
1050 }
1051
414d1234 1052 Int_t iz = 0, il = 0;
1053 Int_t im = 0, ix = 0, iy = 0;
1054 Int_t nn = 0;
4d59e381 1055 for (im=0; im<fgkTotUM; im++)
1056 {
1057 for (ix=0; ix<fgkRow; ix++)
1058 {
1059 for (iy=0; iy<fgkCol; iy++)
1060 {
1061 nn = fCPVCounter[im][ix][iy];
1062 if (nn > 1)
1063 {
1064 // This block handles if a cell is fired
1065 // many times by many tracks
1066 status1 = new Int_t[nn];
18418801 1067 status2 = new Int_t[2*nn];
4d59e381 1068 trnarray = new Int_t[nn];
1069 for (iz = 0; iz < nn; iz++)
1070 {
1071 status1[iz] = cpvTrack[im][ix][iy][iz];
1072 }
1073 TMath::Sort(nn,status1,status2,jsort);
1074 Int_t trackOld = -99999;
1075 Int_t track, trCount = 0;
1076 for (iz = 0; iz < nn; iz++)
1077 {
1078 track = status1[status2[iz]];
1079 if (trackOld != track)
1080 {
1081 trnarray[trCount] = track;
1082 trCount++;
1083 }
1084 trackOld = track;
1085 }
1086 delete [] status1;
1087 delete [] status2;
1088 Float_t totEdp = 0.;
1089 trEdp = new Float_t[trCount];
1090 fracEdp = new Float_t[trCount];
1091 for (il = 0; il < trCount; il++)
1092 {
1093 trEdp[il] = 0.;
1094 track = trnarray[il];
1095 for (iz = 0; iz < nn; iz++)
1096 {
1097 if (track == cpvTrack[im][ix][iy][iz])
1098 {
1099 trEdp[il] += cpvEdep[im][ix][iy][iz];
1100 }
1101 }
1102 totEdp += trEdp[il];
1103 }
1104 Int_t ilOld = 0;
1105 Float_t fracOld = 0.;
1106
1107 for (il = 0; il < trCount; il++)
1108 {
1109 fracEdp[il] = trEdp[il]/totEdp;
1110 if (fracOld < fracEdp[il])
1111 {
1112 fracOld = fracEdp[il];
1113 ilOld = il;
1114 }
1115 }
1116 fCPVTrackNo[im][ix][iy] = trnarray[ilOld];
1117 delete [] fracEdp;
1118 delete [] trEdp;
1119 delete [] trnarray;
1120 }
1121 else if (nn == 1)
1122 {
1123 // This only handles if a cell is fired
1124 // by only one track
1125
1126 fCPVTrackNo[im][ix][iy] = cpvTrack[im][ix][iy][0];
1127
1128 }
1129 else if (nn ==0)
1130 {
1131 // This is if no cell is fired
1132 fCPVTrackNo[im][ix][iy] = -999;
1133 }
1134 } // end of iy
1135 } // end of ix
1136 } // end of im
1137
1138 // Delete all the pointers
1139
1140 for (i = 0; i < fgkTotUM; i++)
1141 {
1142 for (j = 0; j < fgkRow; j++)
1143 {
1144 for (k = 0; k < fgkCol; k++)
1145 {
1146 delete []cpvTrack[i][j][k];
1147 delete []cpvEdep[i][j][k];
1148 }
1149 }
1150 }
1151
1152 for (i = 0; i < fgkTotUM; i++)
1153 {
1154 for (j = 0; j < fgkRow; j++)
1155 {
1156 delete [] cpvTrack[i][j];
1157 delete [] cpvEdep[i][j];
1158 }
1159 }
1160
1161 for (i = 0; i < fgkTotUM; i++)
1162 {
1163 delete [] cpvTrack[i];
1164 delete [] cpvEdep[i];
1165 }
1166 delete [] cpvTrack;
1167 delete [] cpvEdep;
1168
1169 //
1170 // End of the cell id assignment
1171 //
1172}
1173//____________________________________________________________________________
7088dbe9 1174
1175void AliPMDDigitizer::MergeSDigits(Int_t filenumber, Int_t troffset)
1176{
1177 // merging sdigits
f21fc003 1178 fRunLoader = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(filenumber));
7088dbe9 1179 fPMDLoader = fRunLoader->GetLoader("PMDLoader");
1180 fPMDLoader->LoadSDigits("read");
ebd83c56 1181 TTree* treeS = fPMDLoader->TreeS();
7088dbe9 1182 AliPMDsdigit *pmdsdigit;
ebd83c56 1183 TBranch *branch = treeS->GetBranch("PMDSDigit");
1184 if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
7088dbe9 1185 branch->SetAddress(&fSDigits);
1186
414d1234 1187 Int_t itrackno = 1, itrackpid = 0, idet = 0, ism = 0;
1188 Int_t ixp = 0, iyp = 0;
1189 Float_t edep = 0.;
ebd83c56 1190 Int_t nmodules = (Int_t) treeS->GetEntries();
ecee2a1a 1191 AliDebug(1,Form("Number of Modules in the treeS = %d",nmodules));
1192 AliDebug(1,Form("Track Offset = %d",troffset));
7088dbe9 1193 for (Int_t imodule = 0; imodule < nmodules; imodule++)
1194 {
5e2cc5bc 1195 treeS->GetEntry(imodule);
7088dbe9 1196 Int_t nentries = fSDigits->GetLast();
ecee2a1a 1197 AliDebug(2,Form("Number of Entries per Module = %d",nentries));
7088dbe9 1198 for (Int_t ient = 0; ient < nentries+1; ient++)
1199 {
1200 pmdsdigit = (AliPMDsdigit*)fSDigits->UncheckedAt(ient);
1201 itrackno = pmdsdigit->GetTrackNumber();
920e13db 1202 itrackpid = pmdsdigit->GetTrackPid();
7088dbe9 1203 idet = pmdsdigit->GetDetector();
1204 ism = pmdsdigit->GetSMNumber();
1205 ixp = pmdsdigit->GetRow();
1206 iyp = pmdsdigit->GetColumn();
1207 edep = pmdsdigit->GetCellEdep();
7088dbe9 1208 if (idet == 0)
1209 {
1210 if (fPRE[ism][ixp][iyp] < edep)
1211 {
1212 fPRETrackNo[ism][ixp][iyp] = troffset + itrackno;
920e13db 1213 fPRETrackPid[ism][ixp][iyp] = itrackpid;
7088dbe9 1214 }
1215 fPRE[ism][ixp][iyp] += edep;
1216 }
1217 else if (idet == 1)
1218 {
1219 if (fCPV[ism][ixp][iyp] < edep)
1220 {
1221 fCPVTrackNo[ism][ixp][iyp] = troffset + itrackno;
920e13db 1222 fCPVTrackPid[ism][ixp][iyp] = itrackpid;
7088dbe9 1223 }
1224 fCPV[ism][ixp][iyp] += edep;
1225 }
1226 }
1227 }
1228
1229}
1230// ----------------------------------------------------------------------
b2b13ee2 1231void AliPMDDigitizer::TrackAssignment2Cell()
1232{
b2b13ee2 1233 //
75bb524e 1234 // This block assigns the cell id when there are
b2b13ee2 1235 // multiple tracks in a cell according to the
1236 // energy deposition
1237 //
75bb524e 1238 Bool_t jsort = false;
8599803d 1239
414d1234 1240 Int_t i = 0, j = 0, k = 0;
b2b13ee2 1241
4d59e381 1242 Int_t *status1;
1243 Int_t *status2;
1244 Int_t *trnarray;
36031625 1245 Float_t *fracEdp;
1246 Float_t *trEdp;
4d59e381 1247
36031625 1248 Int_t ****pmdTrack;
1249 Float_t ****pmdEdep;
b2b13ee2 1250
36031625 1251 pmdTrack = new Int_t ***[fgkTotUM];
1252 pmdEdep = new Float_t ***[fgkTotUM];
1253 for (i=0; i<fgkTotUM; i++)
b2b13ee2 1254 {
36031625 1255 pmdTrack[i] = new Int_t **[fgkRow];
1256 pmdEdep[i] = new Float_t **[fgkRow];
b2b13ee2 1257 }
1258
36031625 1259 for (i = 0; i < fgkTotUM; i++)
b2b13ee2 1260 {
36031625 1261 for (j = 0; j < fgkRow; j++)
b2b13ee2 1262 {
36031625 1263 pmdTrack[i][j] = new Int_t *[fgkCol];
1264 pmdEdep[i][j] = new Float_t *[fgkCol];
b2b13ee2 1265 }
1266 }
1267
36031625 1268 for (i = 0; i < fgkTotUM; i++)
b2b13ee2 1269 {
36031625 1270 for (j = 0; j < fgkRow; j++)
b2b13ee2 1271 {
36031625 1272 for (k = 0; k < fgkCol; k++)
b2b13ee2 1273 {
36031625 1274 Int_t nn = fPRECounter[i][j][k];
b2b13ee2 1275 if(nn > 0)
1276 {
36031625 1277 pmdTrack[i][j][k] = new Int_t[nn];
1278 pmdEdep[i][j][k] = new Float_t[nn];
b2b13ee2 1279 }
1280 else
1281 {
1282 nn = 1;
36031625 1283 pmdTrack[i][j][k] = new Int_t[nn];
1284 pmdEdep[i][j][k] = new Float_t[nn];
5e2cc5bc 1285 }
36031625 1286 fPRECounter[i][j][k] = 0;
b2b13ee2 1287 }
1288 }
1289 }
1290
1291
ebd83c56 1292 Int_t nentries = fCell.GetEntries();
b2b13ee2 1293
1294 Int_t mtrackno, ism, ixp, iyp;
1295 Float_t edep;
1296
1297 for (i = 0; i < nentries; i++)
1298 {
ebd83c56 1299 AliPMDcell* cell = (AliPMDcell*)fCell.UncheckedAt(i);
b2b13ee2 1300
ebd83c56 1301 mtrackno = cell->GetTrackNumber();
1302 ism = cell->GetSMNumber();
1303 ixp = cell->GetX();
1304 iyp = cell->GetY();
1305 edep = cell->GetEdep();
36031625 1306 Int_t nn = fPRECounter[ism][ixp][iyp];
36031625 1307 pmdTrack[ism][ixp][iyp][nn] = (Int_t) mtrackno;
1308 pmdEdep[ism][ixp][iyp][nn] = edep;
1309 fPRECounter[ism][ixp][iyp]++;
b2b13ee2 1310 }
1311
414d1234 1312 Int_t iz = 0, il = 0;
1313 Int_t im = 0, ix = 0, iy = 0;
1314 Int_t nn = 0;
b2b13ee2 1315
36031625 1316 for (im=0; im<fgkTotUM; im++)
b2b13ee2 1317 {
36031625 1318 for (ix=0; ix<fgkRow; ix++)
b2b13ee2 1319 {
36031625 1320 for (iy=0; iy<fgkCol; iy++)
b2b13ee2 1321 {
36031625 1322 nn = fPRECounter[im][ix][iy];
b2b13ee2 1323 if (nn > 1)
1324 {
1325 // This block handles if a cell is fired
1326 // many times by many tracks
75bb524e 1327 status1 = new Int_t[nn];
18418801 1328 status2 = new Int_t[2*nn];
75bb524e 1329 trnarray = new Int_t[nn];
b2b13ee2 1330 for (iz = 0; iz < nn; iz++)
1331 {
36031625 1332 status1[iz] = pmdTrack[im][ix][iy][iz];
b2b13ee2 1333 }
75bb524e 1334 TMath::Sort(nn,status1,status2,jsort);
36031625 1335 Int_t trackOld = -99999;
1336 Int_t track, trCount = 0;
b2b13ee2 1337 for (iz = 0; iz < nn; iz++)
1338 {
75bb524e 1339 track = status1[status2[iz]];
36031625 1340 if (trackOld != track)
b2b13ee2 1341 {
36031625 1342 trnarray[trCount] = track;
1343 trCount++;
5e2cc5bc 1344 }
36031625 1345 trackOld = track;
b2b13ee2 1346 }
6eb80a7c 1347 delete [] status1;
1348 delete [] status2;
36031625 1349 Float_t totEdp = 0.;
1350 trEdp = new Float_t[trCount];
1351 fracEdp = new Float_t[trCount];
1352 for (il = 0; il < trCount; il++)
b2b13ee2 1353 {
36031625 1354 trEdp[il] = 0.;
75bb524e 1355 track = trnarray[il];
b2b13ee2 1356 for (iz = 0; iz < nn; iz++)
1357 {
36031625 1358 if (track == pmdTrack[im][ix][iy][iz])
b2b13ee2 1359 {
36031625 1360 trEdp[il] += pmdEdep[im][ix][iy][iz];
b2b13ee2 1361 }
1362 }
36031625 1363 totEdp += trEdp[il];
b2b13ee2 1364 }
36031625 1365 Int_t ilOld = 0;
1366 Float_t fracOld = 0.;
b2b13ee2 1367
36031625 1368 for (il = 0; il < trCount; il++)
b2b13ee2 1369 {
36031625 1370 fracEdp[il] = trEdp[il]/totEdp;
1371 if (fracOld < fracEdp[il])
b2b13ee2 1372 {
36031625 1373 fracOld = fracEdp[il];
1374 ilOld = il;
b2b13ee2 1375 }
1376 }
36031625 1377 fPRETrackNo[im][ix][iy] = trnarray[ilOld];
6eb80a7c 1378 delete [] fracEdp;
1379 delete [] trEdp;
1380 delete [] trnarray;
b2b13ee2 1381 }
1382 else if (nn == 1)
1383 {
1384 // This only handles if a cell is fired
1385 // by only one track
1386
36031625 1387 fPRETrackNo[im][ix][iy] = pmdTrack[im][ix][iy][0];
b2b13ee2 1388
1389 }
1390 else if (nn ==0)
1391 {
1392 // This is if no cell is fired
36031625 1393 fPRETrackNo[im][ix][iy] = -999;
b2b13ee2 1394 }
1395 } // end of iy
1396 } // end of ix
1397 } // end of im
1398
1399 // Delete all the pointers
1400
36031625 1401 for (i = 0; i < fgkTotUM; i++)
b2b13ee2 1402 {
36031625 1403 for (j = 0; j < fgkRow; j++)
b2b13ee2 1404 {
36031625 1405 for (k = 0; k < fgkCol; k++)
b2b13ee2 1406 {
36031625 1407 delete [] pmdTrack[i][j][k];
1408 delete [] pmdEdep[i][j][k];
b2b13ee2 1409 }
1410 }
1411 }
1412
36031625 1413 for (i = 0; i < fgkTotUM; i++)
b2b13ee2 1414 {
36031625 1415 for (j = 0; j < fgkRow; j++)
b2b13ee2 1416 {
36031625 1417 delete [] pmdTrack[i][j];
1418 delete [] pmdEdep[i][j];
b2b13ee2 1419 }
1420 }
1421
36031625 1422 for (i = 0; i < fgkTotUM; i++)
b2b13ee2 1423 {
36031625 1424 delete [] pmdTrack[i];
1425 delete [] pmdEdep[i];
b2b13ee2 1426 }
6eb80a7c 1427 delete [] pmdTrack;
1428 delete [] pmdEdep;
b2b13ee2 1429 //
1430 // End of the cell id assignment
1431 //
1432}
7088dbe9 1433//____________________________________________________________________________
a7d110b8 1434void AliPMDDigitizer::MeV2ADC(Float_t mev, Float_t & adc) const
b2b13ee2 1435{
36031625 1436 // This converts the simulated edep to ADC according to the
1437 // Test Beam Data
ab67948f 1438 // PS Test in June 2010, Voltage @ 1300 V
4d59e381 1439 // KeV - ADC conversion for 12bit ADC
7ba814f3 1440 // MPV data used for the fit and taken here
1441
414d1234 1442 // constants are from Test Beam 2010
4d59e381 1443
414d1234 1444 const Float_t kConstant = 0.612796;
1445 const Float_t kSlope = 130.158;
4d59e381 1446
414d1234 1447 Float_t adc12bit = kSlope*mev*0.001 + kConstant;
7ba814f3 1448 if (adc12bit < 0.) adc12bit = 0.;
414d1234 1449
1450 //Introducing Readout Resolution for ALICE-PMD
1451
1452 Float_t sigrr = 0.605016 - 0.000273*adc12bit + 6.54e-8*adc12bit*adc12bit;
1453 Float_t adcwithrr = gRandom->Gaus(adc12bit,sigrr);
1454
f564dd5a 1455 if(adcwithrr < 0.)
b8e69f1f 1456 {
f564dd5a 1457 adc = 0.;
1458 }
1459 else if(adcwithrr >= 0. && adcwithrr < 1600.0)
1460 {
1461 adc = adcwithrr;
b8e69f1f 1462 }
414d1234 1463 else if (adcwithrr >= 1600.0)
b8e69f1f 1464 {
4d59e381 1465 adc = 1600.0;
b8e69f1f 1466 }
414d1234 1467
b2b13ee2 1468}
7088dbe9 1469//____________________________________________________________________________
920e13db 1470void AliPMDDigitizer::AddSDigit(Int_t trnumber, Int_t trpid, Int_t det,
1471 Int_t smnumber, Int_t irow, Int_t icol,
1472 Float_t adc)
b2b13ee2 1473{
36031625 1474 // Add SDigit
1475 //
ebd83c56 1476 if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
b2b13ee2 1477 TClonesArray &lsdigits = *fSDigits;
920e13db 1478 new(lsdigits[fNsdigit++]) AliPMDsdigit(trnumber,trpid,det,smnumber,irow,icol,adc);
b2b13ee2 1479}
7088dbe9 1480//____________________________________________________________________________
b2b13ee2 1481
920e13db 1482void AliPMDDigitizer::AddDigit(Int_t trnumber, Int_t trpid, Int_t det,
1483 Int_t smnumber, Int_t irow, Int_t icol,
1484 Float_t adc)
b2b13ee2 1485{
36031625 1486 // Add Digit
1487 //
ebd83c56 1488 if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
b2b13ee2 1489 TClonesArray &ldigits = *fDigits;
920e13db 1490 new(ldigits[fNdigit++]) AliPMDdigit(trnumber,trpid, det,smnumber,irow,icol,adc);
b2b13ee2 1491}
7088dbe9 1492//____________________________________________________________________________
b2b13ee2 1493
b2b13ee2 1494void AliPMDDigitizer::SetZPosition(Float_t zpos)
1495{
1496 fZPos = zpos;
1497}
7088dbe9 1498//____________________________________________________________________________
b2b13ee2 1499Float_t AliPMDDigitizer::GetZPosition() const
1500{
1501 return fZPos;
1502}
7088dbe9 1503//____________________________________________________________________________
b2b13ee2 1504
1505void AliPMDDigitizer::ResetCell()
1506{
36031625 1507 // clears the cell array and also the counter
1508 // for each cell
1509 //
4d59e381 1510 fCPVCell.Delete();
ebd83c56 1511 fCell.Delete();
36031625 1512 for (Int_t i = 0; i < fgkTotUM; i++)
b2b13ee2 1513 {
36031625 1514 for (Int_t j = 0; j < fgkRow; j++)
b2b13ee2 1515 {
36031625 1516 for (Int_t k = 0; k < fgkCol; k++)
b2b13ee2 1517 {
4d59e381 1518 fCPVCounter[i][j][k] = 0;
5e2cc5bc 1519 fPRECounter[i][j][k] = 0;
b2b13ee2 1520 }
1521 }
1522 }
1523}
7088dbe9 1524//____________________________________________________________________________
b2b13ee2 1525void AliPMDDigitizer::ResetSDigit()
1526{
36031625 1527 // Clears SDigits
b2b13ee2 1528 fNsdigit = 0;
ebd83c56 1529 if (fSDigits) fSDigits->Delete();
b2b13ee2 1530}
7088dbe9 1531//____________________________________________________________________________
b2b13ee2 1532void AliPMDDigitizer::ResetDigit()
1533{
36031625 1534 // Clears Digits
b2b13ee2 1535 fNdigit = 0;
ebd83c56 1536 if (fDigits) fDigits->Delete();
b2b13ee2 1537}
7088dbe9 1538//____________________________________________________________________________
b2b13ee2 1539
1540void AliPMDDigitizer::ResetCellADC()
1541{
7088dbe9 1542 // Clears individual cells edep and track number
36031625 1543 for (Int_t i = 0; i < fgkTotUM; i++)
b2b13ee2 1544 {
36031625 1545 for (Int_t j = 0; j < fgkRow; j++)
b2b13ee2 1546 {
36031625 1547 for (Int_t k = 0; k < fgkCol; k++)
b2b13ee2 1548 {
920e13db 1549 fCPV[i][j][k] = 0.;
1550 fPRE[i][j][k] = 0.;
1551 fCPVTrackNo[i][j][k] = 0;
1552 fPRETrackNo[i][j][k] = 0;
1553 fCPVTrackPid[i][j][k] = -1;
1554 fPRETrackPid[i][j][k] = -1;
b2b13ee2 1555 }
1556 }
1557 }
1558}
7088dbe9 1559//____________________________________________________________________________
b2b13ee2 1560
1561void AliPMDDigitizer::UnLoad(Option_t *option)
1562{
36031625 1563 // Unloads all the root files
1564 //
b2b13ee2 1565 const char *cS = strstr(option,"S");
1566 const char *cD = strstr(option,"D");
1567
1568 fRunLoader->UnloadgAlice();
1569 fRunLoader->UnloadHeader();
1570 fRunLoader->UnloadKinematics();
1571
1572 if (cS)
1573 {
36031625 1574 fPMDLoader->UnloadHits();
b2b13ee2 1575 }
1576 if (cD)
1577 {
36031625 1578 fPMDLoader->UnloadHits();
1579 fPMDLoader->UnloadSDigits();
b2b13ee2 1580 }
1581}
09a06455 1582
1583//----------------------------------------------------------------------
1584Float_t AliPMDDigitizer::Gain(Int_t det, Int_t smn, Int_t row, Int_t col) const
1585{
1586 // returns of the gain of the cell
1587 // Added this method by ZA
1588
1589 //cout<<" I am here in gain "<<fCalibData<< "smn,row, col "<<smn
1590 //<<" "<<row<<" "<<col<<endl;
1591
35535af7 1592 if(!fCalibGain) {
09a06455 1593 AliError("No calibration data loaded from CDB!!!");
50555ba1 1594 return 1;
09a06455 1595 }
1596
1597 Float_t GainFact;
35535af7 1598 GainFact = fCalibGain->GetGainFact(det,smn,row,col);
09a06455 1599 return GainFact;
1600}
1601//----------------------------------------------------------------------
35535af7 1602AliPMDCalibData* AliPMDDigitizer::GetCalibGain() const
09a06455 1603{
1604 // The run number will be centralized in AliCDBManager,
1605 // you don't need to set it here!
1606 // Added this method by ZA
0dd3d6f9 1607 // Cleaned up by Alberto
35535af7 1608 AliCDBEntry *entry = AliCDBManager::Instance()->Get("PMD/Calib/Gain");
09a06455 1609
0dd3d6f9 1610 if(!entry) AliFatal("Calibration object retrieval failed!");
09a06455 1611
1612 AliPMDCalibData *calibdata=0;
1613 if (entry) calibdata = (AliPMDCalibData*) entry->GetObject();
1614
0dd3d6f9 1615 if (!calibdata) AliFatal("No calibration data from calibration database !");
09a06455 1616
1617 return calibdata;
1618}
35535af7 1619//----------------------------------------------------------------------
1620AliPMDPedestal* AliPMDDigitizer::GetCalibPed() const
1621{
1622 // The run number will be centralized in AliCDBManager,
1623 // you don't need to set it here!
1624
1625 AliCDBEntry *entry = AliCDBManager::Instance()->Get("PMD/Calib/Ped");
1626
1627 if(!entry) AliFatal("Pedestal object retrieval failed!");
1628
1629 AliPMDPedestal *pedestal=0;
1630 if (entry) pedestal = (AliPMDPedestal*) entry->GetObject();
1631
1632 if (!pedestal) AliFatal("No pedestal data from calibration database !");
1633
1634 return pedestal;
1635}