1 //-----------------------------------------------------//
3 // Source File : PMDDigitizer.cxx, Version 00 //
5 // Date : September 20 2002 //
7 //-----------------------------------------------------//
12 #include <TGeometry.h>
13 #include <TObjArray.h>
14 #include <TClonesArray.h>
17 #include <TParticle.h>
22 #include "AliDetector.h"
23 #include "AliRunLoader.h"
24 #include "AliLoader.h"
25 #include "AliConfig.h"
27 #include "AliRunDigitizer.h"
28 #include "AliHeader.h"
30 #include "AliPMDcell.h"
31 #include "AliPMDsdigit.h"
32 #include "AliPMDdigit.h"
33 #include "AliPMDDigitizer.h"
34 #include "AliPMDClustering.h"
37 ClassImp(AliPMDDigitizer)
41 AliPMDDigitizer::AliPMDDigitizer()
43 if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
45 if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
48 for (Int_t i = 0; i < fTotUM; i++)
50 for (Int_t j = 0; j < fRow; j++)
52 for (Int_t k = 0; k < fCol; k++)
60 if (!fCell) fCell = new TObjArray();
62 fZPos = 361.5; // in units of cm, This is the default position of PMD
64 AliPMDDigitizer::~AliPMDDigitizer()
73 void AliPMDDigitizer::OpengAliceFile(Char_t *file, Option_t *option)
76 fRunLoader = AliRunLoader::Open(file,AliConfig::fgkDefaultEventFolderName,
81 Error("Open","Can not open session for file %s.",file);
84 fRunLoader->LoadgAlice();
85 fRunLoader->LoadHeader();
86 fRunLoader->LoadKinematics();
88 gAlice = fRunLoader->GetAliRun();
92 printf("<AliPMDdigitizer::Open> ");
93 printf("AliRun object found on file.\n");
97 printf("<AliPMDdigitizer::Open> ");
98 printf("Could not find AliRun object.\n");
101 PMD = (AliPMD*)gAlice->GetDetector("PMD");
102 pmdloader = fRunLoader->GetLoader("PMDLoader");
103 if (pmdloader == 0x0)
105 cerr<<"Hits2Digits : Can not find PMD or PMDLoader\n";
108 const char *cHS = strstr(option,"HS");
109 const char *cHD = strstr(option,"HD");
110 const char *cSD = strstr(option,"SD");
114 pmdloader->LoadHits("READ");
115 pmdloader->LoadSDigits("recreate");
119 pmdloader->LoadHits("READ");
120 pmdloader->LoadDigits("recreate");
124 pmdloader->LoadSDigits("READ");
125 pmdloader->LoadDigits("recreate");
129 void AliPMDDigitizer::Hits2SDigits(Int_t ievt)
131 cout << " -------- Beginning of Hits2SDigits ----------- " << endl;
142 Float_t xPos, yPos, zPos;
143 Int_t xpad = -1, ypad = -1;
145 Float_t vx = -999.0, vy = -999.0, vz = -999.0;
150 printf("Event Number = %d \n",ievt);
151 Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
152 printf("Number of Particles = %d \n", nparticles);
153 fRunLoader->GetEvent(ievt);
154 Particles = gAlice->GetMCApp()->Particles();
155 // ------------------------------------------------------- //
156 // Pointer to specific detector hits.
157 // Get pointers to Alice detectors and Hits containers
159 treeH = pmdloader->TreeH();
161 Int_t ntracks = (Int_t) treeH->GetEntries();
162 printf("Number of Tracks in the TreeH = %d \n", ntracks);
164 treeS = pmdloader->TreeS();
167 pmdloader->MakeTree("S");
168 treeS = pmdloader->TreeS();
170 Int_t bufsize = 16000;
171 treeS->Branch("PMDSDigit", &fSDigits, bufsize);
173 if (PMD) PMDhits = PMD->Hits();
175 // Start loop on tracks in the hits containers
178 for (Int_t track=0; track<ntracks;track++)
181 treeH->GetEvent(track);
185 npmd = PMDhits->GetEntriesFast();
186 for (int ipmd = 0; ipmd < npmd; ipmd++)
188 pmdHit = (AliPMDhit*) PMDhits->UncheckedAt(ipmd);
189 trackno = pmdHit->GetTrack();
191 // get kinematics of the particles
193 particle = gAlice->GetMCApp()->Particle(trackno);
194 trackpid = particle->GetPdgCode();
203 TParticle* mparticle = particle;
204 Int_t trackno_old=0, trackpid_old=0, status_old = 0;
205 if (mparticle->GetFirstMother() == -1)
207 trackno_old = trackno;
208 trackpid_old = trackpid;
213 while((imo = mparticle->GetFirstMother()) >= 0)
216 mparticle = gAlice->GetMCApp()->Particle(imo);
217 id_mo = mparticle->GetPdgCode();
219 vx = mparticle->Vx();
220 vy = mparticle->Vy();
221 vz = mparticle->Vz();
223 //printf("==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, id_mo, vx, vy, vz);
224 //fprintf(fpw1,"==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, id_mo, vx, vy, vz);
225 if ((id_mo == kGamma || id_mo == -11 || id_mo == 11) && vx == 0. && vy == 0. && vz == 0.)
233 if (id_mo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
242 if (id_mo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
252 if (status_old == -1)
254 mtrackno = trackno_old;
255 mtrackpid = trackpid_old;
261 edep = pmdHit->fEnergy;
262 Int_t Vol1 = pmdHit->fVolume[1]; // Column
263 Int_t Vol2 = pmdHit->fVolume[2]; // Row
264 Int_t Vol3 = pmdHit->fVolume[3]; // UnitModule
265 Int_t Vol6 = pmdHit->fVolume[6]; // SuperModule
266 // -----------------------------------------//
267 // For Super Module 1 & 2 //
268 // nrow = 96, ncol = 48 //
269 // For Super Module 3 & 4 //
270 // nrow = 48, ncol = 96 //
271 // -----------------------------------------//
273 smnumber = (Vol6-1)*6 + Vol3;
275 if (Vol6 == 1 || Vol6 == 2)
280 else if (Vol6 == 3 || Vol6 == 4)
286 //cout << "zpos = " << zPos << " edep = " << edep << endl;
288 Float_t zposition = TMath::Abs(zPos);
289 if (zposition < fZPos)
294 else if (zposition > fZPos)
299 Int_t smn = smnumber - 1;
300 Int_t ixx = xpad - 1;
301 Int_t iyy = ypad - 1;
304 fPMD[smn][ixx][iyy] += edep;
305 fPMDCounter[smn][ixx][iyy]++;
307 pmdcell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
313 fCPV[smn][ixx][iyy] += edep;
314 fCPVTrackNo[smn][ixx][iyy] = mtrackno;
318 } // Track Loop ended
320 TrackAssignment2Cell();
328 for (Int_t idet = 0; idet < 2; idet++)
330 for (Int_t ism = 0; ism < fTotUM; ism++)
332 for (Int_t jrow = 0; jrow < fRow; jrow++)
334 for (Int_t kcol = 0; kcol < fCol; kcol++)
336 cellno = jrow*fCol + kcol;
339 deltaE = fPMD[ism][jrow][kcol];
340 trno = fPMDTrackNo[ism][jrow][kcol];
345 deltaE = fCPV[ism][jrow][kcol];
346 trno = fCPVTrackNo[ism][jrow][kcol];
351 AddSDigit(trno,detno,ism,cellno,deltaE);
359 pmdloader->WriteSDigits("OVERWRITE");
362 // cout << " -------- End of Hits2SDigit ----------- " << endl;
365 void AliPMDDigitizer::Hits2Digits(Int_t ievt)
376 Float_t xPos, yPos, zPos;
377 Int_t xpad = -1, ypad = -1;
379 Float_t vx = -999.0, vy = -999.0, vz = -999.0;
384 printf("Event Number = %d \n",ievt);
386 Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
387 printf("Number of Particles = %d \n", nparticles);
388 fRunLoader->GetEvent(ievt);
389 Particles = gAlice->GetMCApp()->Particles();
390 // ------------------------------------------------------- //
391 // Pointer to specific detector hits.
392 // Get pointers to Alice detectors and Hits containers
394 PMD = (AliPMD*)gAlice->GetDetector("PMD");
395 pmdloader = fRunLoader->GetLoader("PMDLoader");
397 if (pmdloader == 0x0)
399 cerr<<"Hits2Digits method : Can not find PMD or PMDLoader\n";
401 treeH = pmdloader->TreeH();
402 Int_t ntracks = (Int_t) treeH->GetEntries();
403 printf("Number of Tracks in the TreeH = %d \n", ntracks);
404 pmdloader->LoadDigits("recreate");
405 treeD = pmdloader->TreeD();
408 pmdloader->MakeTree("D");
409 treeD = pmdloader->TreeD();
411 Int_t bufsize = 16000;
412 treeD->Branch("PMDDigit", &fDigits, bufsize);
414 if (PMD) PMDhits = PMD->Hits();
416 // Start loop on tracks in the hits containers
418 for (Int_t track=0; track<ntracks;track++)
421 treeH->GetEvent(track);
425 npmd = PMDhits->GetEntriesFast();
426 for (int ipmd = 0; ipmd < npmd; ipmd++)
428 pmdHit = (AliPMDhit*) PMDhits->UncheckedAt(ipmd);
429 trackno = pmdHit->GetTrack();
431 // get kinematics of the particles
433 particle = gAlice->GetMCApp()->Particle(trackno);
434 trackpid = particle->GetPdgCode();
443 TParticle* mparticle = particle;
444 Int_t trackno_old=0, trackpid_old=0, status_old = 0;
445 if (mparticle->GetFirstMother() == -1)
447 trackno_old = trackno;
448 trackpid_old = trackpid;
453 while((imo = mparticle->GetFirstMother()) >= 0)
456 mparticle = gAlice->GetMCApp()->Particle(imo);
457 id_mo = mparticle->GetPdgCode();
459 vx = mparticle->Vx();
460 vy = mparticle->Vy();
461 vz = mparticle->Vz();
463 //printf("==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, id_mo, vx, vy, vz);
464 //fprintf(fpw1,"==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, id_mo, vx, vy, vz);
465 if ((id_mo == kGamma || id_mo == -11 || id_mo == 11) && vx == 0. && vy == 0. && vz == 0.)
473 if (id_mo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
482 if (id_mo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
492 if (status_old == -1)
494 mtrackno = trackno_old;
495 mtrackpid = trackpid_old;
501 edep = pmdHit->fEnergy;
502 Int_t Vol1 = pmdHit->fVolume[1]; // Column
503 Int_t Vol2 = pmdHit->fVolume[2]; // Row
504 Int_t Vol3 = pmdHit->fVolume[3]; // UnitModule
505 Int_t Vol6 = pmdHit->fVolume[6]; // SuperModule
507 // -----------------------------------------//
508 // For Super Module 1 & 2 //
509 // nrow = 96, ncol = 48 //
510 // For Super Module 3 & 4 //
511 // nrow = 48, ncol = 96 //
512 // -----------------------------------------//
514 smnumber = (Vol6-1)*6 + Vol3;
516 if (Vol6 == 1 || Vol6 == 2)
521 else if (Vol6 == 3 || Vol6 == 4)
527 //cout << "-zpos = " << -zPos << endl;
529 Float_t zposition = TMath::Abs(zPos);
531 if (zposition < fZPos)
536 else if (zposition > fZPos)
542 Int_t smn = smnumber - 1;
543 Int_t ixx = xpad - 1;
544 Int_t iyy = ypad - 1;
547 fPMD[smn][ixx][iyy] += edep;
548 fPMDCounter[smn][ixx][iyy]++;
550 pmdcell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
556 fCPV[smn][ixx][iyy] += edep;
557 fCPVTrackNo[smn][ixx][iyy] = mtrackno;
561 } // Track Loop ended
563 TrackAssignment2Cell();
571 for (Int_t idet = 0; idet < 2; idet++)
573 for (Int_t ism = 0; ism < fTotUM; ism++)
575 for (Int_t jrow = 0; jrow < fRow; jrow++)
577 for (Int_t kcol = 0; kcol < fCol; kcol++)
579 cellno = jrow*fCol + kcol;
582 deltaE = fPMD[ism][jrow][kcol];
583 trno = fPMDTrackNo[ism][jrow][kcol];
588 deltaE = fCPV[ism][jrow][kcol];
589 trno = fCPVTrackNo[ism][jrow][kcol];
594 AddDigit(trno,detno,ism,cellno,deltaE);
598 } // supermodule loop
603 pmdloader->WriteDigits("OVERWRITE");
606 // cout << " -------- End of Hits2Digit ----------- " << endl;
610 void AliPMDDigitizer::SDigits2Digits(Int_t ievt)
612 // cout << " -------- Beginning of SDigits2Digit ----------- " << endl;
613 fRunLoader->GetEvent(ievt);
615 treeS = pmdloader->TreeS();
616 AliPMDsdigit *pmdsdigit;
617 TBranch *branch = treeS->GetBranch("PMDSDigit");
618 branch->SetAddress(&fSDigits);
620 treeD = pmdloader->TreeD();
623 pmdloader->MakeTree("D");
624 treeD = pmdloader->TreeD();
626 Int_t bufsize = 16000;
627 treeD->Branch("PMDDigit", &fDigits, bufsize);
629 Int_t trno, det, smn;
633 Int_t nmodules = (Int_t) treeS->GetEntries();
635 for (Int_t imodule = 0; imodule < nmodules; imodule++)
637 treeS->GetEntry(imodule);
638 Int_t nentries = fSDigits->GetLast();
639 //cout << " nentries = " << nentries << endl;
640 for (Int_t ient = 0; ient < nentries+1; ient++)
642 pmdsdigit = (AliPMDsdigit*)fSDigits->UncheckedAt(ient);
643 trno = pmdsdigit->GetTrackNumber();
644 det = pmdsdigit->GetDetector();
645 smn = pmdsdigit->GetSMNumber();
646 cellno = pmdsdigit->GetCellNumber();
647 edep = pmdsdigit->GetCellEdep();
650 AddDigit(trno,det,smn,cellno,adc);
655 pmdloader->WriteDigits("OVERWRITE");
656 // cout << " -------- End of SDigits2Digit ----------- " << endl;
659 void AliPMDDigitizer::TrackAssignment2Cell()
662 // This block assigns the cell id when there are
663 // multiple tracks in a cell according to the
666 Bool_t jsort = false;
678 PMDTrack = new Int_t ***[fTotUM];
679 PMDEdep = new Float_t ***[fTotUM];
680 for (i=0; i<fTotUM; i++)
682 PMDTrack[i] = new Int_t **[fRow];
683 PMDEdep[i] = new Float_t **[fRow];
686 for (i = 0; i < fTotUM; i++)
688 for (j = 0; j < fRow; j++)
690 PMDTrack[i][j] = new Int_t *[fCol];
691 PMDEdep[i][j] = new Float_t *[fCol];
695 for (i = 0; i < fTotUM; i++)
697 for (j = 0; j < fRow; j++)
699 for (k = 0; k < fCol; k++)
701 Int_t nn = fPMDCounter[i][j][k];
704 PMDTrack[i][j][k] = new Int_t[nn];
705 PMDEdep[i][j][k] = new Float_t[nn];
710 PMDTrack[i][j][k] = new Int_t[nn];
711 PMDEdep[i][j][k] = new Float_t[nn];
713 fPMDCounter[i][j][k] = 0;
719 Int_t nentries = fCell->GetEntries();
721 Int_t mtrackno, ism, ixp, iyp;
724 for (i = 0; i < nentries; i++)
726 pmdcell = (AliPMDcell*)fCell->UncheckedAt(i);
728 mtrackno = pmdcell->GetTrackNumber();
729 ism = pmdcell->GetSMNumber();
730 ixp = pmdcell->GetX();
731 iyp = pmdcell->GetY();
732 edep = pmdcell->GetEdep();
733 Int_t nn = fPMDCounter[ism][ixp][iyp];
734 // cout << " nn = " << nn << endl;
735 PMDTrack[ism][ixp][iyp][nn] = (Int_t) mtrackno;
736 PMDEdep[ism][ixp][iyp][nn] = edep;
737 fPMDCounter[ism][ixp][iyp]++;
744 for (im=0; im<fTotUM; im++)
746 for (ix=0; ix<fRow; ix++)
748 for (iy=0; iy<fCol; iy++)
750 nn = fPMDCounter[im][ix][iy];
753 // This block handles if a cell is fired
754 // many times by many tracks
755 status1 = new Int_t[nn];
756 status2 = new Int_t[nn];
757 trnarray = new Int_t[nn];
758 for (iz = 0; iz < nn; iz++)
760 status1[iz] = PMDTrack[im][ix][iy][iz];
762 TMath::Sort(nn,status1,status2,jsort);
763 Int_t track_old = -99999;
764 Int_t track, tr_count = 0;
765 for (iz = 0; iz < nn; iz++)
767 track = status1[status2[iz]];
768 if (track_old != track)
770 trnarray[tr_count] = track;
777 Float_t tot_edp = 0.;
778 tr_edp = new Float_t[tr_count];
779 frac_edp = new Float_t[tr_count];
780 for (il = 0; il < tr_count; il++)
783 track = trnarray[il];
784 for (iz = 0; iz < nn; iz++)
786 if (track == PMDTrack[im][ix][iy][iz])
788 tr_edp[il] += PMDEdep[im][ix][iy][iz];
791 tot_edp += tr_edp[il];
794 Float_t frac_old = 0.;
796 for (il = 0; il < tr_count; il++)
798 frac_edp[il] = tr_edp[il]/tot_edp;
799 if (frac_old < frac_edp[il])
801 frac_old = frac_edp[il];
805 fPMDTrackNo[im][ix][iy] = trnarray[il_old];
812 // This only handles if a cell is fired
815 fPMDTrackNo[im][ix][iy] = PMDTrack[im][ix][iy][0];
820 // This is if no cell is fired
821 fPMDTrackNo[im][ix][iy] = -999;
827 // Delete all the pointers
829 for (i = 0; i < fTotUM; i++)
831 for (j = 0; j < fRow; j++)
833 for (k = 0; k < fCol; k++)
835 delete [] PMDTrack[i][j][k];
836 delete [] PMDEdep[i][j][k];
841 for (i = 0; i < fTotUM; i++)
843 for (j = 0; j < fRow; j++)
845 delete [] PMDTrack[i][j];
846 delete [] PMDEdep[i][j];
850 for (i = 0; i < fTotUM; i++)
852 delete [] PMDTrack[i];
853 delete [] PMDEdep[i];
858 // End of the cell id assignment
863 void AliPMDDigitizer::MeV2ADC(Float_t mev, Float_t & adc)
869 void AliPMDDigitizer::AddSDigit(Int_t trnumber, Int_t det, Int_t smnumber,
870 Int_t cellnumber, Float_t adc)
872 TClonesArray &lsdigits = *fSDigits;
873 AliPMDsdigit *newcell;
874 newcell = new AliPMDsdigit(trnumber,det,smnumber,cellnumber,adc);
875 new(lsdigits[fNsdigit++]) AliPMDsdigit(newcell);
879 void AliPMDDigitizer::AddDigit(Int_t trnumber, Int_t det, Int_t smnumber,
880 Int_t cellnumber, Float_t adc)
882 TClonesArray &ldigits = *fDigits;
883 AliPMDdigit *newcell;
884 newcell = new AliPMDdigit(trnumber,det,smnumber,cellnumber,adc);
885 new(ldigits[fNdigit++]) AliPMDdigit(newcell);
889 Int_t AliPMDDigitizer::Convert2RealSMNumber(Int_t smnumber1)
891 Int_t smnumber = -999;
893 if (smnumber1==1) smnumber = 1;
894 if (smnumber1==2) smnumber = 10;
895 if (smnumber1==3) smnumber = 19;
896 if (smnumber1==4) smnumber = 1;
897 if (smnumber1==5) smnumber = 10;
898 if (smnumber1==6) smnumber = 19;
899 if (smnumber1==7) smnumber = 2;
900 if (smnumber1==8) smnumber = 3;
901 if (smnumber1==9) smnumber = 4;
902 if (smnumber1==10) smnumber = 5;
903 if (smnumber1==11) smnumber = 6;
904 if (smnumber1==12) smnumber = 7;
905 if (smnumber1==13) smnumber = 8;
906 if (smnumber1==14) smnumber = 9;
907 if (smnumber1==15) smnumber = 11;
908 if (smnumber1==16) smnumber = 12;
909 if (smnumber1==17) smnumber = 13;
910 if (smnumber1==18) smnumber = 14;
911 if (smnumber1==19) smnumber = 15;
912 if (smnumber1==20) smnumber = 16;
913 if (smnumber1==21) smnumber = 17;
914 if (smnumber1==22) smnumber = 18;
915 if (smnumber1==23) smnumber = 20;
916 if (smnumber1==24) smnumber = 21;
917 if (smnumber1==25) smnumber = 22;
918 if (smnumber1==26) smnumber = 23;
919 if (smnumber1==27) smnumber = 24;
920 if (smnumber1==28) smnumber = 25;
921 if (smnumber1==29) smnumber = 26;
922 if (smnumber1==30) smnumber = 27;
926 void AliPMDDigitizer::SetZPosition(Float_t zpos)
930 Float_t AliPMDDigitizer::GetZPosition() const
935 void AliPMDDigitizer::ResetCell()
938 for (Int_t i = 0; i < fTotUM; i++)
940 for (Int_t j = 0; j < fRow; j++)
942 for (Int_t k = 0; k < fCol; k++)
944 fPMDCounter[i][j][k] = 0;
949 void AliPMDDigitizer::ResetSDigit()
952 if (fSDigits) fSDigits->Clear();
954 void AliPMDDigitizer::ResetDigit()
957 if (fDigits) fDigits->Clear();
960 void AliPMDDigitizer::ResetCellADC()
962 for (Int_t i = 0; i < fTotUM; i++)
964 for (Int_t j = 0; j < fRow; j++)
966 for (Int_t k = 0; k < fCol; k++)
975 void AliPMDDigitizer::UnLoad(Option_t *option)
977 const char *cS = strstr(option,"S");
978 const char *cD = strstr(option,"D");
980 fRunLoader->UnloadgAlice();
981 fRunLoader->UnloadHeader();
982 fRunLoader->UnloadKinematics();
986 pmdloader->UnloadHits();
990 pmdloader->UnloadHits();
991 pmdloader->UnloadSDigits();