parameters have been redistributed; Hits2SDigits etc ... introduce
authorschutz <schutz@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 5 Feb 2001 09:03:17 +0000 (09:03 +0000)
committerschutz <schutz@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 5 Feb 2001 09:03:17 +0000 (09:03 +0000)
12 files changed:
PHOS/AliPHOS.cxx
PHOS/AliPHOS.h
PHOS/AliPHOSAnalyze.cxx
PHOS/AliPHOSDigit.cxx
PHOS/AliPHOSDigit.h
PHOS/AliPHOSHit.cxx
PHOS/AliPHOSHit.h
PHOS/AliPHOSv1.cxx
PHOS/AliPHOSv1.h
PHOS/AliPHOSv2.cxx
PHOS/AliPHOSv2.h
PHOS/AliPHOSv3.cxx

index 0c180e494a28918f726ae1ac8a697530152e401b..af959462ba55d2477ed1717858b53ef2119948c2 100644 (file)
@@ -32,6 +32,7 @@ class TFile;
 
 // --- Standard library ---
 
+#include <strstream.h>
 // --- AliRoot header files ---
 
 #include "AliPHOS.h"
@@ -348,11 +349,22 @@ void AliPHOS::CreateMaterials()
 //____________________________________________________________________________
 void AliPHOS::SetTreeAddress()
 { 
+
   // TBranch *branch;
   AliDetector::SetTreeAddress();
 
   TBranch * branch ;
 
+  if(fSDigits)
+    fSDigits->Clear();
+  else
+    fSDigits = new TClonesArray("AliPHOSDigit",1) ;
+
+  if (gAlice->TreeS()  && fSDigits ) {
+    branch = gAlice->TreeS()->GetBranch("PHOS");
+    if (branch) branch->SetAddress(&fSDigits) ;
+  } 
+
   if(fDigits)
     fDigits->Clear();
   else
index 60782a82540d745283a113210305dee8d7dcf8f2..1761e48f9854c0c9766c427fd1dc69c375692b61 100644 (file)
@@ -70,7 +70,7 @@ class AliPHOS : public AliDetector {
   }
  
 protected:
-  
+  TClonesArray                           *fSDigits      ; // List of summable digits
   AliPHOSRecPoint::RecPointsList         *fEmcRecPoints ; // The RecPoints (clusters) list in EMC 
   AliPHOSRecPoint::RecPointsList         *fPpsdRecPoints ;// The RecPoints (clusters) list in PPSD (veto)
   AliPHOSTrackSegment::TrackSegmentsList *fTrackSegments ;// The TrackSegment list in PHOS
index ab1c91b15fc407d0c4289ee362b97dc20d534e79..413a5c4b04c272f739e6954fa4812a3919a67e4b 100644 (file)
@@ -325,10 +325,9 @@ void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod){
        (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
       cout <<  "======= Analyze ======> Event " << ievent+1 << endl ;
     
-    //=========== Connects the various Tree's for evt
-    Int_t tracks = gAlice->GetEvent(ievent);
+    fPHOS->Enable() ;
     
-    fPHOS->Hit2Digit(tracks) ;
+    gAlice->Hits2Digits() ;
     
     //=========== Do the reconstruction
     fPHOS->Reconstruction(fRec);    
@@ -428,183 +427,183 @@ void AliPHOSAnalyze::ReadAndPrintCPV(Int_t EvFirst, Int_t EvLast)
 //____________________________________________________________________________
 void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
 {
-  //
-  // Analyzes CPV characteristics
-  // Author: Yuri Kharlov
-  // 9 October 2000
-  //
-
-  // Book histograms
-
-  TH1F *hDx   = new TH1F("hDx"  ,"CPV x-resolution@reconstruction",100,-5. , 5.);
-  TH1F *hDz   = new TH1F("hDz"  ,"CPV z-resolution@reconstruction",100,-5. , 5.);
-  TH1F *hDr   = new TH1F("hDr"  ,"CPV r-resolution@reconstruction",100, 0. , 5.);
-  TH1S *hNrp  = new TH1S("hNrp" ,"CPV rec.point multiplicity",      21,-0.5,20.5);
-  TH1S *hNrpX = new TH1S("hNrpX","CPV rec.point Phi-length"  ,      21,-0.5,20.5);
-  TH1S *hNrpZ = new TH1S("hNrpZ","CPV rec.point Z-length"    ,      21,-0.5,20.5);
-
-  cout << "Start CPV Analysis"<< endl ;
-  for ( Int_t ievent=0; ievent<Nevents; ievent++) {  
+//   //
+//   // Analyzes CPV characteristics
+//   // Author: Yuri Kharlov
+//   // 9 October 2000
+//   //
+
+//   // Book histograms
+
+//   TH1F *hDx   = new TH1F("hDx"  ,"CPV x-resolution@reconstruction",100,-5. , 5.);
+//   TH1F *hDz   = new TH1F("hDz"  ,"CPV z-resolution@reconstruction",100,-5. , 5.);
+//   TH1F *hDr   = new TH1F("hDr"  ,"CPV r-resolution@reconstruction",100, 0. , 5.);
+//   TH1S *hNrp  = new TH1S("hNrp" ,"CPV rec.point multiplicity",      21,-0.5,20.5);
+//   TH1S *hNrpX = new TH1S("hNrpX","CPV rec.point Phi-length"  ,      21,-0.5,20.5);
+//   TH1S *hNrpZ = new TH1S("hNrpZ","CPV rec.point Z-length"    ,      21,-0.5,20.5);
+
+//   cout << "Start CPV Analysis"<< endl ;
+//   for ( Int_t ievent=0; ievent<Nevents; ievent++) {  
       
-    //========== Event Number>         
-//      if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
-      cout << endl <<  "==== AnalyzeCPV ====> Event is " << ievent+1 << endl ;
+//     //========== Event Number>         
+// //      if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
+//       cout << endl <<  "==== AnalyzeCPV ====> Event is " << ievent+1 << endl ;
     
-    //=========== Connects the various Tree's for evt
-    Int_t ntracks = gAlice->GetEvent(ievent);
+//     //=========== Connects the various Tree's for evt
+//     Int_t ntracks = gAlice->GetEvent(ievent);
     
-    //========== Creating branches ===================================
-    AliPHOSRecPoint::RecPointsList ** emcRecPoints = fPHOS->EmcRecPoints() ;
-    gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , emcRecPoints  ) ;
+//     //========== Creating branches ===================================
+//     AliPHOSRecPoint::RecPointsList ** emcRecPoints = fPHOS->EmcRecPoints() ;
+//     gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , emcRecPoints  ) ;
     
-    AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->PpsdRecPoints() ;
-    gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", cpvRecPoints ) ;
+//     AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->PpsdRecPoints() ;
+//     gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", cpvRecPoints ) ;
 
-    // Create and fill arrays of hits for each CPV module
+//     // Create and fill arrays of hits for each CPV module
       
-    Int_t nOfModules = fGeom->GetNModules();
-    TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
-    Int_t iModule = 0;         
-    for (iModule=0; iModule < nOfModules; iModule++)
-      hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100);
-
-    AliPHOSCPVModule cpvModule;
-    TClonesArray    *cpvHits;
-    Int_t           nCPVhits;
-    AliPHOSCPVHit   *cpvHit;
-    TLorentzVector   p;
-    Float_t          xzgen[2];
-    Int_t            ipart;
-
-    // First go through all primary tracks and fill the arrays
-    // of hits per each CPV module
-
-    for (Int_t itrack=0; itrack<ntracks; itrack++) {
-      // Get the Hits Tree for the Primary track itrack
-      gAlice->ResetHits();
-      gAlice->TreeH()->GetEvent(itrack);
-      for (Int_t iModule=0; iModule < nOfModules; iModule++) {
-       cpvModule = fPHOS->GetCPVModule(iModule);
-       cpvHits   = cpvModule.Hits();
-       nCPVhits  = cpvHits->GetEntriesFast();
-       for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
-         cpvHit   = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
-         p        = cpvHit->GetMomentum();
-         xzgen[0] = cpvHit->X();
-         xzgen[1] = cpvHit->Y();
-         ipart    = cpvHit->GetIpart();
-         TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
-         new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*cpvHit);
-       }
-       cpvModule.Clear();
-      }
-    }
-    for (iModule=0; iModule < nOfModules; iModule++) {
-      Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
-      printf("Module %d has %d hits\n",iModule,nsum);
-    }
+//     Int_t nOfModules = fGeom->GetNModules();
+//     TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
+//     Int_t iModule = 0;      
+//     for (iModule=0; iModule < nOfModules; iModule++)
+//       hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100);
 
-    // Then go through reconstructed points and for each find
-    // the closeset hit
-    // The distance from the rec.point to the closest hit
-    // gives the coordinate resolution of the CPV
-
-    // Get the Reconstruction Tree
-    gAlice->TreeR()->GetEvent(0) ;
-    TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
-    AliPHOSCpvRecPoint *cpvRecPoint ;
-    Float_t xgen, zgen;
-    while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) {
-      TVector3  locpos;
-      cpvRecPoint->GetLocalPosition(locpos);
-      Int_t phosModule = cpvRecPoint->GetPHOSMod();
-      Int_t rpMult     = cpvRecPoint->GetDigitsMultiplicity();
-      Int_t rpMultX, rpMultZ;
-      cpvRecPoint->GetClusterLengths(rpMultX,rpMultZ);
-      Float_t xrec  = locpos.X();
-      Float_t zrec  = locpos.Z();
-      Float_t dxmin = 1.e+10;
-      Float_t dzmin = 1.e+10;
-      Float_t r2min = 1.e+10;
-      Float_t r2;
-
-      cpvHits = hitsPerModule[phosModule-1];
-      Int_t nCPVhits  = cpvHits->GetEntriesFast();
-      for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
-       cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
-       xgen   = cpvHit->X();
-       zgen   = cpvHit->Y();
-       r2 = TMath::Power((xgen-xrec),2) + TMath::Power((zgen-zrec),2);
-       if ( r2 < r2min ) {
-         r2min = r2;
-         dxmin = xgen - xrec;
-         dzmin = zgen - zrec;
-       }
-      }
-      hDx  ->Fill(dxmin);
-      hDz  ->Fill(dzmin);
-      hDr  ->Fill(TMath::Sqrt(r2min));
-      hNrp ->Fill(rpMult);
-      hNrpX->Fill(rpMultX);
-      hNrpZ->Fill(rpMultZ);
-    }
-    delete [] hitsPerModule;
-  }
-  // Save histograms
+//     AliPHOSCPVModule cpvModule;
+//     TClonesArray    *cpvHits;
+//     Int_t           nCPVhits;
+//     AliPHOSCPVHit   *cpvHit;
+//     TLorentzVector   p;
+//     Float_t          xzgen[2];
+//     Int_t            ipart;
 
-  Text_t outputname[80] ;
-  sprintf(outputname,"%s.analyzed",fRootFile->GetName());
-  TFile output(outputname,"RECREATE");
-  output.cd();
+//     // First go through all primary tracks and fill the arrays
+//     // of hits per each CPV module
 
