]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PMD/AliPMDDigitizer.cxx
update from Marco
[u/mrichter/AliRoot.git] / PMD / AliPMDDigitizer.cxx
index f8d0c5ff16f0afc7ec0035b11dc5ccc2c0732974..d77c940d4097a40d4e2ba42434259466fcbee49c 100644 (file)
@@ -40,7 +40,7 @@
 #include "AliLoader.h"
 #include "AliConfig.h"
 #include "AliMagF.h"
-#include "AliRunDigitizer.h"
+#include "AliDigitizationInput.h"
 #include "AliDigitizer.h"
 #include "AliHeader.h"
 #include "AliCDBManager.h"
@@ -128,8 +128,8 @@ AliPMDDigitizer & AliPMDDigitizer::operator=(const AliPMDDigitizer& /*digitizer*
   return *this;
 }
 //____________________________________________________________________________
-AliPMDDigitizer::AliPMDDigitizer(AliRunDigitizer* manager):
-  AliDigitizer(manager),
+AliPMDDigitizer::AliPMDDigitizer(AliDigitizationInput* digInput):
+  AliDigitizer(digInput),
   fRunLoader(0),
   fPMDHit(0),
   fPMD(0),
@@ -261,16 +261,16 @@ void AliPMDDigitizer::Hits2SDigits(Int_t ievt)
 
   const Int_t kPi0 = 111;
   const Int_t kGamma = 22;
-  Int_t npmd;
-  Int_t trackno;
-  Int_t smnumber;
-  Int_t trackpid;
-  Int_t mtrackno;
-  Int_t mtrackpid;
-
-  Float_t xPos, yPos, zPos;
+  Int_t npmd = 0;
+  Int_t trackno = 0;
+  Int_t smnumber = 0;
+  Int_t trackpid = 0;
+  Int_t mtrackno = 0;
+  Int_t mtrackpid = 0;
+
+  Float_t xPos = 0., yPos = 0., zPos = 0.;
   Int_t xpad = -1, ypad = -1;
-  Float_t edep;
+  Float_t edep = 0.;
   Float_t vx = -999.0, vy = -999.0, vz = -999.0;
 
 
@@ -280,15 +280,7 @@ void AliPMDDigitizer::Hits2SDigits(Int_t ievt)
   AliDebug(1,Form("Event Number = %d",ievt));
   Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
   AliDebug(1,Form("Number of Particles = %d",nparticles));
-  //
-  // Ajay on 24th May 2009
-  //No of primary
-  Int_t nprimary   = fRunLoader->GetHeader()->GetNprimary();
-  Int_t  * mtraPid  = new Int_t [nprimary]; 
-  for(Int_t i = 0; i < nprimary; i++)
-    {
-      mtraPid[i] = -1; 
-    }
+
   //
 
   fRunLoader->GetEvent(ievt);
@@ -340,7 +332,8 @@ void AliPMDDigitizer::Hits2SDigits(Int_t ievt)
                  statusOld   = -1;
                }
              Int_t igstatus = 0;
-             //------------------modified by Mriganka ----------------------
+
+             Int_t trnotemp = trackno;    // Modified on 25th Nov 2009
              if(ks==1||(imo = mparticle->GetFirstMother())<0 ){
                vx = mparticle->Vx();
                vy = mparticle->Vy();
@@ -361,9 +354,19 @@ void AliPMDDigitizer::Hits2SDigits(Int_t ievt)
                  vy = mparticle->Vy();
                  vz = mparticle->Vz();
                  
-                 trackno=imo;
-                               
+                 // Modified on 25th Nov 2009
+
+                 trnotemp = trackno;
+                  if(trackpid == 111)
+                   {
+                     trackno = trnotemp;
                    }
+                  if(trackpid != 111)
+                   {
+                     trackno=imo;
+                   }
+                 // end of modification on 25th Nov 2009
+               }
         
              if(trackpid==kGamma||trackpid==11||trackpid==-11||
                 trackpid==kPi0)igstatus=1;
@@ -372,13 +375,14 @@ void AliPMDDigitizer::Hits2SDigits(Int_t ievt)
              trackpid=trackpidOld;
              trackno=tracknoOld;
              
-             mtraPid[trackno] = mtrackpid;
-             
              //-----------------end of modification----------------
+             Float_t ptime = fPMDHit->GetTime()*1e6;   // time in microsec
+             if (ptime < 0. || ptime > 1.2) continue;  
+
              xPos = fPMDHit->X();
              yPos = fPMDHit->Y();
              zPos = fPMDHit->Z();
-             
+
              edep       = fPMDHit->GetEnergy();
              Int_t vol1 = fPMDHit->GetVolume(1); // Column
              Int_t vol2 = fPMDHit->GetVolume(2); // Row
@@ -415,18 +419,18 @@ void AliPMDDigitizer::Hits2SDigits(Int_t ievt)
                }
 
              AliDebug(2,Form("Zposition = %f Edeposition = %f",zPos,edep));
-             //Float_t zposition = TMath::Abs(zPos);
-             if (zPos < fZPos)
+
+             if (vol7 < 24)
                {
-                 // CPV
-                 fDetNo = 1;
+                 // PRE
+                 fDetNo = 0;
                }
-             else if (zPos > fZPos)
+             else
                {
-                 // PMD
-                 fDetNo = 0;
+                 // CPV
+                 fDetNo = 1;
                }
-             //Int_t smn = smnumber - 1;
+
              Int_t smn = smnumber;
              Int_t ixx = xpad     - 1;
              Int_t iyy = ypad     - 1;
@@ -469,18 +473,19 @@ void AliPMDDigitizer::Hits2SDigits(Int_t ievt)
                    {
                      deltaE = fPRE[ism][jrow][kcol];
                      trno   = fPRETrackNo[ism][jrow][kcol];
-                     trpid  = mtraPid[trno]; // added
                      detno = 0;
                    }
                  else if (idet == 1)
                    {
                      deltaE = fCPV[ism][jrow][kcol];
                      trno   = fCPVTrackNo[ism][jrow][kcol];
-                     trpid  = mtraPid[trno]; // added
                      detno = 1;
                    }
                  if (deltaE > 0.)
                    {
+                     // Natasha
+                     TParticle *mparticle = gAlice->GetMCApp()->Particle(trno);
+                     trpid = mparticle->GetPdgCode();
                      AddSDigit(trno,trpid,detno,ism,jrow,kcol,deltaE);
                    }
                }
@@ -491,7 +496,6 @@ void AliPMDDigitizer::Hits2SDigits(Int_t ievt)
     }
   fPMDLoader->WriteSDigits("OVERWRITE");
   ResetCellADC();
-  delete [] mtraPid;
 }
 //____________________________________________________________________________
 
@@ -502,16 +506,16 @@ void AliPMDDigitizer::Hits2Digits(Int_t ievt)
   //
   const Int_t kPi0 = 111;
   const Int_t kGamma = 22;
-  Int_t npmd;
-  Int_t trackno;
-  Int_t smnumber;
-  Int_t trackpid;
-  Int_t mtrackno;
-  Int_t mtrackpid;
-
-  Float_t xPos, yPos, zPos;
+  Int_t npmd      = 0;
+  Int_t trackno   = 0;
+  Int_t smnumber  = 0;
+  Int_t trackpid  = 0;
+  Int_t mtrackno  = 0;
+  Int_t mtrackpid = 0;
+
+  Float_t xPos = 0., yPos = 0., zPos = 0.;
   Int_t xpad = -1, ypad = -1;
-  Float_t edep;
+  Float_t edep = 0.;
   Float_t vx = -999.0, vy = -999.0, vz = -999.0;
 
   if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
@@ -521,17 +525,6 @@ void AliPMDDigitizer::Hits2Digits(Int_t ievt)
   Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
   AliDebug(1,Form("Number of Particles = %d", nparticles));
 
-  //
-  // Ajay on 24th May 2009
-  //No of primary
-  Int_t nprimary   = fRunLoader->GetHeader()->GetNprimary();
-  Int_t  * mtraPid  = new Int_t [nprimary]; 
-  for(Int_t i = 0; i < nprimary; i++)
-    {
-      mtraPid[i] = -1; 
-    }
-  //
-
   fRunLoader->GetEvent(ievt);
   // ------------------------------------------------------- //
   // Pointer to specific detector hits.
@@ -590,7 +583,8 @@ void AliPMDDigitizer::Hits2Digits(Int_t ievt)
                }
 
              Int_t igstatus = 0;
-             //-----------------------modified by Mriganka ------------------
+
+             Int_t trnotemp = trackno;   // modified on 25th Nov 2009
              if(ks==1||(imo = mparticle->GetFirstMother())<0 ){
                vx = mparticle->Vx();
                vy = mparticle->Vy();
@@ -611,9 +605,18 @@ void AliPMDDigitizer::Hits2Digits(Int_t ievt)
                  vy = mparticle->Vy();
                  vz = mparticle->Vz();
                  
-                 trackno=imo;
-                               
+                 // Modified on 25th Nov 2009
+                 
+                 trnotemp = trackno;
+                  if(trackpid == 111)
+                   {
+                     trackno = trnotemp;
+                   }
+                  if(trackpid != 111)
+                   {
+                     trackno=imo;
                    }
+               }
         
              if(trackpid==kGamma||trackpid==11||trackpid==-11||trackpid==kPi0)
                igstatus=1;
@@ -621,10 +624,10 @@ void AliPMDDigitizer::Hits2Digits(Int_t ievt)
              mtrackno=trackno;
              trackpid=trackpidOld;
              trackno=tracknoOld;
+
+             Float_t ptime = fPMDHit->GetTime()*1e6;
+             if (ptime < 0. || ptime > 1.2) continue;
              
-             mtraPid[mtrackno] = mtrackpid;  // added by Ajay
-             
-             //-----------------end of modification----------------
              xPos = fPMDHit->X();
              yPos = fPMDHit->Y();
              zPos = fPMDHit->Z();
@@ -662,21 +665,18 @@ void AliPMDDigitizer::Hits2Digits(Int_t ievt)
                  ypad = vol2;
                }
 
-             AliDebug(2,Form("ZPosition = %f Edeposition = %d",zPos,edep));
-             //Float_t zposition = TMath::Abs(zPos);
+             AliDebug(2,Form("ZPosition = %f Edeposition = %f",zPos,edep));
 
-             if (zPos < fZPos)
+             if (vol7 < 24)
                {
-                 // CPV
-                 fDetNo = 1;
+                 // PRE
+                 fDetNo = 0;
                }
-             else if (zPos > fZPos)
+             else
                {
-                 // PMD
-                 fDetNo = 0;
+                 fDetNo = 1;
                }
 
-             //Int_t smn = smnumber - 1;
              Int_t smn = smnumber;
              Int_t ixx = xpad     - 1;
              Int_t iyy = ypad     - 1;
@@ -703,12 +703,12 @@ void AliPMDDigitizer::Hits2Digits(Int_t ievt)
   TrackAssignment2Cell();
   ResetCell();
 
-  Float_t gain1;
-  Float_t adc;
+  Float_t gain1  = 1.;
+  Float_t adc    = 0. ;
   Float_t deltaE = 0.;
-  Int_t detno = 0;
-  Int_t trno = 1;
-  Int_t trpid = -99;
+  Int_t detno    = 0;
+  Int_t trno     = 1;
+  Int_t trpid    = -99;
 
   for (Int_t idet = 0; idet < 2; idet++)
   {
@@ -722,14 +722,12 @@ void AliPMDDigitizer::Hits2Digits(Int_t ievt)
                  {
                      deltaE = fPRE[ism][jrow][kcol];
                      trno   = fPRETrackNo[ism][jrow][kcol];
-                     trpid  = mtraPid[trno];   // added
                      detno  = 0;
                  }
                  else if (idet == 1)
                  {
                      deltaE = fCPV[ism][jrow][kcol];
                      trno   = fCPVTrackNo[ism][jrow][kcol];
-                     trpid  = mtraPid[trno];   // added
                      detno  = 1;
                  }
                  if (deltaE > 0.)
@@ -759,6 +757,10 @@ void AliPMDDigitizer::Hits2Digits(Int_t ievt)
                      if (adc > 0.)
                      {
                          adc += (pedmean + 3.0*pedrms);
+                         TParticle *mparticle
+                           = gAlice->GetMCApp()->Particle(trno);
+                         trpid = mparticle->GetPdgCode();
+                         
                          AddDigit(trno,trpid,detno,ism,jrow,kcol,adc);
                      }
                  }
@@ -772,7 +774,6 @@ void AliPMDDigitizer::Hits2Digits(Int_t ievt)
   fPMDLoader->WriteDigits("OVERWRITE");
   ResetCellADC();
 
-  delete [] mtraPid;
 }
 //____________________________________________________________________________
 
@@ -806,9 +807,9 @@ void AliPMDDigitizer::SDigits2Digits(Int_t ievt)
   if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
   treeD->Branch("PMDDigit", &fDigits, bufsize);
 
-  Int_t   trno, trpid, det, smn;
-  Int_t   irow, icol;
-  Float_t edep, adc;
+  Int_t   trno = 1, trpid = 0, det = 0, smn = 0;
+  Int_t   irow = 0, icol = 0;
+  Float_t edep = 0., adc = 0.;
 
   Int_t nmodules = (Int_t) treeS->GetEntries();
   AliDebug(1,Form("Number of modules = %d",nmodules));
@@ -862,7 +863,7 @@ void AliPMDDigitizer::SDigits2Digits(Int_t ievt)
 
 }
 //____________________________________________________________________________
-void AliPMDDigitizer::Exec(Option_t *option)
+void AliPMDDigitizer::Digitize(Option_t *option)
 {
   // Does the event merging and digitization
   const char *cdeb = strstr(option,"deb");
@@ -871,17 +872,17 @@ void AliPMDDigitizer::Exec(Option_t *option)
       AliDebug(100," *** PMD Exec is called ***");
     }
 
-  Int_t ninputs = fManager->GetNinputs();
+  Int_t ninputs = fDigInput->GetNinputs();
   AliDebug(1,Form("Number of files to be processed = %d",ninputs));
   ResetCellADC();
 
   for (Int_t i = 0; i < ninputs; i++)
     {
-      Int_t troffset = fManager->GetMask(i);
+      Int_t troffset = fDigInput->GetMask(i);
       MergeSDigits(i, troffset);
     }
 
-  fRunLoader = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
+  fRunLoader = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
   fPMD  = (AliPMD*)gAlice->GetDetector("PMD");
   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
   if (fPMDLoader == 0x0)
@@ -899,7 +900,7 @@ void AliPMDDigitizer::Exec(Option_t *option)
   if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
   treeD->Branch("PMDDigit", &fDigits, bufsize);
 
-  Float_t adc;
+  Float_t adc    = 0.;
   Float_t deltaE = 0.;
   Int_t detno    = 0;
   Int_t trno     = 1;
@@ -978,7 +979,7 @@ void AliPMDDigitizer::TrackAssignment2CPVCell()
   // This method added by Ajay
   Bool_t jsort = false;
 
-  Int_t i, j, k;
+  Int_t i = 0, j = 0, k = 0;
 
   Int_t   *status1;
   Int_t   *status2;
@@ -1031,8 +1032,8 @@ void AliPMDDigitizer::TrackAssignment2CPVCell()
 
   Int_t nentries = fCPVCell.GetEntries();
  
-  Int_t   mtrackno, ism, ixp, iyp;
-  Float_t edep;
+  Int_t   mtrackno = 0, ism = 0, ixp = 0, iyp = 0;
+  Float_t edep = 0.;
   for (i = 0; i < nentries; i++)
     {
       AliPMDcell* cpvcell = (AliPMDcell*)fCPVCell.UncheckedAt(i);
@@ -1048,9 +1049,9 @@ void AliPMDDigitizer::TrackAssignment2CPVCell()
       fCPVCounter[ism][ixp][iyp]++;
     }
   
-  Int_t iz, il;
-  Int_t im, ix, iy;
-  Int_t nn;
+  Int_t iz = 0, il = 0;
+  Int_t im = 0, ix = 0, iy = 0;
+  Int_t nn = 0;
   for (im=0; im<fgkTotUM; im++)
     {
       for (ix=0; ix<fgkRow; ix++)
@@ -1063,7 +1064,7 @@ void AliPMDDigitizer::TrackAssignment2CPVCell()
                  // This block handles if a cell is fired
                  // many times by many tracks
                  status1  = new Int_t[nn];
-                 status2  = new Int_t[nn];
+                 status2  = new Int_t[2*nn];
                  trnarray = new Int_t[nn];
                  for (iz = 0; iz < nn; iz++)
                    {
@@ -1174,7 +1175,7 @@ void AliPMDDigitizer::TrackAssignment2CPVCell()
 void AliPMDDigitizer::MergeSDigits(Int_t filenumber, Int_t troffset)
 {
   // merging sdigits
-  fRunLoader = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(filenumber));
+  fRunLoader = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(filenumber));
   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
   fPMDLoader->LoadSDigits("read");
   TTree* treeS = fPMDLoader->TreeS();
@@ -1183,9 +1184,9 @@ void AliPMDDigitizer::MergeSDigits(Int_t filenumber, Int_t troffset)
   if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
   branch->SetAddress(&fSDigits);
 
-  Int_t   itrackno, itrackpid, idet, ism;
-  Int_t   ixp, iyp;
-  Float_t edep;
+  Int_t   itrackno = 1, itrackpid = 0, idet = 0, ism = 0;
+  Int_t   ixp = 0, iyp = 0;
+  Float_t edep = 0.;
   Int_t nmodules = (Int_t) treeS->GetEntries();
   AliDebug(1,Form("Number of Modules in the treeS = %d",nmodules));
   AliDebug(1,Form("Track Offset = %d",troffset));
@@ -1236,7 +1237,7 @@ void AliPMDDigitizer::TrackAssignment2Cell()
   //
   Bool_t jsort = false;
 
-  Int_t i, j, k;
+  Int_t i = 0, j = 0, k = 0;
 
   Int_t   *status1;
   Int_t   *status2;
@@ -1308,9 +1309,9 @@ void AliPMDDigitizer::TrackAssignment2Cell()
       fPRECounter[ism][ixp][iyp]++;
     }
   
-  Int_t iz, il;
-  Int_t im, ix, iy;
-  Int_t nn;
+  Int_t iz = 0, il = 0;
+  Int_t im = 0, ix = 0, iy = 0;
+  Int_t nn = 0;
   
   for (im=0; im<fgkTotUM; im++)
     {
@@ -1324,7 +1325,7 @@ void AliPMDDigitizer::TrackAssignment2Cell()
                  // This block handles if a cell is fired
                  // many times by many tracks
                  status1  = new Int_t[nn];
-                 status2  = new Int_t[nn];
+                 status2  = new Int_t[2*nn];
                  trnarray = new Int_t[nn];
                  for (iz = 0; iz < nn; iz++)
                    {
@@ -1434,29 +1435,36 @@ void AliPMDDigitizer::MeV2ADC(Float_t mev, Float_t & adc) const
 {
   // This converts the simulated edep to ADC according to the
   // Test Beam Data
-  //PS Test in September 2003 and 2006
+  // PS Test in June 2010, Voltage @ 1300 V
   // KeV - ADC conversion for 12bit ADC
-  // Modified by Ajay
-  const Float_t kConstant   =  0.07;
-  const Float_t kErConstant =  0.1;
-  const Float_t kSlope      = 76.0;
-  const Float_t kErSlope    =  5.0;
-  
-  Float_t cons   = gRandom->Gaus(kConstant,kErConstant);
-  Float_t slop   = gRandom->Gaus(kSlope,kErSlope);
+  // MPV data used for the fit and taken here
+
+  // constants are from Test Beam 2010
   
-  Float_t adc12bit = slop*mev*0.001 + cons;
+  const Float_t kConstant   = 0.612796;
+  const Float_t kSlope      = 130.158;
   
-                          
-  if(adc12bit < 1600.0)
+  Float_t adc12bit = kSlope*mev*0.001 + kConstant;
+  if (adc12bit < 0.) adc12bit = 0.;
+
+  //Introducing Readout Resolution for ALICE-PMD
+
+  Float_t sigrr     = 0.605016 - 0.000273*adc12bit + 6.54e-8*adc12bit*adc12bit;
+  Float_t adcwithrr = gRandom->Gaus(adc12bit,sigrr);
+
+  if(adcwithrr < 0.)
     {
-      adc = (Float_t) adc12bit;
+      adc = 0.;
     }
-  else if (adc12bit >= 1600.0)
+  else if(adcwithrr >= 0. && adcwithrr < 1600.0)
+    {
+      adc = adcwithrr;
+    }
+  else if (adcwithrr >= 1600.0)
     {
       adc = 1600.0;
     }
+
 }
 //____________________________________________________________________________
 void AliPMDDigitizer::AddSDigit(Int_t trnumber, Int_t trpid, Int_t det,