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