-  hDx  ->Write() ;
-  hDz  ->Write() ;
-  hDr  ->Write() ;
-  hNrp ->Write() ;
-  hNrpX->Write() ;
-  hNrpZ->Write() ;
-
-  // Plot histograms
-
-  TCanvas *cpvCanvas = new TCanvas("CPV","CPV analysis",20,20,800,400);
-  gStyle->SetOptStat(111111);
-  gStyle->SetOptFit(1);
-  gStyle->SetOptDate(1);
-  cpvCanvas->Divide(3,2);
-
-  cpvCanvas->cd(1);
-  gPad->SetFillColor(10);
-  hNrp->SetFillColor(16);
-  hNrp->Draw();
-
-  cpvCanvas->cd(2);
-  gPad->SetFillColor(10);
-  hNrpX->SetFillColor(16);
-  hNrpX->Draw();
-
-  cpvCanvas->cd(3);
-  gPad->SetFillColor(10);
-  hNrpZ->SetFillColor(16);
-  hNrpZ->Draw();
-
-  cpvCanvas->cd(4);
-  gPad->SetFillColor(10);
-  hDx->SetFillColor(16);
-  hDx->Fit("gaus");
-  hDx->Draw();
-
-  cpvCanvas->cd(5);
-  gPad->SetFillColor(10);
-  hDz->SetFillColor(16);
-  hDz->Fit("gaus");
-  hDz->Draw();
-
-  cpvCanvas->cd(6);
-  gPad->SetFillColor(10);
-  hDr->SetFillColor(16);
-  hDr->Draw();
-
-  cpvCanvas->Print("CPV.ps");
+//     for (Int_t itrack=0; itrack<ntracks; itrack++) {
+//       // Get the Hits Tree for the Primary track itrack
+//       gAlice->ResetHits();
+//       gAlice->TreeH()->GetEvent(itrack);
+//       for (Int_t iModule=0; iModule < nOfModules; iModule++) {
+//     cpvModule = fPHOS->GetCPVModule(iModule);
+//     cpvHits   = cpvModule.Hits();
+//     nCPVhits  = cpvHits->GetEntriesFast();
+//     for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
+//       cpvHit   = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
+//       p        = cpvHit->GetMomentum();
+//       xzgen[0] = cpvHit->X();
+//       xzgen[1] = cpvHit->Y();
+//       ipart    = cpvHit->GetIpart();
+//       TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
+//       new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*cpvHit);
+//     }
+//     cpvModule.Clear();
+//       }
+//     }
+//     for (iModule=0; iModule < nOfModules; iModule++) {
+//       Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
+//       printf("Module %d has %d hits\n",iModule,nsum);
+//     }
+
+//     // Then go through reconstructed points and for each find
+//     // the closeset hit
+//     // The distance from the rec.point to the closest hit
+//     // gives the coordinate resolution of the CPV
+
+//     // Get the Reconstruction Tree
+//     gAlice->TreeR()->GetEvent(0) ;
+//     TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
+//     AliPHOSCpvRecPoint *cpvRecPoint ;
+//     Float_t xgen, zgen;
+//     while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) {
+//       TVector3  locpos;
+//       cpvRecPoint->GetLocalPosition(locpos);
+//       Int_t phosModule = cpvRecPoint->GetPHOSMod();
+//       Int_t rpMult     = cpvRecPoint->GetDigitsMultiplicity();
+//       Int_t rpMultX, rpMultZ;
+//       cpvRecPoint->GetClusterLengths(rpMultX,rpMultZ);
+//       Float_t xrec  = locpos.X();
+//       Float_t zrec  = locpos.Z();
+//       Float_t dxmin = 1.e+10;
+//       Float_t dzmin = 1.e+10;
+//       Float_t r2min = 1.e+10;
+//       Float_t r2;
+
+//       cpvHits = hitsPerModule[phosModule-1];
+//       Int_t nCPVhits  = cpvHits->GetEntriesFast();
+//       for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
+//     cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
+//     xgen   = cpvHit->X();
+//     zgen   = cpvHit->Y();
+//     r2 = TMath::Power((xgen-xrec),2) + TMath::Power((zgen-zrec),2);
+//     if ( r2 < r2min ) {
+//       r2min = r2;
+//       dxmin = xgen - xrec;
+//       dzmin = zgen - zrec;
+//     }
+//       }
+//       hDx  ->Fill(dxmin);
+//       hDz  ->Fill(dzmin);
+//       hDr  ->Fill(TMath::Sqrt(r2min));
+//       hNrp ->Fill(rpMult);
+//       hNrpX->Fill(rpMultX);
+//       hNrpZ->Fill(rpMultZ);
+//     }
+//     delete [] hitsPerModule;
+//   }
+//   // Save histograms
+
+//   Text_t outputname[80] ;
+//   sprintf(outputname,"%s.analyzed",fRootFile->GetName());
+//   TFile output(outputname,"RECREATE");
+//   output.cd();
+
+//   hDx  ->Write() ;
+//   hDz  ->Write() ;
+//   hDr  ->Write() ;
+//   hNrp ->Write() ;
+//   hNrpX->Write() ;
+//   hNrpZ->Write() ;
+
+//   // Plot histograms
+
+//   TCanvas *cpvCanvas = new TCanvas("CPV","CPV analysis",20,20,800,400);
+//   gStyle->SetOptStat(111111);
+//   gStyle->SetOptFit(1);
+//   gStyle->SetOptDate(1);
+//   cpvCanvas->Divide(3,2);
+
+//   cpvCanvas->cd(1);
+//   gPad->SetFillColor(10);
+//   hNrp->SetFillColor(16);
+//   hNrp->Draw();
+
+//   cpvCanvas->cd(2);
+//   gPad->SetFillColor(10);
+//   hNrpX->SetFillColor(16);
+//   hNrpX->Draw();
+
+//   cpvCanvas->cd(3);
+//   gPad->SetFillColor(10);
+//   hNrpZ->SetFillColor(16);
+//   hNrpZ->Draw();
+
+//   cpvCanvas->cd(4);
+//   gPad->SetFillColor(10);
+//   hDx->SetFillColor(16);
+//   hDx->Fit("gaus");
+//   hDx->Draw();
+
+//   cpvCanvas->cd(5);
+//   gPad->SetFillColor(10);
+//   hDz->SetFillColor(16);
+//   hDz->Fit("gaus");
+//   hDz->Draw();
+
+//   cpvCanvas->cd(6);
+//   gPad->SetFillColor(10);
+//   hDr->SetFillColor(16);
+//   hDr->Draw();
+
+//   cpvCanvas->Print("CPV.ps");
 
 }
 
@@ -723,257 +722,257 @@ void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
 //____________________________________________________________________________
  void AliPHOSAnalyze::ReadAndPrintEMC(Int_t EvFirst, Int_t EvLast)    
 {
-  //
-  // Read and print generated and reconstructed hits in EMC
-  // for events from EvFirst to Nevent.
-  // If only EvFirst is defined, print only this one event.
-  // Author: Yuri Kharlov
-  // 24 November 2000
-  //
-
-  if (EvFirst!=0 && EvLast==0) EvLast=EvFirst;
-  Int_t ievent;
-  for (ievent=EvFirst; ievent<=EvLast; ievent++) {  
+//   //
+//   // Read and print generated and reconstructed hits in EMC
+//   // for events from EvFirst to Nevent.
+//   // If only EvFirst is defined, print only this one event.
+//   // Author: Yuri Kharlov
+//   // 24 November 2000
+//   //
+
+//   if (EvFirst!=0 && EvLast==0) EvLast=EvFirst;
+//   Int_t ievent;
+//   for (ievent=EvFirst; ievent<=EvLast; ievent++) {  
     
-    //========== Event Number>
-    cout << endl <<  "==== ReadAndPrintEMC ====> Event is " << ievent+1 << endl ;
+//     //========== Event Number>
+//     cout << endl <<  "==== ReadAndPrintEMC ====> Event is " << ievent+1 << endl ;
 
-    //=========== Connects the various Tree's for evt
-    Int_t ntracks = gAlice->GetEvent(ievent);
-    fPHOS->SetTreeAddress() ;
+//     //=========== Connects the various Tree's for evt
+//     Int_t ntracks = gAlice->GetEvent(ievent);
+//     fPHOS->SetTreeAddress() ;
     
-    gAlice->TreeD()->GetEvent(0) ;
-    gAlice->TreeR()->GetEvent(0) ;
+//     gAlice->TreeD()->GetEvent(0) ;
+//     gAlice->TreeR()->GetEvent(0) ;
 
-    // Loop over reconstructed particles
+//     // Loop over reconstructed particles
       
-    TClonesArray ** recParticleList  = fPHOS->RecParticles() ;     
-    AliPHOSRecParticle * recParticle ;
-    Int_t iRecParticle ;
-    Int_t *primList;
-    Int_t nPrimary;
-    for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ ) {
-      recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
-      Float_t recE = recParticle->Energy();
-      primList     = recParticle->GetPrimaries(nPrimary);
-      Int_t moduleNumberRec ;
-      Double_t recX, recZ ;
-      fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
-      printf("Rec point: module %d, (X,Z) = (%8.4f,%8.4f) cm, E = %.3f GeV, primary = %d\n",
-            moduleNumberRec,recX,recZ,recE,*primList);
-    }
+//     TClonesArray ** recParticleList  = fPHOS->RecParticles() ;     
+//     AliPHOSRecParticle * recParticle ;
+//     Int_t iRecParticle ;
+//     Int_t *primList;
+//     Int_t nPrimary;
+//     for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ ) {
+//       recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
+//       Float_t recE = recParticle->Energy();
+//       primList     = recParticle->GetPrimaries(nPrimary);
+//       Int_t moduleNumberRec ;
+//       Double_t recX, recZ ;
+//       fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
+//       printf("Rec point: module %d, (X,Z) = (%8.4f,%8.4f) cm, E = %.3f GeV, primary = %d\n",
+//          moduleNumberRec,recX,recZ,recE,*primList);
+//     }
 
-    // Read and print EMC hits from EMCn branches
+//     // Read and print EMC hits from EMCn branches
       
-    AliPHOSCPVModule emcModule;
-    TClonesArray    *emcHits;
-    Int_t           nEMChits;
-    AliPHOSCPVHit   *emcHit;
-    TLorentzVector   p;
-    Float_t          xgen, zgen;
-    Int_t            ipart, primary;
-    Int_t            nGenHits = 0;
-    for (Int_t itrack=0; itrack<ntracks; itrack++) {
-      //=========== Get the Hits Tree for the Primary track itrack
-      gAlice->ResetHits();
-      gAlice->TreeH()->GetEvent(itrack);
-      Int_t iModule = 0 ;
-      for (iModule=0; iModule < fGeom->GetNModules(); iModule++) {
-       emcModule = fPHOS->GetEMCModule(iModule);
-       emcHits   = emcModule.Hits();
-       nEMChits  = emcHits->GetEntriesFast();
-       for (Int_t ihit=0; ihit<nEMChits; ihit++) {
-         nGenHits++;
-         emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
-         p      = emcHit->GetMomentum();
-         xgen   = emcHit->X();
-         zgen   = emcHit->Y();
-         ipart  = emcHit->GetIpart();
-         primary= emcHit->GetTrack();
-         printf("EMC hit A: module %d, ",iModule+1);
-         printf("    p = (%f .4, %f .4, %f .4, %f .4) GeV,\n",
-                p.Px(),p.Py(),p.Pz(),p.Energy());
-         printf("                     (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
-                xgen,zgen,ipart,primary);
-       }
-      }
-    }
+//     AliPHOSCPVModule emcModule;
+//     TClonesArray    *emcHits;
+//     Int_t           nEMChits;
+//     AliPHOSCPVHit   *emcHit;
+//     TLorentzVector   p;
+//     Float_t          xgen, zgen;
+//     Int_t            ipart, primary;
+//     Int_t            nGenHits = 0;
+//     for (Int_t itrack=0; itrack<ntracks; itrack++) {
+//       //=========== Get the Hits Tree for the Primary track itrack
+//       gAlice->ResetHits();
+//       gAlice->TreeH()->GetEvent(itrack);
+//       Int_t iModule = 0 ;
+//       for (iModule=0; iModule < fGeom->GetNModules(); iModule++) {
+//     emcModule = fPHOS->GetEMCModule(iModule);
+//     emcHits   = emcModule.Hits();
+//     nEMChits  = emcHits->GetEntriesFast();
+//     for (Int_t ihit=0; ihit<nEMChits; ihit++) {
+//       nGenHits++;
+//       emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
+//       p      = emcHit->GetMomentum();
+//       xgen   = emcHit->X();
+//       zgen   = emcHit->Y();
+//       ipart  = emcHit->GetIpart();
+//       primary= emcHit->GetTrack();
+//       printf("EMC hit A: module %d, ",iModule+1);
+//       printf("    p = (%f .4, %f .4, %f .4, %f .4) GeV,\n",
+//              p.Px(),p.Py(),p.Pz(),p.Energy());
+//       printf("                     (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
+//              xgen,zgen,ipart,primary);
+//     }
+//       }
+//     }
 
-//      // Read and print EMC hits from PHOS branch
-
-//      for (Int_t itrack=0; itrack<ntracks; itrack++) {
-//        //=========== Get the Hits Tree for the Primary track itrack
-//        gAlice->ResetHits();
-//        gAlice->TreeH()->GetEvent(itrack);
-//        TClonesArray *hits = fPHOS->Hits();
-//        AliPHOSHit   *hit ;
-//        Int_t ihit;
-//        for ( ihit = 0 ; ihit < hits->GetEntries() ; ihit++ ) {
-//     hit = (AliPHOSHit*)hits->At(ihit) ;
-//     Float_t hitXYZ[3];
-//     hitXYZ[0]   = hit->X();
-//     hitXYZ[1]   = hit->Y();
-//     hitXYZ[2]   = hit->Z();
-//     ipart       = hit->GetPid();
-//     primary     = hit->GetPrimary();
-//     Int_t absId = hit->GetId();
-//     Int_t relId[4];
-//     fGeom->AbsToRelNumbering(absId, relId) ;
-//     Int_t module = relId[0];
-//     if (relId[1]==0 && !(hitXYZ[0]==0 && hitXYZ[2]==0))
-//       printf("EMC hit B: module %d, (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
-//              module,hitXYZ[0],hitXYZ[2],ipart,primary);
-//        }
-//      }
+// //      // Read and print EMC hits from PHOS branch
+
+// //      for (Int_t itrack=0; itrack<ntracks; itrack++) {
+// //        //=========== Get the Hits Tree for the Primary track itrack
+// //        gAlice->ResetHits();
+// //        gAlice->TreeH()->GetEvent(itrack);
+// //        TClonesArray *hits = fPHOS->Hits();
+// //        AliPHOSHit   *hit ;
+// //        Int_t ihit;
+// //        for ( ihit = 0 ; ihit < hits->GetEntries() ; ihit++ ) {
+// //          hit = (AliPHOSHit*)hits->At(ihit) ;
+// //          Float_t hitXYZ[3];
+// //          hitXYZ[0]   = hit->X();
+// //          hitXYZ[1]   = hit->Y();
+// //          hitXYZ[2]   = hit->Z();
+// //          ipart       = hit->GetPid();
+// //          primary     = hit->GetPrimary();
+// //          Int_t absId = hit->GetId();
+// //          Int_t relId[4];
+// //          fGeom->AbsToRelNumbering(absId, relId) ;
+// //          Int_t module = relId[0];
+// //          if (relId[1]==0 && !(hitXYZ[0]==0 && hitXYZ[2]==0))
+// //            printf("EMC hit B: module %d, (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
+// //                   module,hitXYZ[0],hitXYZ[2],ipart,primary);
+// //        }
+// //      }
 
-  }
+//   }
 }
 
 //____________________________________________________________________________
  void AliPHOSAnalyze::AnalyzeEMC(Int_t Nevents)
 {
-  //
-  // Read generated and reconstructed hits in EMC for Nevents events.
-  // Plots the coordinate and energy resolution histograms.
-  // Coordinate resolution is a difference between the reconstructed
-  // coordinate and the exact coordinate on the face of the PHOS
-  // Author: Yuri Kharlov
-  // 27 November 2000
-  //
-
-  // Book histograms
-
-  TH1F *hDx1   = new TH1F("hDx1"  ,"EMC x-resolution", 100,-5. , 5.);
-  TH1F *hDz1   = new TH1F("hDz1"  ,"EMC z-resolution", 100,-5. , 5.);
-  TH1F *hDE1   = new TH1F("hDE1"  ,"EMC E-resolution", 100,-2. , 2.);
-
-  TH2F *hDx2   = new TH2F("hDx2"  ,"EMC x-resolution", 100, 0., 10., 100,-5. , 5.);
-  TH2F *hDz2   = new TH2F("hDz2"  ,"EMC z-resolution", 100, 0., 10., 100,-5. , 5.);
-  TH2F *hDE2   = new TH2F("hDE2"  ,"EMC E-resolution", 100, 0., 10., 100, 0. , 5.);
-
-  cout << "Start EMC Analysis"<< endl ;
-  for (Int_t ievent=0; ievent<Nevents; ievent++) {  
+//   //
+//   // Read generated and reconstructed hits in EMC for Nevents events.
+//   // Plots the coordinate and energy resolution histograms.
+//   // Coordinate resolution is a difference between the reconstructed
+//   // coordinate and the exact coordinate on the face of the PHOS
+//   // Author: Yuri Kharlov
+//   // 27 November 2000
+//   //
+
+//   // Book histograms
+
+//   TH1F *hDx1   = new TH1F("hDx1"  ,"EMC x-resolution", 100,-5. , 5.);
+//   TH1F *hDz1   = new TH1F("hDz1"  ,"EMC z-resolution", 100,-5. , 5.);
+//   TH1F *hDE1   = new TH1F("hDE1"  ,"EMC E-resolution", 100,-2. , 2.);
+
+//   TH2F *hDx2   = new TH2F("hDx2"  ,"EMC x-resolution", 100, 0., 10., 100,-5. , 5.);
+//   TH2F *hDz2   = new TH2F("hDz2"  ,"EMC z-resolution", 100, 0., 10., 100,-5. , 5.);
+//   TH2F *hDE2   = new TH2F("hDE2"  ,"EMC E-resolution", 100, 0., 10., 100, 0. , 5.);
+
+//   cout << "Start EMC Analysis"<< endl ;
+//   for (Int_t ievent=0; ievent<Nevents; ievent++) {  
       
-    //========== Event Number>         
-    if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
-      cout << "==== AnalyzeEMC ====> Event is " << ievent+1 << endl ;
+//     //========== Event Number>         
+//     if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
+//       cout << "==== AnalyzeEMC ====> Event is " << ievent+1 << endl ;
     
-    //=========== Connects the various Tree's for evt
-    Int_t ntracks = gAlice->GetEvent(ievent);
+//     //=========== Connects the various Tree's for evt
+//     Int_t ntracks = gAlice->GetEvent(ievent);
 
-    fPHOS->SetTreeAddress() ;
+//     fPHOS->SetTreeAddress() ;
     
-    gAlice->TreeD()->GetEvent(0) ;
-    gAlice->TreeR()->GetEvent(0) ;
+//     gAlice->TreeD()->GetEvent(0) ;
+//     gAlice->TreeR()->GetEvent(0) ;
 
-    // Create and fill arrays of hits for each EMC module
+//     // Create and fill arrays of hits for each EMC module
       
-    Int_t nOfModules = fGeom->GetNModules();
-    TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
-    Int_t iModule;
-    for (iModule=0; iModule < nOfModules; iModule++)
-      hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100);
-
-    AliPHOSCPVModule emcModule;
-    TClonesArray    *emcHits;
-    Int_t           nEMChits;
-    AliPHOSCPVHit   *emcHit;
-
-    // First go through all primary tracks and fill the arrays
-    // of hits per each EMC module
-
-    for (Int_t itrack=0; itrack<ntracks; itrack++) {
-      // Get the Hits Tree for the Primary track itrack
-      gAlice->ResetHits();
-      gAlice->TreeH()->GetEvent(itrack);
-      for (Int_t iModule=0; iModule < nOfModules; iModule++) {
-       emcModule = fPHOS->GetEMCModule(iModule);
-       emcHits   = emcModule.Hits();
-       nEMChits  = emcHits->GetEntriesFast();
-       for (Int_t ihit=0; ihit<nEMChits; ihit++) {
-         emcHit   = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
-         TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
-         new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*emcHit);
-       }
-       emcModule.Clear();
-      }
-    }
+//     Int_t nOfModules = fGeom->GetNModules();
+//     TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
+//     Int_t iModule;
+//     for (iModule=0; iModule < nOfModules; iModule++)
+//       hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100);
 
-    // Loop over reconstructed particles
-      
-    TClonesArray ** recParticleList  = fPHOS->RecParticles() ;     
-    AliPHOSRecParticle * recParticle ;
-    Int_t nEMCrecs = (*recParticleList)->GetEntries();
-    if (nEMCrecs == 1) {
-      recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(0) ;
-      Float_t recE = recParticle->Energy();
-      Int_t phosModule;
-      Double_t recX, recZ ;
-      fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), phosModule, recX, recZ) ;
-
-      // for this rec.point take the hit list in the same PHOS module
-
-      emcHits = hitsPerModule[phosModule-1];
-      Int_t nEMChits  = emcHits->GetEntriesFast();
-      if (nEMChits == 1) {
-       Float_t genX, genZ, genE;
-       for (Int_t ihit=0; ihit<nEMChits; ihit++) {
-         emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
-         genX   = emcHit->X();
-         genZ   = emcHit->Y();
-         genE   = emcHit->GetMomentum().E();
-       }
-       Float_t dx = recX - genX;
-       Float_t dz = recZ - genZ;
-       Float_t de = recE - genE;
-       hDx1  ->Fill(dx);
-       hDz1  ->Fill(dz);
-       hDE1  ->Fill(de);
-       hDx2  ->Fill(genE,dx);
-       hDz2  ->Fill(genE,dz);
-       hDE2  ->Fill(genE,recE);
-      }
-    }
-    delete [] hitsPerModule;
-  }
-  // Save histograms
+//     AliPHOSCPVModule emcModule;
+//     TClonesArray    *emcHits;
+//     Int_t           nEMChits;
+//     AliPHOSCPVHit   *emcHit;
 
-  Text_t outputname[80] ;
-  sprintf(outputname,"%s.analyzed",fRootFile->GetName());
-  TFile output(outputname,"RECREATE");
-  output.cd();
+//     // First go through all primary tracks and fill the arrays
+//     // of hits per each EMC module
 
-  hDx1  ->Write() ;
-  hDz1  ->Write() ;
-  hDE1  ->Write() ;
-  hDx2  ->Write() ;
-  hDz2  ->Write() ;
-  hDE2  ->Write() ;
-
-  // Plot histograms
-
-  TCanvas *emcCanvas = new TCanvas("EMC","EMC analysis",20,20,700,300);
-  gStyle->SetOptStat(111111);
-  gStyle->SetOptFit(1);
-  gStyle->SetOptDate(1);
-  emcCanvas->Divide(3,1);
-
-  emcCanvas->cd(1);
-  gPad->SetFillColor(10);
-  hDx1->SetFillColor(16);
-  hDx1->Draw();
-
-  emcCanvas->cd(2);
-  gPad->SetFillColor(10);
-  hDz1->SetFillColor(16);
-  hDz1->Draw();
-
-  emcCanvas->cd(3);
-  gPad->SetFillColor(10);
-  hDE1->SetFillColor(16);
-  hDE1->Draw();
-
-  emcCanvas->Print("EMC.ps");
+//     for (Int_t itrack=0; itrack<ntracks; itrack++) {
+//       // Get the Hits Tree for the Primary track itrack
+//       gAlice->ResetHits();
+//       gAlice->TreeH()->GetEvent(itrack);
+//       for (Int_t iModule=0; iModule < nOfModules; iModule++) {
+//     emcModule = fPHOS->GetEMCModule(iModule);
+//     emcHits   = emcModule.Hits();
+//     nEMChits  = emcHits->GetEntriesFast();
+//     for (Int_t ihit=0; ihit<nEMChits; ihit++) {
+//       emcHit   = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
+//       TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
+//       new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*emcHit);
+//     }
+//     emcModule.Clear();
+//       }
+//     }
+
+//     // Loop over reconstructed particles
+      
+//     TClonesArray ** recParticleList  = fPHOS->RecParticles() ;     
+//     AliPHOSRecParticle * recParticle ;
+//     Int_t nEMCrecs = (*recParticleList)->GetEntries();
+//     if (nEMCrecs == 1) {
+//       recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(0) ;
+//       Float_t recE = recParticle->Energy();
+//       Int_t phosModule;
+//       Double_t recX, recZ ;
+//       fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), phosModule, recX, recZ) ;
+
+//       // for this rec.point take the hit list in the same PHOS module
+
+//       emcHits = hitsPerModule[phosModule-1];
+//       Int_t nEMChits  = emcHits->GetEntriesFast();
+//       if (nEMChits == 1) {
+//     Float_t genX, genZ, genE;
+//     for (Int_t ihit=0; ihit<nEMChits; ihit++) {
+//       emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
+//       genX   = emcHit->X();
+//       genZ   = emcHit->Y();
+//       genE   = emcHit->GetMomentum().E();
+//     }
+//     Float_t dx = recX - genX;
+//     Float_t dz = recZ - genZ;
+//     Float_t de = recE - genE;
+//     hDx1  ->Fill(dx);
+//     hDz1  ->Fill(dz);
+//     hDE1  ->Fill(de);
+//     hDx2  ->Fill(genE,dx);
+//     hDz2  ->Fill(genE,dz);
+//     hDE2  ->Fill(genE,recE);
+//       }
+//     }
+//     delete [] hitsPerModule;
+//   }
+//   // Save histograms
+
+//   Text_t outputname[80] ;
+//   sprintf(outputname,"%s.analyzed",fRootFile->GetName());
+//   TFile output(outputname,"RECREATE");
+//   output.cd();
+
+//   hDx1  ->Write() ;
+//   hDz1  ->Write() ;
+//   hDE1  ->Write() ;
+//   hDx2  ->Write() ;
+//   hDz2  ->Write() ;
+//   hDE2  ->Write() ;
+
+//   // Plot histograms
+
+//   TCanvas *emcCanvas = new TCanvas("EMC","EMC analysis",20,20,700,300);
+//   gStyle->SetOptStat(111111);
+//   gStyle->SetOptFit(1);
+//   gStyle->SetOptDate(1);
+//   emcCanvas->Divide(3,1);
+
+//   emcCanvas->cd(1);
+//   gPad->SetFillColor(10);
+//   hDx1->SetFillColor(16);
+//   hDx1->Draw();
+
+//   emcCanvas->cd(2);
+//   gPad->SetFillColor(10);
+//   hDz1->SetFillColor(16);
+//   hDz1->Draw();
+
+//   emcCanvas->cd(3);
+//   gPad->SetFillColor(10);
+//   hDE1->SetFillColor(16);
+//   hDE1->Draw();
+
+//   emcCanvas->Print("EMC.ps");
 
 }
 
index c137b5896e74c1e2193692eaffb9b0d2cb167351..7ad00684d193f5de5236f687635b530f34220175 100644 (file)
@@ -95,7 +95,7 @@ AliPHOSDigit::~AliPHOSDigit()
 }
 
 //____________________________________________________________________________
-Int_t AliPHOSDigit::Compare(TObject * obj)
+Int_t AliPHOSDigit::Compare(const TObject * obj) const
 {
   // Compares two digits with respect to its Id
   // to sort according increasing Id
index 5d4c5ccb2cc0fc9394dbac9ecd7d028072d0549b..c8bafc6483b8b1e7b0fbb589b397afc9f7042aff 100644 (file)
@@ -39,7 +39,7 @@ class AliPHOSDigit : public AliDigitNew {
   Bool_t operator==(const AliPHOSDigit &rValue) const;
   AliPHOSDigit& operator+(AliPHOSDigit const &rValue) ;
     
-  Int_t   Compare(TObject * obj) ;  
+  Int_t   Compare(const TObject * obj) const ;  
   Int_t   GetNprimary() const { 
     // returns the number of primaries
     return fNprimary ; }
index 762689547b02a6cd98dee7a8f826cfea50a2d865..669553ea341b777c1234433ee65c95eb8cf0506c 100644 (file)
@@ -54,27 +54,30 @@ AliPHOSHit::AliPHOSHit(const AliPHOSHit & hit)
 } 
 
 //____________________________________________________________________________
-AliPHOSHit::AliPHOSHit(Int_t Shunt, Int_t primary, Int_t Track, Int_t id, Float_t *hits, Int_t pid) : AliHit(Shunt, Track)
+AliPHOSHit::AliPHOSHit(Int_t shunt, Int_t primary, Int_t track, Int_t id, Float_t *hits, Int_t pid, TLorentzVector p, Float_t *xy): AliHit(shunt, track)
 {
-  // ctor
-  
-   fId         = id ;
-   fX          = hits[0] ;
-   fY          = hits[1] ;
-   fZ          = hits[2] ;
-   fELOS       = hits[3] ;
-   fPrimary    = primary ;
-   fPid        = pid ; 
+  //
+  // Create a CPV hit object
+  //
+
+  fId         = id ;
+  fELOS       = hits[3] ;
+  fPrimary    = primary ;
+  fPid        = pid ; 
+  fMomentum  = p;
+  fX         = xy[0];  //position of particle first entering cristall/pad
+  fY         = xy[1];
+  fZ         = 9999.;  //Fake Z to avoid FPE
 }
 
 //____________________________________________________________________________
 Bool_t AliPHOSHit::operator==(AliPHOSHit const &rValue) const
 { 
-  // Two hits are identical if they have the same Id and originate from the same primary
+  // Two hits are identical if they have the same Id and originate from the same primary 
 
   Bool_t rv = kFALSE ; 
 
-  if ( fId == rValue.GetId() && fPrimary == rValue.GetPrimary() )
+  if ( (fId == rValue.GetId()) && ( fPrimary == rValue.GetPrimary() ) && (fPid*rValue.fPid ! = 0) )
     rv = kTRUE;
   
   return rv;
@@ -94,6 +97,8 @@ AliPHOSHit AliPHOSHit::operator+(const AliPHOSHit &rValue) const
 
    added.fELOS += rValue.GetEnergy() ;
     
+   if(fPid == 0) fPid = rValue.fPid ;
+
    return added;
 
 }
index 936da55fcc82d59db4d986d8857cd429ab812049..f2b0beddc92d4513094b8914151f1ae3a66cde51 100644 (file)
@@ -20,6 +20,8 @@
 
 #include <iostream.h>
 
+class TLorentzVector ;
+
 class AliPHOSHit : public AliHit {
 
   friend ostream& operator << (ostream&, const AliPHOSHit&) ;
@@ -30,7 +32,7 @@ class AliPHOSHit : public AliHit {
     // default ctor 
   }
   AliPHOSHit(const AliPHOSHit & hit) ; 
-  AliPHOSHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t id, Float_t *hits, Int_t pid;
+  AliPHOSHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t id, Float_t *hits, Int_t pid, TLorentzVector p, Float_t *xy);
   virtual ~AliPHOSHit(void) {
     // dtor 
   }  
@@ -51,6 +53,8 @@ class AliPHOSHit : public AliHit {
     // returns the primary particle id at the origine of this hit 
     return fPrimary ; 
   }
+  TLorentzVector GetMomentum()  { return  fMomentum; }   
+    // momentum of the particle which initiated this hit
 
   Bool_t operator == (AliPHOSHit const &rValue) const ;
   AliPHOSHit operator + (const AliPHOSHit& rValue) const ;
@@ -58,10 +62,12 @@ class AliPHOSHit : public AliHit {
 
  private:
 
-  Int_t     fId ;        // Absolute Id number of PHOS Xtal or PPSD pad
-  Float_t   fELOS ;      // Energy deposited
-  Int_t     fPid ;       // type of the particle that initiates that hit 
-  Int_t     fPrimary ;   // Primary particles at the origine of the hit
+  Int_t          fId ;        // Absolute Id number of PHOS Xtal or PPSD pad
+  Float_t        fELOS ;      // Energy deposited
+  Int_t          fPid ;       // type of the particle that initiates that hit 
+  Int_t          fPrimary ;   // Primary particles at the origine of the hit
+
+  TLorentzVector fMomentum;   // 4-momentum of the particle
 
   ClassDef(AliPHOSHit,1)  // Hit for PHOS
 
index 30a4bdf97814c885dfca31b753b419b43e9a8d8e..6e50ea0da962093d4b24d3cbc30ab6e143067974 100644 (file)
@@ -67,16 +67,18 @@ AliPHOSv1::AliPHOSv1()
  
   fReconstructioner  = 0;
   fTrackSegmentMaker = 0;
+
+  fHits = new TClonesArray("AliPHOSHit",1000) ;
  
-  if ( 0==(fEMCModules=new TClonesArray("AliPHOSCPVModule",0)) ) {
-    Error("AliPHOSv1","Can not create array of EMC modules");
-    exit(1);
-  }
+  //  if ( 0==(fEMCModules=new TClonesArray("AliPHOSCPVModule",0)) ) {
+  //    Error("AliPHOSv1","Can not create array of EMC modules");
+  //    exit(1);
+  //  }
 
-  if ( 0==(fCPVModules=new TClonesArray("AliPHOSCPVModule",0)) ) {
-    Error("AliPHOSv1","Can not create array of CPV modules");
-    exit(1);
-  }
+  //  if ( 0==(fCPVModules=new TClonesArray("AliPHOSCPVModule",0)) ) {
+  //    Error("AliPHOSv1","Can not create array of CPV modules");
+  //    exit(1);
+  //  }
 
 }
 
@@ -96,7 +98,10 @@ AliPHOSv0(name,title)
   //
 
   fPinElectronicNoise = 0.010 ;
-  fDigitThreshold      = 0.1 ;   // 1 GeV 
+  fDigitThreshold      = 0.01 ;   // 1 GeV 
+  fDigitizeA= 0. ;             
+  fDigitizeB = 10000000. ;    
+
 
   // We do not want to save in TreeH the raw hits
   // But save the cumulated hits instead (need to create the branch myself)
@@ -114,32 +119,32 @@ AliPHOSv0(name,title)
   
   // Create array of EMC modules of the size of PHOS modules number
   
-  if ( 0==(fEMCModules=new TClonesArray("AliPHOSCPVModule",fGeom->GetNModules())) ) {
-    Error("AliPHOSv1","Can not create array of EMC modules");
-    exit(1);
-  }
-  TClonesArray &lemcmodule = *fEMCModules;
-  for (Int_t i=0; i<fGeom->GetNModules(); i++) new(lemcmodule[i]) AliPHOSCPVModule();
+//  if ( 0==(fEMCModules=new TClonesArray("AliPHOSCPVModule",fGeom->GetNModules())) ) {
+//    Error("AliPHOSv1","Can not create array of EMC modules");
+//    exit(1);
+//  }
+//  TClonesArray &lemcmodule = *fEMCModules;
+//  for (Int_t i=0; i<fGeom->GetNModules(); i++) new(lemcmodule[i]) AliPHOSCPVModule();
 
   // Create array of CPV modules for the IHEP's version of CPV
 
-  if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
-    // Create array of CPV modules of the size of PHOS modules number
+//  if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
+//    // Create array of CPV modules of the size of PHOS modules number
 
-    if ( 0==(fCPVModules=new TClonesArray("AliPHOSCPVModule",fGeom->GetNCPVModules())) ) {
-      Error("AliPHOSv1","Can not create array of CPV modules");
-      exit(1);
-    }
-    TClonesArray &lcpvmodule = *fCPVModules;
-    for (Int_t i=0; i<fGeom->GetNCPVModules(); i++) new(lcpvmodule[i]) AliPHOSCPVModule();
-  }
-  else {
-    // Create an empty array of AliPHOSCPVModule to satisfy
-    // AliPHOSv1::Streamer when writing root file
+//    if ( 0==(fCPVModules=new TClonesArray("AliPHOSCPVModule",fGeom->GetNCPVModules())) ) {
+//      Error("AliPHOSv1","Can not create array of CPV modules");
+//      exit(1);
+//    }
+//    TClonesArray &lcpvmodule = *fCPVModules;
+//    for (Int_t i=0; i<fGeom->GetNCPVModules(); i++) new(lcpvmodule[i]) AliPHOSCPVModule();
+//  }
+//  else {
+//    // Create an empty array of AliPHOSCPVModule to satisfy
+//    // AliPHOSv1::Streamer when writing root file
 
-    fCPVModules=new TClonesArray("AliPHOSCPVModule",0);
+//    fCPVModules=new TClonesArray("AliPHOSCPVModule",0);
 
-  }
+//  }
 }
 
 //____________________________________________________________________________
@@ -194,6 +199,12 @@ AliPHOSv1::~AliPHOSv1()
     fHits = 0 ; 
   }
 
+  if ( fSDigits) {
+    fSDigits->Delete() ; 
+    delete fSDigits ;
+    fSDigits = 0 ; 
+  }
+
   if ( fDigits) {
     fDigits->Delete() ; 
     delete fDigits ;
@@ -221,7 +232,7 @@ AliPHOSv1::~AliPHOSv1()
 }
 
 //____________________________________________________________________________
-void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id, Float_t * hits, Int_t trackpid)
+void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id, Float_t * hits, Int_t trackpid, TLorentzVector p, Float_t * lpos)
 {
   // Add a hit to the hit list.
   // A PHOS hit is the sum of all hits in a single crystal
@@ -232,7 +243,7 @@ void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id,
   AliPHOSHit *curHit ;
   Bool_t deja = kFALSE ;
 
-  newHit = new AliPHOSHit(shunt, primary, tracknumber, Id, hits, trackpid) ;
+  newHit = new AliPHOSHit(shunt, primary, tracknumber, Id, hits, trackpid, p, lpos) ;
 
   for ( hitCounter = 0 ; hitCounter < fNhits && !deja ; hitCounter++ ) {
     curHit = (AliPHOSHit*) (*fHits)[hitCounter] ;
@@ -250,49 +261,26 @@ void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id,
   delete newHit;
 }
 
-//___________________________________________________________________________
-Int_t AliPHOSv1::Digitize(Float_t Energy)
-{
-  // Applies the energy calibration
-  
-  Float_t fB = 10000000. ;
-  Float_t fA = 0. ;
-  Int_t chan = Int_t(fA + Energy*fB ) ;
-  return chan ;
-}
-
 //____________________________________________________________________________
-void AliPHOSv1::Hit2Digit(Int_t ntracks){
-  //Collects all hits in the same active volume into digits
-  
-  if(fDigits!= 0)
-    fDigits->Clear() ;
-  else
-    fDigits = new TClonesArray("AliPHOSDigit",1000) ;
-  
-  // Branch address for digit tree
-  char branchname[20];
-  sprintf(branchname,"%s",GetName());
-  gAlice->TreeD()->Branch(branchname,&fDigits,fBufferSize);  
-  
-  gAlice->TreeD()->GetEvent(0);
+void AliPHOSv1::Hits2SDigits(){
+  //Collects all hits in the same active volume into digit
 
-  
   Int_t i ;
-  Int_t relid[4];
   Int_t j ; 
   AliPHOSHit   * hit ;
   AliPHOSDigit * newdigit ;
   AliPHOSDigit * curdigit ;
   Bool_t deja = kFALSE ; 
   
+
   Int_t itrack ;
-  for (itrack=0; itrack<ntracks; itrack++){
-    
+  for (itrack=0; itrack<gAlice->GetNtrack(); itrack++){
+        
     //=========== Get the Hits Tree for the Primary track itrack
     gAlice->ResetHits();    
     gAlice->TreeH()->GetEvent(itrack);
       
+
     for ( i = 0 ; i < fHits->GetEntries() ; i++ ) {
       hit = (AliPHOSHit*)fHits->At(i) ;
     
@@ -303,9 +291,8 @@ void AliPHOSv1::Hit2Digit(Int_t ntracks){
        newdigit = new AliPHOSDigit( -1               , hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
       deja =kFALSE ;
 
-
-      for ( j = 0 ; j < fNdigits ;  j++) { 
-       curdigit = (AliPHOSDigit*) fDigits->At(j) ;
+      for ( j = 0 ; j < fnSdigits ;  j++) { 
+       curdigit = (AliPHOSDigit*) fSDigits->At(j) ;
        if ( *curdigit == *newdigit) {
          *curdigit = *curdigit + *newdigit ; 
          deja = kTRUE ; 
@@ -313,44 +300,122 @@ void AliPHOSv1::Hit2Digit(Int_t ntracks){
       }
 
       if ( !deja ) {
-       new((*fDigits)[fNdigits]) AliPHOSDigit(* newdigit) ;
-       fNdigits++ ;  
+       new((*fSDigits)[fnSdigits]) AliPHOSDigit(* newdigit) ;
+       fnSdigits++ ;  
       }
       
       delete newdigit ;    
     } 
 
   } // loop over tracks
-    
-  // Noise induced by the PIN diode of the PbWO crystals
+
+  fSDigits->Sort() ;
+
+  fnSdigits = fSDigits->GetEntries() ;
+  fSDigits->Expand(fnSdigits) ;
+
+  for (i = 0 ; i < fnSdigits ; i++) { 
+    AliPHOSDigit * digit = (AliPHOSDigit *) fSDigits->At(i) ; 
+    digit->SetIndexInList(i) ;     
+  }
+
+  gAlice->TreeS()->Fill() ;
   
-  Float_t energyandnoise ;
-  for ( i = 0 ; i < fNdigits ; i++ ) {
-    newdigit =  (AliPHOSDigit * ) fDigits->At(i) ;
+  TTree *ts = gAlice->TreeS();
+  TDirectory *dir = ts->GetDirectory() ; 
+  cout << "dir name is " << dir->GetName() << endl ;  
+  //  ts->Write(0,TObject::kOverwrite) ;
+
+gAlice->TreeS()->GetBranch("PHOS")->Fill();
+gAlice->TreeS()->GetBranch("PHOS")->Write();
+
+
+}
+//____________________________________________________________________________
+void AliPHOSv1::SDigits2Digits(){
+  //Adds noise to the summable digits and removes everething below thresholds
+  //Note, that sDigits should be SORTED in accordance with abs ID.
+
+
+  gAlice->TreeS()->GetEvent(0) ;
+
+  cout << "fSdigits " <<  fSDigits << " " << fSDigits->GetEntries() << endl ;
+  
+
+  // Noise induced by the PIN diode of the PbWO crystals
+  Int_t iCurSDigit = 0 ;
+  //we assume, that there is al least one EMC digit...
+  Int_t idCurSDigit = ((AliPHOSDigit *)fSDigits->At(0))->GetId() ;
+
+  cout << "fDigits " << fDigits << " " <<idCurSDigit<<  endl ;
+
+  Int_t absID ;
+  for(absID = 1; absID < fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ(); absID++){
     
-    fGeom->AbsToRelNumbering(newdigit->GetId(), relid) ;
+    cout << "absID " << absID << " " << idCurSDigit << endl ; 
+
+    Float_t noise = gRandom->Gaus(0., fPinElectronicNoise) ; 
+    if(absID < idCurSDigit ){ 
+      cout << "In < idC  " << noise << " " << fDigitThreshold << endl ;
+      if(noise >fDigitThreshold ){
+       cout << "noise " << absID << " " << noise << endl;
+       new((*fDigits)[fNdigits]) AliPHOSDigit( -1,absID,Digitize(noise) ) ;
+       cout << "FN " << fNdigits << endl ;
+       fNdigits++ ;  
+      }
+    }
+    else{ //add noise and may be remove the true hit
+      cout << "correcting digit " << iCurSDigit << endl ;
+      Float_t signal = noise + Calibrate(((AliPHOSDigit *)fSDigits->At(iCurSDigit))->GetAmp()) ;
+      cout << "signal " << signal << endl ;
+      if( signal >fDigitThreshold ){
+       cout << "signal " << signal << endl ;
+       AliPHOSDigit * digit = (AliPHOSDigit*) fSDigits->At(iCurSDigit) ;
+       new((*fDigits)[fNdigits]) AliPHOSDigit( *digit ) ;
+       ((AliPHOSDigit *)fDigits->At(fNdigits))->SetAmp(Digitize(signal));
+       cout << "fNdigits " << fNdigits << endl ;
+       fNdigits++ ;  
+      }
+
+      if(iCurSDigit < fSDigits->GetEntries()-1){
+       iCurSDigit++ ;
+       idCurSDigit = ((AliPHOSDigit*)fSDigits->At(iCurSDigit))->GetId() ;
+      }
+      else
+       idCurSDigit = 10000000; //no real hits left    
+    }
     
-    if (relid[1]==0){   // Digits belong to EMC (PbW0_4 crystals)
-      energyandnoise = newdigit->GetAmp() + Digitize(gRandom->Gaus(0., fPinElectronicNoise)) ;
-      
-      if (energyandnoise < 0 ) 
-       energyandnoise = 0 ;
-       
-      if ( newdigit->GetAmp() < fDigitThreshold ) // if threshold not surpassed, remove digit from list
-       fDigits->RemoveAt(i) ; 
+  }  
+
+  //remove PPSD/CPV digits below thresholds
+  Int_t idigit ;
+  for(idigit = iCurSDigit; idigit < fSDigits->GetEntries() ; idigit++){  //loop over CPV/PPSD digits
+    
+    AliPHOSDigit * digit = (AliPHOSDigit *) fSDigits->At(idigit) ; 
+    Float_t ene = Calibrate(digit->GetAmp()) ;
+    
+    Int_t relid[4] ; 
+    fGeom->AbsToRelNumbering(digit->GetId(), relid) ; 
+    if ( relid[0] > fGeom->GetNCPVModules() ){ //ppsd
+      if ( ( (relid[1] > 0) && (ene > fPpsdEnergyThreshold)) ||    //PPSD digit
+          ( (relid[1] < 0) && (ene > fCpvEnergyThreshold ) ) )    //CPV digit 
+       new((*fDigits)[fNdigits]) AliPHOSDigit( *digit ) ;
+       fNdigits++ ;
     }
-  }
+  }    
     
   fDigits->Compress() ;  
   
   fNdigits = fDigits->GetEntries() ;
   fDigits->Expand(fNdigits) ;
 
+  Int_t i ;
   for (i = 0 ; i < fNdigits ; i++) { 
-    newdigit = (AliPHOSDigit *) fDigits->At(i) ; 
-    newdigit->SetIndexInList(i) ;     
+    AliPHOSDigit * digit = (AliPHOSDigit *) fDigits->At(i) ; 
+    digit->SetIndexInList(i) ;     
   }
 
+
   gAlice->TreeD()->Fill() ;
 
   gAlice->TreeD()->Write(0,TObject::kOverwrite) ;
@@ -359,20 +424,81 @@ void AliPHOSv1::Hit2Digit(Int_t ntracks){
 
 //___________________________________________________________________________
 void AliPHOSv1::MakeBranch(Option_t* opt, char *file)
-{  
+{ 
+
+
+  char *cH ; 
   // Create new branche in the current Root Tree in the digit Tree
   AliDetector::MakeBranch(opt) ;
-  
-  // Create new branches EMC<i> for hits in EMC modules
-  
-  for( Int_t i=0; i<fGeom->GetNModules(); i++ ) GetEMCModule(i).MakeBranch("EMC",i+1,file);
-  
-  // Create new branches CPV<i> for hits in CPV modules for IHEP geometry
-  
-  if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
-    for( Int_t i=0; i<fGeom->GetNCPVModules(); i++ ) GetCPVModule(i).MakeBranch("CPV",i+1,file);
+
+
+  cH = strstr(opt,"S");
+  //Create a branch for SDigits
+  if( cH ){
+    char branchname[20];
+    sprintf(branchname,"%s",GetName());  
+    if(fSDigits)
+      fSDigits->Clear();
+    else
+      fSDigits = new TClonesArray("AliPHOSDigit",1000);
+    fnSdigits = 0 ;
+    cout << " AliPHOSv1::MakeBranch : " << file << endl ; 
+    gAlice->MakeBranchInTree(gAlice->TreeS(),branchname,&fSDigits,fBufferSize,file);  
   }
   
+  cH = strstr(opt,"D");
+  //Create a branch for Digits
+  if( cH ){
+    char branchname[20];
+    sprintf(branchname,"%s",GetName());  
+    if(fSDigits)
+      fDigits->Clear();
+    else
+      fDigits = new TClonesArray("AliPHOSDigit",1000);
+    fNdigits = 0 ;
+    
+    gAlice->MakeBranchInTree(gAlice->TreeD(),branchname,&fSDigits,fBufferSize,file);  
+  }
+
+  cH = strstr(opt,"R");
+  //Create a branch for Reconstruction
+  if( cH ){
+    char branchname[20];
+
+    Int_t splitlevel = 0 ; 
+    
+    fEmcRecPoints->Delete() ; 
+
+    if ( fEmcRecPoints && gAlice->TreeR() ) {
+      sprintf(branchname,"%sEmcRP",GetName()) ;
+      gAlice->TreeR()->Branch(branchname, "TObjArray", &fEmcRecPoints, fBufferSize, splitlevel) ; 
+    }
+
+    fPpsdRecPoints->Delete() ; 
+
+    if ( fPpsdRecPoints && gAlice->TreeR() ) {
+      sprintf(branchname,"%sPpsdRP",GetName()) ;
+      gAlice->TreeR()->Branch(branchname, "TObjArray", &fPpsdRecPoints, fBufferSize, splitlevel) ;
+    }
+
+    fTrackSegments->Delete() ; 
+    
+    if ( fTrackSegments && gAlice->TreeR() ) { 
+      sprintf(branchname,"%sTS",GetName()) ;
+      gAlice->TreeR()->Branch(branchname, &fTrackSegments, fBufferSize) ;
+    }
+    
+    fRecParticles->Delete() ; 
+    
+    if ( fRecParticles && gAlice->TreeR() ) { 
+      sprintf(branchname,"%sRP",GetName()) ;
+      gAlice->TreeR()->Branch(branchname, &fRecParticles, fBufferSize) ;
+    }
+    
+  }
+
+
+
 }
 
 //_____________________________________________________________________________
@@ -385,41 +511,12 @@ void AliPHOSv1::Reconstruction(AliPHOSReconstructioner * Reconstructioner)
   // 5. Write the Tree to File
   
   fReconstructioner = Reconstructioner ;
-  
-  char branchname[10] ;
-  
+    
   // 1.
 
   //  gAlice->MakeTree("R") ; 
-  Int_t splitlevel = 0 ; 
   
-  fEmcRecPoints->Delete() ; 
-
-  if ( fEmcRecPoints && gAlice->TreeR() ) {
-    sprintf(branchname,"%sEmcRP",GetName()) ;
-    gAlice->TreeR()->Branch(branchname, "TObjArray", &fEmcRecPoints, fBufferSize, splitlevel) ; 
-  }
-
-  fPpsdRecPoints->Delete() ; 
-
-  if ( fPpsdRecPoints && gAlice->TreeR() ) {
-    sprintf(branchname,"%sPpsdRP",GetName()) ;
-    gAlice->TreeR()->Branch(branchname, "TObjArray", &fPpsdRecPoints, fBufferSize, splitlevel) ;
-  }
-
-  fTrackSegments->Delete() ; 
-
-  if ( fTrackSegments && gAlice->TreeR() ) { 
-    sprintf(branchname,"%sTS",GetName()) ;
-    gAlice->TreeR()->Branch(branchname, &fTrackSegments, fBufferSize) ;
-  }
-
-  fRecParticles->Delete() ; 
-
-  if ( fRecParticles && gAlice->TreeR() ) { 
-    sprintf(branchname,"%sRP",GetName()) ;
-    gAlice->TreeR()->Branch(branchname, &fRecParticles, fBufferSize) ;
-  }
+  MakeBranch("R") ;
   
   // 3.
 
@@ -462,10 +559,10 @@ void AliPHOSv1::ResetHits()
   // Yuri Kharlov, 28 September 2000
 
   AliDetector::ResetHits();
-  for (Int_t i=0; i<fGeom->GetNModules(); i++) ((AliPHOSCPVModule*)(*fEMCModules)[i]) -> Clear();
-  if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
-    for (Int_t i=0; i<fGeom->GetNCPVModules(); i++) ((AliPHOSCPVModule*)(*fCPVModules)[i]) -> Clear();
-  }
+  //  for (Int_t i=0; i<fGeom->GetNModules(); i++) ((AliPHOSCPVModule*)(*fEMCModules)[i]) -> Clear();
+  //  if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
+  //    for (Int_t i=0; i<fGeom->GetNCPVModules(); i++) ((AliPHOSCPVModule*)(*fCPVModules)[i]) -> Clear();
+  //  }
   
 }  
 //____________________________________________________________________________
@@ -480,6 +577,12 @@ void AliPHOSv1::ResetReconstruction()
   
 }
 
+//____________________________________________________________________________
+//void AliPHOSv1::SDigits2Digits() {
+//  // Adds the noise to the summable digits and keeps digits above a threshold 
+//  // To make a digit
+//}
+
 //____________________________________________________________________________
 void AliPHOSv1::SetTreeAddress()
 { 
@@ -496,14 +599,14 @@ void AliPHOSv1::SetTreeAddress()
   // Set branch address for the Hits Tree for hits in EMC modules
   // Yuri Kharlov, 23 November 2000.
 
-  for( Int_t i=0; i<fGeom->GetNModules(); i++ ) GetEMCModule(i).SetTreeAddress("EMC",i+1);
+  //  for( Int_t i=0; i<fGeom->GetNModules(); i++ ) GetEMCModule(i).SetTreeAddress("EMC",i+1);
 
   // Set branch address for the Hits Tree for hits in CPV modules for IHEP geometry
   // Yuri Kharlov, 28 September 2000.
 
-  if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
-    for( Int_t i=0; i<fGeom->GetNCPVModules(); i++ ) GetCPVModule(i).SetTreeAddress("CPV",i+1);
-  }
+  //  if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
+  //    for( Int_t i=0; i<fGeom->GetNCPVModules(); i++ ) GetCPVModule(i).SetTreeAddress("CPV",i+1);
+  //  }
 
 }
 
@@ -522,13 +625,43 @@ void AliPHOSv1::StepManager(void)
   Int_t          absid    ;      // absolute cell ID number
   Float_t        xyze[4]={0,0,0,0}  ; // position wrt MRS and energy deposited
   TLorentzVector pos      ;      // Lorentz vector of the track current position
+  TLorentzVector pmom     ;      //momentum of the particle initiated hit
+  Float_t        xyd[2]   ;      //local posiiton of the entering
+  Bool_t         entered = kFALSE    ;  
   Int_t          copy     ;
-  Int_t          i        ;
 
   Int_t tracknumber =  gAlice->CurrentTrack() ; 
   Int_t primary     =  gAlice->GetPrimary( gAlice->CurrentTrack() ); 
   TString name      =  fGeom->GetName() ; 
-  Int_t trackpid    =  gMC->TrackPid() ; 
+  Int_t trackpid    =  0  ; 
+
+  if( gMC->IsTrackEntering() ){ // create hit with position and momentum of new particle, 
+                                // but may be without energy deposition
+
+    // Current position of the hit in the local ref. system
+      gMC -> TrackPosition(pos);
+      Float_t xyzm[3], xyzd[3] ;
+      Int_t i;
+      for (i=0; i<3; i++) xyzm[i] = pos[i];
+      gMC -> Gmtod (xyzm, xyzd, 1);    // transform coordinate from master to daughter system
+      xyd[0]  = xyzd[0];
+      xyd[1]  =-xyzd[2];
+      
+      // Current momentum of the hit's track in the local ref. system
+      gMC -> TrackMomentum(pmom);
+      Float_t pm[3], pd[3];
+      for (i=0; i<3; i++) pm[i]   = pmom[i];
+      gMC -> Gmtod (pm, pd, 2);        // transform 3-momentum from master to daughter system
+      pmom[0] = pd[0];
+      pmom[1] =-pd[1];
+      pmom[2] =-pd[2];
+
+      trackpid = gMC->TrackPid();
+      entered = kTRUE ;      // Mark to create hit even withou energy deposition
+
+  }
+
+
   if ( name == "GPS2" || name == "MIXT" ) {            // ======> CPV is a GPS' PPSD
 
     if( gMC->CurrentVolID(copy) == gMC->VolId("GCEL") ) // We are inside a gas cell 
@@ -539,7 +672,7 @@ void AliPHOSv1::StepManager(void)
       xyze[2] = pos[2] ;
       xyze[3] = gMC->Edep() ; 
 
-      if ( xyze[3] != 0 ) { // there is deposited energy 
+      if ( (xyze[3] != 0) || entered ) { // there is deposited energy or new particle entering  PPSD
                gMC->CurrentVolOffID(5, relid[0]) ;  // get the PHOS Module number
        if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(5),"PHO1") == 0 ){
          relid[0] += fGeom->GetNModules() - fGeom->GetNPPSDModules();
@@ -555,7 +688,8 @@ void AliPHOSv1::StepManager(void)
                fGeom->RelToAbsNumbering(relid, absid) ; 
 
        // add current hit to the hit list      
-       AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid);
+         AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid, pmom, xyd);
+
 
       } // there is deposited energy 
     } // We are inside the gas of the CPV  
@@ -566,57 +700,18 @@ void AliPHOSv1::StepManager(void)
     // Yuri Kharlov, 28 September 2000
 
     if( gMC->CurrentVolID(copy) == gMC->VolId("CPVQ") &&
-       gMC->IsTrackEntering() &&
-       gMC->TrackCharge() != 0) {
-
-      // Charged track has just entered to the CPV sensitive plane
+       entered &&
+       gMC->TrackCharge() != 0) {      
       
-      AliPHOSv1 &phos = *(AliPHOSv1*)gAlice->GetModule("PHOS");
+      // Digitize the current CPV hit:
+
+      // 1. find pad response and
       
       Int_t moduleNumber;
       gMC->CurrentVolOffID(3,moduleNumber);
       moduleNumber--;
-      
-      // Current position of the hit in the CPV module ref. system
-
-      gMC -> TrackPosition(pos);
-      Float_t xyzm[3], xyzd[3], xyd[3]={0,0,0};
-      Int_t i;
-      for (i=0; i<3; i++) xyzm[i] = pos[i];
-      gMC -> Gmtod (xyzm, xyzd, 1);    // transform coordinate from master to daughter system
-      xyd[0]  = xyzd[0];
-      xyd[1]  =-xyzd[2];
-      
-      // Current momentum of the hit's track in the CPV module ref. system
-      
-      TLorentzVector  pmom;
-      gMC -> TrackMomentum(pmom);
-      Float_t pm[3], pd[3];
-      for (i=0; i<3; i++) pm[i]   = pmom[i];
-      gMC -> Gmtod (pm, pd, 2);        // transform 3-momentum from master to daughter system
-      pmom[0] = pd[0];
-      pmom[1] =-pd[1];
-      pmom[2] =-pd[2];
-
-      // Current particle type of the hit's track
 
-      Int_t ipart = gMC->TrackPid();
-
-      // Add the current particle in the list of the CPV hits.
-
-      phos.GetCPVModule(moduleNumber).AddHit(fIshunt,primary,pmom,xyd,ipart);
-
-      if (fDebug == 1) {
-       printf("CPV hit added to module #%2d: p = (% .4f, % .4f, % .4f, % .4f) GeV,\n",
-              moduleNumber+1,pmom.Px(),pmom.Py(),pmom.Pz(),pmom.E());
-       printf( "                            xy = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
-               xyd[0],xyd[1],ipart,primary);
-      }
-
-      // Digitize the current CPV hit:
 
-      // 1. find pad response and
-      
       TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0);   // array of digits for current hit
       CPVDigitize(pmom,xyd,moduleNumber,cpvDigits);
       
@@ -664,7 +759,7 @@ void AliPHOSv1::StepManager(void)
        xyze[2] = 0. ;
        xyze[3] = cpvDigit->GetQpad() ;                           // amplitude in a pad
        primary = -1;                                             // No need in primary for CPV
-       AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid);
+       AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid, pmom, xyd);
 
        if (cpvDigit->GetQpad() > 0.02) {
          xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5);
@@ -676,6 +771,7 @@ void AliPHOSv1::StepManager(void)
     }
   } // end of IHEP configuration
   
+
   if(gMC->CurrentVolID(copy) == gMC->VolId("PXTL") ) { //  We are inside a PBWO crystal
     gMC->TrackPosition(pos) ;
     xyze[0] = pos[0] ;
@@ -683,62 +779,11 @@ void AliPHOSv1::StepManager(void)
     xyze[2] = pos[2] ;
     xyze[3] = gMC->Edep() ;
 
-    // Track enters to the crystal from the top edge
-
-    if (gMC->IsTrackEntering()) {
-      Float_t posloc[3];
-      gMC -> Gmtod (xyze, posloc, 1);
-      if (posloc[1] > fGeom->GetCrystalSize(1)/2-0.01) {
-       Int_t row,cel;
-       Float_t xyd[3]={0,0,0};
-       AliPHOSv1 &phos = *(AliPHOSv1*)gAlice->GetModule("PHOS");
-
-       Int_t moduleNumber;
-       gMC->CurrentVolOffID(10,moduleNumber);
-
-       if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(10),"PHO1") == 0 )
-         moduleNumber += fGeom->GetNModules() - fGeom->GetNPPSDModules();
-       moduleNumber--;
-
-       gMC->CurrentVolOffID(4, row) ; // get the row  number inside the module
-       gMC->CurrentVolOffID(3, cel) ; // get the cell number inside the module
-       xyd[0] = -(posloc[2] + (cel-0.5-fGeom->GetNZ()  /2) *
-         (fGeom->GetCrystalSize(2) + 2 * fGeom->GetGapBetweenCrystals()));
-       xyd[1] =   posloc[0] + (row-0.5-fGeom->GetNPhi()/2) *
-         (fGeom->GetCrystalSize(0) + 2 * fGeom->GetGapBetweenCrystals());
-
-       // Current momentum of the hit's track in the CPV module ref. system
-       
-       TLorentzVector  pmom;
-       gMC -> TrackMomentum(pmom);
-       Float_t pm[3], pd[3];
-       for (i=0; i<3; i++) pm[i]   = pmom[i];
-       gMC -> Gmtod (pm, pd, 2);        // transform 3-momentum from master to daughter system
-       pmom[0] = pd[0];
-       pmom[1] =-pd[1];
-       pmom[2] =-pd[2];
-       
-       // Current particle type of the hit's track
-       
-       Int_t ipart = gMC->TrackPid();
-
-       // Add the current particle in the list of the EMC hits.
-
-       phos.GetEMCModule(moduleNumber).AddHit(fIshunt,primary,pmom,xyd,ipart);
-
-       if (fDebug == 1) {
-         printf("EMC hit added to module #%2d: p = (% .4f, % .4f, % .4f, % .4f) GeV,\n",
-                moduleNumber+1,pmom.Px(),pmom.Py(),pmom.Pz(),pmom.E());
-         printf( "                            xy = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
-                 xyd[0],xyd[1],ipart,primary);
-       }
-      }
-    }
-
-    // Track is inside the crystal and deposits some energy
+  
+    if ( (xyze[3] != 0) || entered ) {  // Track is inside the crystal and deposits some energy or just entered 
 
-    if ( xyze[3] != 0 ) {
       gMC->CurrentVolOffID(10, relid[0]) ; // get the PHOS module number ;
+
       if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(10),"PHO1") == 0 )
        relid[0] += fGeom->GetNModules() - fGeom->GetNPPSDModules();      
 
@@ -750,10 +795,13 @@ void AliPHOSv1::StepManager(void)
       fGeom->RelToAbsNumbering(relid, absid) ; 
 
       // add current hit to the hit list
-      AddHit(fIshunt, primary,tracknumber, absid, xyze, trackpid);
+       AddHit(fIshunt, primary,tracknumber, absid, xyze, trackpid,pmom, xyd);
+
 
     } // there is deposited energy
   } // we are inside a PHOS Xtal
+
+
 }
 
 //____________________________________________________________________________
index c51912854968408ad5460bad3df1ad2cddf12923..741ae58a1e157f0db99c5efcb4d8c29efb7af5f5 100644 (file)
@@ -39,23 +39,30 @@ public:
   }
   virtual ~AliPHOSv1(void) ;
 
-  virtual void   AddHit( Int_t shunt, Int_t primary, Int_t track, Int_t id, Float_t *hits, Int_t pid) ; 
-  Int_t          Digitize(Float_t Energy);
-  virtual void   Hit2Digit(Int_t event) ;
-  virtual Int_t  IsVersion(void) const {
-    // Gives the version number 
-    return 1 ; 
-  }
+  virtual void   AddHit( Int_t shunt, Int_t primary, Int_t track, Int_t id, Float_t *hits, Int_t pid, TLorentzVector p, Float_t *pos) ; 
+  Float_t        Calibrate(Int_t amp){ return (amp - fDigitizeA)/fDigitizeB ; }
+  Int_t          Digitize(Float_t Energy){ return (Int_t ) (fDigitizeA + Energy*fDigitizeB); }
+  //  virtual void   Hit2Digit(Int_t event) ;
+  virtual void   Hits2SDigits() ;
   virtual void   MakeBranch(Option_t* opt, char *file=0 ) ;
   void           Reconstruction(AliPHOSReconstructioner * Reconstructioner) ;
   void           ResetClusters(){} ;
   virtual void   ResetHits() ; 
+  virtual void   SDigits2Digits() ;  
+  virtual Int_t  IsVersion(void) const {
+    // Gives the version number 
+    return 1 ; 
+  }
+
   virtual void   ResetReconstruction() ; // Reset reconstructed objects
   void           SetReconstructioner(AliPHOSReconstructioner& Reconstructioner) {
     // sets the reconstructionner object to be used
     fReconstructioner = &Reconstructioner ;
   }  
   void           SetDigitThreshold(Float_t th) { fDigitThreshold = th ; } 
+  void           SetPpsdEnergyThreshold(Float_t enth)      { fPpsdEnergyThreshold = enth ; } 
+  void           SetCpvEnergyThreshold(Float_t enth)       { fCpvEnergyThreshold = enth ; } 
+
   virtual void   SetTreeAddress(); 
   virtual void   StepManager(void) ;                              
   virtual TString Version(void){ 
@@ -71,8 +78,8 @@ public:
 
   // IHEP's CPV specific functions
 
-  AliPHOSCPVModule &GetEMCModule(int n) { return *(AliPHOSCPVModule*)fEMCModules->operator[](n); }
-  AliPHOSCPVModule &GetCPVModule(int n) { return *(AliPHOSCPVModule*)fCPVModules->operator[](n); }
+  //  AliPHOSCPVModule &GetEMCModule(int n) { return *(AliPHOSCPVModule*)fEMCModules->operator[](n); }
+  //  AliPHOSCPVModule &GetCPVModule(int n) { return *(AliPHOSCPVModule*)fCPVModules->operator[](n); }
 
   void       CPVDigitize (TLorentzVector p, Float_t *xy, Int_t moduleNumber, TClonesArray *digits) ;
   Float_t    CPVPadResponseFunction(Float_t qhit, Float_t zg, Float_t xg) ;
@@ -81,11 +88,16 @@ public:
 protected:
 
   Float_t fDigitThreshold ;                       // Threshold for the digit registration 
+  Float_t fPpsdEnergyThreshold;                   //PPSD
+  Float_t fCpvEnergyThreshold;                    //CPV  
   Float_t fPinElectronicNoise  ;                  // Electronic Noise in the PIN
-  AliPHOSReconstructioner  * fReconstructioner ;  // Reconstrutioner of the PHOS event: Clusterization and subtracking procedures
+  Float_t fDigitizeA ;                            //Parameters of the 
+  Float_t fDigitizeB ;                            //digitization 
+  Int_t   fnSdigits ; 
+  AliPHOSReconstructioner  * fReconstructioner ;  // Clusterization and subtracking procedures
   AliPHOSTrackSegmentMaker * fTrackSegmentMaker ; // Reconstructioner of the PHOS track segment: 2 x PPSD + 1 x EMC
-  TClonesArray             * fEMCModules;         // Array of EMC modules
-  TClonesArray             * fCPVModules;         // Array of CPV modules for the IHEP's version of CPV
+  //  TClonesArray             * fEMCModules;         // Array of EMC modules
+  //  TClonesArray             * fCPVModules;         // Array of CPV modules for the IHEP's version of CPV
 
   ClassDef(AliPHOSv1,1)  // Implementation of PHOS manager class for layout EMC+PPSD
 
index 6de2a51fb44afb100557ea46f2ffe8460bb4f001..838143c38be11b73afbf2b8cd904cc86956da596 100644 (file)
@@ -69,13 +69,13 @@ AliPHOSv2::~AliPHOSv2()
 }
 
 //____________________________________________________________________________
-void AliPHOSv2::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id, Float_t * hits, Int_t pid)
+void AliPHOSv2::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id, Float_t * hits, Int_t pid, TLorentzVector p, Float_t * lpos)
 {
   // Add a hit to the hit list.
 
   AliPHOSHit *newHit ;
 
-  newHit = new AliPHOSHit(shunt, primary, tracknumber, Id, hits, pid) ;
+  newHit = new AliPHOSHit(shunt, primary, tracknumber, Id, hits, pid, p, lpos) ;
 
   new((*fHits)[fNhits]) AliPHOSHit(*newHit) ;    
   fNhits++ ;
index 340887600f524f9d62d63112f7675c93de0ba6dc..6e6e01ae260aedc039f0b637d41c8f9d1ad1410f 100644 (file)
@@ -24,7 +24,7 @@ public:
   AliPHOSv2(const char *name, const char *title="") ;
   virtual ~AliPHOSv2(void) ;
 
-  virtual void    AddHit( Int_t shunt, Int_t primary, Int_t track, Int_t id, Float_t *hits, Int_t pid ) ; 
+  virtual void    AddHit( Int_t shunt, Int_t primary, Int_t track, Int_t id, Float_t *hits, Int_t pid, TLorentzVector p, Float_t *pos ) ; 
   virtual Int_t   IsVersion(void) const { 
     // Gives the version number 
     return 2 ; 
index 49877304ede3d2254ef5afbf5b1f7d551228d51a..87a718c13fcf2ca6c5c3ba457d51419911bb52e8 100644 (file)
@@ -92,123 +92,225 @@ AliPHOSv3::AliPHOSv3(AliPHOSReconstructioner * Reconstructioner, const char *nam
   fRecalibrationFactor = 6.2 / fLightYieldMean ;
   fElectronsPerGeV = 2.77e+8 ;
 }
-
 //____________________________________________________________________________
+
 void AliPHOSv3::StepManager(void)
 {
   // Accumulates hits as long as the track stays in a single crystal or PPSD gas Cell
-  // Adds the energy deposited in the PIN diode
+
+//    if (gMC->IsTrackEntering())
+//      cout << "Track enters the volume " << gMC->CurrentVolName() << endl;
+//    if (gMC->IsTrackExiting())
+//      cout << "Track leaves the volume " << gMC->CurrentVolName() << endl;
 
   Int_t          relid[4] ;      // (box, layer, row, column) indices
-  Float_t        xyze[4] ;       // position wrt MRS and energy deposited
-  TLorentzVector pos ;
-  Int_t copy;
-  Float_t        lightyield ;   // Light Yield per GeV
-  Float_t        nElectrons ;   // Number of electrons in the PIN diode
-  TString name = fGeom->GetName() ; 
-  Float_t        global[3] ;
-  Float_t        local[3] ;
-  Float_t        lostenergy ;
+  Int_t          absid    ;      // absolute cell ID number
+  Float_t        xyze[4]={0,0,0,0}  ; // position wrt MRS and energy deposited
+  TLorentzVector pos      ;      // Lorentz vector of the track current position
+  TLorentzVector pmom     ;      //momentum of the particle initiated hit
+  Float_t        xyd[2]   ;      //local posiiton of the entering
+  Bool_t         entered = kFALSE    ;  
+  Int_t          copy     ;
 
   Int_t tracknumber =  gAlice->CurrentTrack() ; 
-  Int_t primary =  gAlice->GetPrimary( gAlice->CurrentTrack() ); 
-  Int_t trackpid    =  gMC->TrackPid() ; 
+  Int_t primary     =  gAlice->GetPrimary( gAlice->CurrentTrack() ); 
+  TString name      =  fGeom->GetName() ; 
+  Int_t trackpid    =  0  ; 
+
+  if( gMC->IsTrackEntering() ){ // create hit with position and momentum of new particle, 
+                                // but may be without energy deposition
+
+    // Current position of the hit in the local ref. system
+      gMC -> TrackPosition(pos);
+      Float_t xyzm[3], xyzd[3] ;
+      Int_t i;
+      for (i=0; i<3; i++) xyzm[i] = pos[i];
+      gMC -> Gmtod (xyzm, xyzd, 1);    // transform coordinate from master to daughter system
+      xyd[0]  = xyzd[0];
+      xyd[1]  =-xyzd[2];
+      
+      // Current momentum of the hit's track in the local ref. system
+      gMC -> TrackMomentum(pmom);
+      Float_t pm[3], pd[3];
+      for (i=0; i<3; i++) pm[i]   = pmom[i];
+      gMC -> Gmtod (pm, pd, 2);        // transform 3-momentum from master to daughter system
+      pmom[0] = pd[0];
+      pmom[1] =-pd[1];
+      pmom[2] =-pd[2];
+
+      trackpid = gMC->TrackPid();
+      entered = kTRUE ;      // Mark to create hit even withou energy deposition
+
+  }
+
+
+  if ( name == "GPS2" || name == "MIXT" ) {            // ======> CPV is a GPS' PPSD
 
-  if ( name == "GPS2" ) { // the CPV is a PPSD
     if( gMC->CurrentVolID(copy) == gMC->VolId("GCEL") ) // We are inside a gas cell 
     {
       gMC->TrackPosition(pos) ;
       xyze[0] = pos[0] ;
       xyze[1] = pos[1] ;
       xyze[2] = pos[2] ;
-      xyze[3] = gMC->Edep() ;
-
+      xyze[3] = gMC->Edep() ; 
 
-      if ( xyze[3] != 0 ) { // there is deposited energy 
+      if ( (xyze[3] != 0) || entered ) { // there is deposited energy or new particle entering  PPSD
                gMC->CurrentVolOffID(5, relid[0]) ;  // get the PHOS Module number
+       if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(5),"PHO1") == 0 ){
+         relid[0] += fGeom->GetNModules() - fGeom->GetNPPSDModules();
+       }
                gMC->CurrentVolOffID(3, relid[1]) ;  // get the Micromegas Module number 
-      // 1-> Geom->GetNumberOfModulesPhi() *  fGeom->GetNumberOfModulesZ() upper                         
-      //  >  fGeom->GetNumberOfModulesPhi()  *  fGeom->GetNumberOfModulesZ() lower
+      // 1-> fGeom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() upper
+      //   > fGeom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() lower
                gMC->CurrentVolOffID(1, relid[2]) ;  // get the row number of the cell
         gMC->CurrentVolID(relid[3]) ;        // get the column number 
 
        // get the absolute Id number
 
-       Int_t absid ; 
-               fGeom->RelToAbsNumbering(relid,absid) ; 
-       
+               fGeom->RelToAbsNumbering(relid, absid) ; 
+
+       // add current hit to the hit list      
+         AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid, pmom, xyd);
 
-       AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid);
 
       } // there is deposited energy 
-     } // We are inside the gas of the CPV  
-   } // GPS2 configuration
+    } // We are inside the gas of the CPV  
+  } // GPS2 configuration
+
+  if ( name == "IHEP" || name == "MIXT" ) {       // ======> CPV is a IHEP's one
+
+    // Yuri Kharlov, 28 September 2000
+
+    if( gMC->CurrentVolID(copy) == gMC->VolId("CPVQ") &&
+       entered &&
+       gMC->TrackCharge() != 0) {      
+      
+      // Digitize the current CPV hit:
+
+      // 1. find pad response and
+      
+      Int_t moduleNumber;
+      gMC->CurrentVolOffID(3,moduleNumber);
+      moduleNumber--;
+
+
+      TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0);   // array of digits for current hit
+      CPVDigitize(pmom,xyd,moduleNumber,cpvDigits);
+      
+      Float_t xmean = 0;
+      Float_t zmean = 0;
+      Float_t qsum  = 0;
+      Int_t   idigit,ndigits;
+
+      // 2. go through the current digit list and sum digits in pads
+
+      ndigits = cpvDigits->GetEntriesFast();
+      for (idigit=0; idigit<ndigits-1; idigit++) {
+       AliPHOSCPVDigit  *cpvDigit1 = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(idigit);
+       Float_t x1 = cpvDigit1->GetXpad() ;
+       Float_t z1 = cpvDigit1->GetYpad() ;
+       for (Int_t jdigit=idigit+1; jdigit<ndigits; jdigit++) {
+         AliPHOSCPVDigit  *cpvDigit2 = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(jdigit);
+         Float_t x2 = cpvDigit2->GetXpad() ;
+         Float_t z2 = cpvDigit2->GetYpad() ;
+         if (x1==x2 && z1==z2) {
+           Float_t qsum = cpvDigit1->GetQpad() + cpvDigit2->GetQpad() ;
+           cpvDigit2->SetQpad(qsum) ;
+           cpvDigits->RemoveAt(idigit) ;
+         }
+       }
+      }
+      cpvDigits->Compress() ;
+
+      // 3. add digits to temporary hit list fTmpHits
+
+      ndigits = cpvDigits->GetEntriesFast();
+      for (idigit=0; idigit<ndigits; idigit++) {
+       AliPHOSCPVDigit  *cpvDigit = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(idigit);
+       relid[0] = moduleNumber + 1 ;                             // CPV (or PHOS) module number
+       relid[1] =-1 ;                                            // means CPV
+       relid[2] = cpvDigit->GetXpad() ;                          // column number of a pad
+       relid[3] = cpvDigit->GetYpad() ;                          // row    number of a pad
+       
+       // get the absolute Id number
+       fGeom->RelToAbsNumbering(relid, absid) ; 
+
+       // add current digit to the temporary hit list
+       xyze[0] = 0. ;
+       xyze[1] = 0. ;
+       xyze[2] = 0. ;
+       xyze[3] = cpvDigit->GetQpad() ;                           // amplitude in a pad
+       primary = -1;                                             // No need in primary for CPV
+       AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid, pmom, xyd);
+
+       if (cpvDigit->GetQpad() > 0.02) {
+         xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5);
+         zmean += cpvDigit->GetQpad() * (cpvDigit->GetYpad() + 0.5);
+         qsum  += cpvDigit->GetQpad();
+       }
+      }
+      delete cpvDigits;
+    }
+  } // end of IHEP configuration
   
-  if(gMC->CurrentVolID(copy) == gMC->VolId("PXTL") )//  We are inside a PBWO crystal 
-     {
-       gMC->TrackPosition(pos) ;
-       xyze[0] = pos[0] ;
-       xyze[1] = pos[1] ;
-       xyze[2] = pos[2] ;
-       lostenergy = gMC->Edep() ; 
-       xyze[3] = gMC->Edep() ;
-
-       global[0] = pos[0] ;
-       global[1] = pos[1] ;
-       global[2] = pos[2] ;
-
-       if ( xyze[3] != 0 ) {
-          gMC->CurrentVolOffID(10, relid[0]) ; // get the PHOS module number ;
-          relid[1] = 0   ;                    // means PW04
-          gMC->CurrentVolOffID(4, relid[2]) ; // get the row number inside the module
-          gMC->CurrentVolOffID(3, relid[3]) ; // get the cell number inside the module
 
-      // get the absolute Id number
+  if(gMC->CurrentVolID(copy) == gMC->VolId("PXTL") ) { //  We are inside a PBWO crystal
+    gMC->TrackPosition(pos) ;
+    xyze[0] = pos[0] ;
+    xyze[1] = pos[1] ;
+    xyze[2] = pos[2] ;
+    xyze[3] = gMC->Edep() ;
 
-          Int_t absid ; 
-          fGeom->RelToAbsNumbering(relid,absid) ; 
-         gMC->Gmtod(global, local, 1) ;
-         
-         // calculating number of electrons in the PIN diode asociated to this hit
-         lightyield = gRandom->Poisson(fLightYieldMean) ;
-         nElectrons = lostenergy * lightyield * fIntrinsicPINEfficiency *
-           exp(-fLightYieldAttenuation * (local[1]+fGeom->GetCrystalSize(1)/2.0 ) ) ;
+  
+    if ( (xyze[3] != 0) || entered ) {  // Track is inside the crystal and deposits some energy or just entered 
 
-         xyze[3] = nElectrons * fRecalibrationFactor ;
-         // add current hit to the hit list
-         AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid);
-
-       } // there is deposited energy
-    } // we are inside a PHOS Xtal
-
-   if(gMC->CurrentVolID(copy) == gMC->VolId("PPIN") ) // We are inside de PIN diode 
-     {
-       gMC->TrackPosition(pos) ;
-       xyze[0] = pos[0] ;
-       xyze[1] = pos[1] ;
-       xyze[2] = pos[2] ;
-       lostenergy = gMC->Edep() ;
-       xyze[3] = gMC->Edep() ;
-
-       if ( xyze[3] != 0 ) {
-          gMC->CurrentVolOffID(11, relid[0]) ; // get the PHOS module number ;
-          relid[1] = 0   ;                    // means PW04 and PIN
-          gMC->CurrentVolOffID(5, relid[2]) ; // get the row number inside the module
-          gMC->CurrentVolOffID(4, relid[3]) ; // get the cell number inside the module
+      gMC->CurrentVolOffID(10, relid[0]) ; // get the PHOS module number ;
 
+      if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(10),"PHO1") == 0 )
+       relid[0] += fGeom->GetNModules() - fGeom->GetNPPSDModules();      
+
+      relid[1] = 0   ;                    // means PBW04
+      gMC->CurrentVolOffID(4, relid[2]) ; // get the row number inside the module
+      gMC->CurrentVolOffID(3, relid[3]) ; // get the cell number inside the module
+      
       // get the absolute Id number
+      fGeom->RelToAbsNumbering(relid, absid) ; 
+
+      // add current hit to the hit list
+       AddHit(fIshunt, primary,tracknumber, absid, xyze, trackpid,pmom, xyd);
 
-          Int_t absid ; 
-          fGeom->RelToAbsNumbering(relid,absid) ;
-         
-         // calculating number of electrons in the PIN diode asociated to this hit
-         nElectrons = lostenergy * fElectronsPerGeV ;
-         xyze[3] = nElectrons * fRecalibrationFactor ;
 
+    } // there is deposited energy
+  } // we are inside a PHOS Xtal
+
+  if(gMC->CurrentVolID(copy) == gMC->VolId("PPIN") ) // We are inside de PIN diode 
+    {
+      gMC->TrackPosition(pos) ;
+      xyze[0] = pos[0] ;
+      xyze[1] = pos[1] ;
+      xyze[2] = pos[2] ;
+      Float_t lostenergy = gMC->Edep() ;
+      xyze[3] = gMC->Edep() ;
+      
+      if ( xyze[3] != 0 ) {
+       gMC->CurrentVolOffID(11, relid[0]) ; // get the PHOS module number ;
+       relid[1] = 0   ;                    // means PW04 and PIN
+       gMC->CurrentVolOffID(5, relid[2]) ; // get the row number inside the module
+       gMC->CurrentVolOffID(4, relid[3]) ; // get the cell number inside the module
+       
+       // get the absolute Id number
+       
+       Int_t absid ; 
+       fGeom->RelToAbsNumbering(relid,absid) ;
+       
+       // calculating number of electrons in the PIN diode asociated to this hit
+         Float_t nElectrons = lostenergy * fElectronsPerGeV ;
+         xyze[3] = nElectrons * fRecalibrationFactor ;
+         
          // add current hit to the hit list
-         AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid);
+         AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid,pmom,xyd);
          //printf("PIN volume is  %d, %d, %d, %d \n",relid[0],relid[1],relid[2],relid[3]);
          //printf("Lost energy in the PIN is %f \n",lostenergy) ;
-       } // there is deposited energy
+      } // there is deposited energy
     } // we are inside a PHOS XtalPHOS PIN diode
 }