make it compile with Federico modifications
authorschutz <schutz@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 30 Nov 2000 11:58:21 +0000 (11:58 +0000)
committerschutz <schutz@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 30 Nov 2000 11:58:21 +0000 (11:58 +0000)
30 files changed:
PHOS/AliPHOS.cxx
PHOS/AliPHOS.h
PHOS/AliPHOSAnalyze.cxx
PHOS/AliPHOSAnalyze.h
PHOS/AliPHOSCPVBaseGeometry.h
PHOS/AliPHOSCPVGeometry.h
PHOS/AliPHOSCPVModule.cxx
PHOS/AliPHOSCPVModule.h
PHOS/AliPHOSClusterizer.h
PHOS/AliPHOSClusterizerv1.cxx
PHOS/AliPHOSClusterizerv1.h
PHOS/AliPHOSCpvRecPoint.h
PHOS/AliPHOSEMCAGeometry.cxx
PHOS/AliPHOSGeometry.cxx
PHOS/AliPHOSGeometry.h
PHOS/AliPHOSHit.cxx
PHOS/AliPHOSHit.h
PHOS/AliPHOSIndexToObject.cxx
PHOS/AliPHOSPIDv1.cxx
PHOS/AliPHOSPPSDGeometry.h
PHOS/AliPHOSRecParticle.cxx
PHOS/AliPHOSRecPoint.h
PHOS/AliPHOSReconstructioner.cxx
PHOS/AliPHOSReconstructioner.h
PHOS/AliPHOSTrackSegmentMakerv1.cxx
PHOS/AliPHOSv0.cxx
PHOS/AliPHOSv0.h
PHOS/AliPHOSv1.cxx
PHOS/AliPHOSv1.h
PHOS/AliPHOSv2.cxx

index 1aef9e8..f42c431 100644 (file)
@@ -26,6 +26,7 @@
 //////////////////////////////////////////////////////////////////////////////
 
 // --- ROOT system ---
+#include "TTree.h"
 
 // --- Standard library ---
 
@@ -43,6 +44,7 @@ AliPHOS::AliPHOS():AliDetector()
 {
   // ctor 
   //We do not create objects, because these pointers will be overwritten durin reading from file.
+  fDigits        = 0 ;
   fEmcRecPoints  = 0 ; 
   fPpsdRecPoints = 0 ;
   fTrackSegments = 0 ;
@@ -53,12 +55,12 @@ AliPHOS::AliPHOS():AliDetector()
 AliPHOS::AliPHOS(const char* name, const char* title): AliDetector(name,title) 
 {
   // ctor
-  
-  fEmcRecPoints  = new AliPHOSRecPoint::RecPointsList(10) ; 
-  fPpsdRecPoints = new AliPHOSRecPoint::RecPointsList(10) ;
-  fTrackSegments = new AliPHOSTrackSegment::TrackSegmentsList("AliPHOSTrackSegment", 10) ;
-  fRecParticles  = new AliPHOSRecParticle::RecParticlesList("AliPHOSRecParticle", 10) ;
-  fDebugLevel    = 0;
+
+  fDigits        = 0 ; 
+  fEmcRecPoints  = 0 ; 
+  fPpsdRecPoints = 0 ;
+  fTrackSegments = 0 ;
+  fRecParticles  = 0 ;
   
 }
 //____________________________________________________________________________
@@ -341,29 +343,71 @@ void AliPHOS::CreateMaterials()
   gMC->Gstpar(idtmed[715], "STRA",2.) ;
 
 }
-
 //____________________________________________________________________________
 void AliPHOS::SetTreeAddress()
 { 
   // TBranch *branch;
   AliDetector::SetTreeAddress();
 
- //  //Branch address for TreeR: RecParticles
-//   TTree *treeR = gAlice->TreeR();
-//   if ( treeR && fRecParticles ) {
-//     branch = treeR->GetBranch("PHOSRP");
-//     if (branch) branch->SetAddress(&fRecParticles) ;
-//   }
-//    //Branch address for TreeR: TrackSegments
-//   if ( treeR && fTrackSegments ) {
-//     branch = treeR->GetBranch("PHOSTS");
-//     if (branch) branch->SetAddress(&fTrackSegments) ;
-//   }
-//   //Branch address for TreeR: EmcRecPoint
-//  if ( treeR && fEmcRecPoints ) {
-//     branch = treeR->GetBranch("PHOSEmcRP");
-//     if (branch) branch->SetAddress(&fEmcRecPoints) ;
-  // }
+  TBranch * branch ;
+
+  if(fDigits)
+    fDigits->Clear();
+  else
+    fDigits = new TClonesArray("AliPHOSDigit",1) ;
+
+  if (gAlice->TreeD()  && fDigits ) {
+    branch = gAlice->TreeD()->GetBranch("PHOS");
+    if (branch) branch->SetAddress(&fDigits) ;
+  } 
+
+   TTree *treeR = gAlice->TreeR();
+   
+   //Branch address for TreeR: EmcRecPoint
+
+  if(fEmcRecPoints)
+    fEmcRecPoints->Delete();
+  else
+    fEmcRecPoints  = new AliPHOSRecPoint::RecPointsList(1) ; 
+
+  if ( treeR && fEmcRecPoints ) {
+    branch = treeR->GetBranch("PHOSEmcRP");
+    if (branch) branch->SetAddress(&fEmcRecPoints) ;
+  }
+
+  //Branch address for TreeR: PPSDRecPoint
+  if(fPpsdRecPoints)
+    fPpsdRecPoints->Delete();
+  else
+    fPpsdRecPoints = new AliPHOSRecPoint::RecPointsList(1) ;
+
+  if ( treeR && fPpsdRecPoints ) {
+    branch = treeR->GetBranch("PHOSPpsdRP");
+    if (branch) branch->SetAddress(&fPpsdRecPoints) ;
+  }
+
+  //Branch address for TreeR: TrackSegments
+  if(fTrackSegments)
+    fTrackSegments->Clear() ;
+  else
+    fTrackSegments = new AliPHOSTrackSegment::TrackSegmentsList("AliPHOSTrackSegment", 1) ;
+  
+  if ( treeR && fTrackSegments ) {
+    branch = treeR->GetBranch("PHOSTS");
+    if (branch) branch->SetAddress(&fTrackSegments) ;
+  }
+  
+  //Branch address for TreeR: RecParticles
+  if(fRecParticles)
+    fRecParticles->Clear() ;
+  else
+    fRecParticles  = new AliPHOSRecParticle::RecParticlesList("AliPHOSRecParticle", 1) ;
+
+  if ( treeR && fRecParticles ) {
+    branch = treeR->GetBranch("PHOSRP");
+    if (branch) branch->SetAddress(&fRecParticles) ;
+  }
+
 }
 
 
index 672293d..02b952e 100644 (file)
@@ -44,6 +44,8 @@ class AliPHOS : public AliDetector {
     return &fEmcRecPoints ;
   }
   virtual  AliPHOSGeometry * GetGeometry() = 0 ;
+  virtual void Hit2Digit(){}  // Will convert hits to digits in versions v1 etc
+
   virtual Int_t   IsVersion(void) const { return -1 ; } 
   virtual  AliPHOSRecPoint::RecPointsList ** PpsdRecPoints() {
     // to be redefined when ppsd is present
@@ -66,7 +68,6 @@ class AliPHOS : public AliDetector {
     assert(0==1) ;
     return *this ; 
   }
-  void SetDebugLevel(Int_t deb) {fDebugLevel = deb;}
  
  protected:
   
@@ -74,8 +75,6 @@ class AliPHOS : public AliDetector {
   AliPHOSRecPoint::RecPointsList * fPpsdRecPoints ;        // The RecPoints (clusters) list in PPSD (veto)
   AliPHOSTrackSegment::TrackSegmentsList * fTrackSegments ;// The TrackSegment list in PHOS
   AliPHOSRecParticle::RecParticlesList * fRecParticles ;   // The reconstructed particles list in PHOS
-  Int_t  fDebugLevel ;                                     // Debug level
-
 
   ClassDef(AliPHOS,2) // Photon Spectrometer Detector (base class)
 
index 52b4927..af91f4e 100644 (file)
@@ -53,6 +53,7 @@
 #include "AliPHOSTrackSegment.h"
 #include "AliPHOSRecParticle.h"
 #include "AliPHOSIndexToObject.h"
+#include "AliPHOSHit.h"
 #include "AliPHOSCPVHit.h"
 #include "AliPHOSCpvRecPoint.h"
 
@@ -126,207 +127,206 @@ AliPHOSAnalyze::~AliPHOSAnalyze()
   if(fTrs)        {delete fTrs      ; fTrs     =0 ;}
 
 }
-
 //____________________________________________________________________________
-void AliPHOSAnalyze::ActivePPSD(Int_t Nevents=1){
-  
-  fhEnergyCorrelations  = new TH2F("hEnergyCorrelations","hEnergyCorrelations",40,  0., 0.15, 30, 0., 3.e-5);
+void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod){
+  //Draws pimary particles and reconstructed 
+  //digits, RecPoints, RecPartices etc 
+  //for event Nevent in the module Nmod.
+
+  TH2F * digitOccupancy  = new TH2F("digitOccupancy","EMC digits", 64,-71.,71.,64,-71.,71.);
+  TH2F * emcOccupancy    = new TH2F("emcOccupancy","EMC RecPoints",64,-71.,71.,64,-71.,71.);
+  TH2F * ppsdUp          = new TH2F("ppsdUp","PPSD Up digits",     128,-71.,71.,128,-71.,71.) ;
+  TH2F * ppsdUpCl        = new TH2F("ppsdUpCl","PPSD Up RecPoints",128,-71.,71.,128,-71.,71.) ;
+  TH2F * ppsdLow         = new TH2F("ppsdLow","PPSD Low digits",     128,-71.,71.,128,-71.,71.) ;
+  TH2F * ppsdLowCl       = new TH2F("ppsdLowCl","PPSD Low RecPoints",128,-71.,71.,128,-71.,71.) ;
+  TH2F * nbar            = new TH2F("nbar","Primary nbar",    64,-71.,71.,64,-71.,71.);
+  TH2F * phot            = new TH2F("phot","Primary Photon",  64,-71.,71.,64,-71.,71.);
+  TH2F * charg           = new TH2F("charg","Primary charged",64,-71.,71.,64,-71.,71.);
+  TH2F * recPhot         = new TH2F("recPhot","RecParticles with primary Photon",64,-71.,71.,64,-71.,71.);
+  TH2F * recNbar         = new TH2F("recNbar","RecParticles with primary Nbar",  64,-71.,71.,64,-71.,71.);
+
   //========== Create the Clusterizer
   fClu = new AliPHOSClusterizerv1() ; 
-  fClu->SetEmcEnergyThreshold(0.01) ; 
-  fClu->SetEmcClusteringThreshold(0.20) ; 
-  fClu->SetPpsdEnergyThreshold    (0.0000002) ; 
-  fClu->SetPpsdClusteringThreshold(0.0000001) ; 
-  fClu->SetLocalMaxCut(0.02) ;
-  fClu->SetCalibrationParameters(0., 0.00000001) ; 
 
-  Int_t ievent;
+  fClu->SetEmcEnergyThreshold(0.05) ;
+  fClu->SetEmcClusteringThreshold(0.20) ;
+  fClu->SetPpsdEnergyThreshold    (0.0000002) ;
+  fClu->SetPpsdClusteringThreshold(0.0000001) ;
+  fClu->SetLocalMaxCut(0.03) ;
+  fClu->SetCalibrationParameters(0., 0.00000001) ;
   
-  for ( ievent=0; ievent<Nevents; ievent++)
-    {  
-      
-      //========== Event Number>         
-      if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. ) 
-       cout <<  "AnalyzeResolutions > " << "Event is " << ievent << endl ;  
-      
-      //=========== Connects the various Tree's for evt
-      gAlice->GetEvent(ievent);
-
-      //=========== Gets the Kine TTree
-      gAlice->TreeK()->GetEvent(0) ;
-      
-      //=========== Get the Digit Tree
-      gAlice->TreeD()->GetEvent(0) ;
-      
-      //========== Creating branches ===================================       
-      AliPHOSRecPoint::RecPointsList ** EmcRecPoints =  fPHOS->EmcRecPoints() ;
-      gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP", EmcRecPoints ) ;
-      
-      AliPHOSRecPoint::RecPointsList ** PpsdRecPoints = fPHOS->PpsdRecPoints() ;
-      gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", PpsdRecPoints ) ;
-      
-      AliPHOSTrackSegment::TrackSegmentsList **  TrackSegmentsList = fPHOS->TrackSegments() ;
-      if( (*TrackSegmentsList) )
-       (*TrackSegmentsList)->Clear() ;
-      gAlice->TreeR()->SetBranchAddress( "PHOSTS", TrackSegmentsList ) ;
-      
-      AliPHOSRecParticle::RecParticlesList ** RecParticleList  = fPHOS->RecParticles() ;
-      if( (*RecParticleList) )
-       (*RecParticleList)->Clear() ;
-      gAlice->TreeR()->SetBranchAddress( "PHOSRP", RecParticleList ) ;
-      
-      
-      //=========== Gets the Reconstraction TTree
-      gAlice->TreeR()->GetEvent(0) ;
-            
-      AliPHOSPpsdRecPoint * RecPoint ;
-      Int_t relid[4] ; 
-      TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
-      while( ( RecPoint = (AliPHOSPpsdRecPoint *)nextRP() ) )
-       {
-         if(!(RecPoint->GetUp()) ) {
-           AliPHOSDigit *digit ;
-           Int_t iDigit ;
-           for(iDigit = 0; iDigit < fPHOS->Digits()->GetEntries(); iDigit++) 
-             {
-               digit = (AliPHOSDigit *) fPHOS->Digits()->At(iDigit) ;
-               fGeom->AbsToRelNumbering(digit->GetId(), relid) ;    
-                if((relid[2]==1)&&(relid[3]==1)&&(relid[0]==RecPoint->GetPHOSMod())){
-                 Float_t ConvertorEnergy = fClu->Calibrate(digit->GetAmp()) ;
-                 fhEnergyCorrelations->Fill(ConvertorEnergy,RecPoint->GetTotalEnergy() );  
-                 break ; 
-               } 
-             }
-           break ;
-         }
-       }
-    }
-  SaveHistograms() ;
-  fhEnergyCorrelations->Draw("BOX") ;
-}
-
+  gAlice->GetEvent(Nevent);
+  
+  TClonesArray * primaryList  = gAlice->Particles();
+  
+  TParticle * primary ;
+  Int_t iPrimary ;
+  for ( iPrimary = 0 ; iPrimary < primaryList->GetEntries() ; iPrimary++)
+    {
+      primary = (TParticle*)primaryList->At(iPrimary) ;
+      Int_t primaryType = primary->GetPdgCode() ;
+      if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212) ) {
+        Int_t moduleNumber ;
+        Double_t primX, primZ ;
+        fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
+        if(moduleNumber==Nmod)
+          charg->Fill(primZ,primX,primary->Energy()) ;
+      }
+      if( primaryType == 22 ) {
+        Int_t moduleNumber ;
+        Double_t primX, primZ ;
+        fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
+        if(moduleNumber==Nmod)
+          phot->Fill(primZ,primX,primary->Energy()) ;
+      }
+      else{
+        if( primaryType == -2112 ) {
+          Int_t moduleNumber ;
+          Double_t primX, primZ ;
+          fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
+          if(moduleNumber==Nmod)
+            nbar->Fill(primZ,primX,primary->Energy()) ;
+        }
+      }
+    }  
 
-//____________________________________________________________________________
-void AliPHOSAnalyze::AnalyzeManyEvents(Int_t Nevents, Int_t module)    
-{
-  // analyzes Nevents events in a single PHOS module  
-  // Events should be reconstructed by Reconstruct()
+  fPHOS->SetTreeAddress() ;
 
-  if ( fRootFile == 0 ) 
-    cout << "AnalyzeManyEvents > " << "Root File not openned" << endl ;  
-  else
+  gAlice->TreeD()->GetEvent(0) ;
+  gAlice->TreeR()->GetEvent(0) ;
+  
+  TObjArray ** emcRecPoints =  fPHOS->EmcRecPoints() ;
+  TObjArray ** ppsdRecPoints = fPHOS->PpsdRecPoints() ;
+  TClonesArray ** recParticleList  = fPHOS->RecParticles() ;
+  
+  Int_t iDigit ;
+  AliPHOSDigit * digit ;
+  
+  for(iDigit = 0; iDigit < fPHOS->Digits()->GetEntries(); iDigit++)
     {
-      //========== Booking Histograms
-      cout << "AnalyzeManyEvents > " << "Booking Histograms" << endl ; 
-      BookingHistograms();
-
-      Int_t ievent;
-      Int_t relid[4] ; 
-      AliPHOSDigit * digit ;
-      AliPHOSEmcRecPoint * emc ;
-      AliPHOSPpsdRecPoint * ppsd ;
-      //      AliPHOSTrackSegment * tracksegment ;
-      AliPHOSRecParticle * recparticle;
-
-      for ( ievent=0; ievent<Nevents; ievent++)
-       {  
-         //========== Event Number>         
-         if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. ) 
-           cout <<  "AnalyzeManyEvents > " << "Event is " << ievent << endl ;  
-
-         //=========== Connects the various Tree's for evt
-         gAlice->GetEvent(ievent);
-
-         //=========== Gets the Digit TTree
-         gAlice->TreeD()->GetEvent(0) ;
-
-         //=========== Gets the number of entries in the Digits array 
-         TIter nextdigit(fPHOS->Digits()) ;
-         while( ( digit = (AliPHOSDigit *)nextdigit() ) )
-           {
-             fGeom->AbsToRelNumbering(digit->GetId(), relid) ;         
-             if (fClu->IsInEmc(digit)) fhEmcDigit->Fill(fClu->Calibrate(digit->GetAmp())) ; 
-             else    
-               {  
-                 if (relid[1]<17) fhVetoDigit->Fill(fClu->Calibrate(digit->GetAmp())); 
-                 if (relid[1]>16) fhConvertorDigit->Fill(fClu->Calibrate(digit->GetAmp()));
-               }
-           }
+      digit = (AliPHOSDigit *) fPHOS->Digits()->At(iDigit) ;
+      Int_t relid[4];
+      fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
+      Float_t x,z ;
+      fGeom->RelPosInModule(relid,x,z) ;
+      Float_t e = fClu->Calibrate(digit->GetAmp()) ;
+      if(relid[0]==Nmod){
+        if(relid[1]==0)  //EMC
+          digitOccupancy->Fill(x,z,e) ;
+        if((relid[1]>0)&&(relid[1]<17))
+          ppsdUp->Fill(x,z,e) ;
+        if(relid[1]>16)
+          ppsdLow->Fill(x,z,e) ;
+      }
+    }
+  
+  Int_t irecp ;
+  TVector3 pos ;
+  
+  for(irecp = 0; irecp < (*emcRecPoints)->GetEntries() ; irecp ++){
+    AliPHOSEmcRecPoint * emc= (AliPHOSEmcRecPoint*)(*emcRecPoints)->At(irecp) ;
+    if(emc->GetPHOSMod()==Nmod){
+      emc->GetLocalPosition(pos) ;
+      emcOccupancy->Fill(pos.X(),pos.Z(),emc->GetEnergy());
+    }
+  }
+  
+  for(irecp = 0; irecp < (*ppsdRecPoints)->GetEntries() ; irecp ++){
+    AliPHOSPpsdRecPoint * ppsd= (AliPHOSPpsdRecPoint *)(*ppsdRecPoints)->At(irecp) ;
+    if(ppsd->GetPHOSMod()==Nmod){
+      ppsd->GetLocalPosition(pos) ;
+      if(ppsd->GetUp())
+        ppsdUpCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy());
+      else
+        ppsdLowCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy());
+    }
+  }
+  
+  AliPHOSRecParticle * recParticle ;
+  Int_t iRecParticle ;
+  for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ )
+    {
+      recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
+      
+      Int_t moduleNumberRec ;
+      Double_t recX, recZ ;
+      fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
+      if(moduleNumberRec == Nmod){
+       
+        Double_t minDistance = 5. ;
+        Int_t closestPrimary = -1 ;
+       
+        Int_t numberofprimaries ;
+        Int_t * listofprimaries  = recParticle->GetPrimaries(numberofprimaries)  ;
+        Int_t index ;
+        TParticle * primary ;
+        Double_t distance = minDistance ;
          
-
-         //=========== Cluster in module
-         TIter nextEmc(*fPHOS->EmcRecPoints()  ) ;
-         while((emc = (AliPHOSEmcRecPoint *)nextEmc())) 
+       for ( index = 0 ; index < numberofprimaries ; index++){
+         primary = (TParticle*)primaryList->At(listofprimaries[index]) ;
+         Int_t moduleNumber ;
+         Double_t primX, primZ ;
+         fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
+         if(moduleNumberRec == moduleNumber)
+           distance = TMath::Sqrt((recX-primX)*(recX-primX)+(recZ-primZ)*(recZ-primZ) ) ;
+         if(minDistance > distance)
            {
-             if ( emc->GetPHOSMod() == module )
-               {  
-                 fhEmcCluster->Fill(  emc->GetTotalEnergy()  ); 
-                 TIter nextPpsd( *fPHOS->PpsdRecPoints()) ;
-                 while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd())) 
-                   {
-                     if ( ppsd->GetPHOSMod() == module )
-                       {                         
-                         if (!ppsd->GetUp()) fhConvertorEmc->Fill(emc->GetTotalEnergy(),ppsd->GetTotalEnergy()) ;
-                       }
-                   } 
-               }
-           }
-
-         //=========== Cluster in module PPSD Down
-         TIter nextPpsd(*fPHOS->PpsdRecPoints() ) ;
-         while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd())) 
-           {
-             if ( ppsd->GetPHOSMod() == module )
-               { 
-                 if (!ppsd->GetUp()) fhConvertorCluster->Fill(ppsd->GetTotalEnergy()) ;
-                 if (ppsd->GetUp())  fhVetoCluster     ->Fill(ppsd->GetTotalEnergy()) ;
-               }
+             minDistance = distance ;
+             closestPrimary = listofprimaries[index] ;
            }
-
-         //========== TRackSegments in the event
-         TIter nextRecParticle(*fPHOS->RecParticles() ) ; 
-         while((recparticle = (AliPHOSRecParticle *)nextRecParticle())) 
-           {
-             if ( recparticle->GetPHOSTrackSegment()->GetPHOSMod() == module )
-               { 
-                 cout << "Particle type is " << recparticle->GetType() << endl ; 
-                 Int_t numberofprimaries = 0 ;
-                 Int_t * listofprimaries = recparticle->GetPrimaries(numberofprimaries) ;
-                 cout << "Number of primaries = " << numberofprimaries << endl ; 
-                 Int_t index ;
-                 for ( index = 0 ; index < numberofprimaries ; index++)
-                   cout << "    primary # " << index << " =  " << listofprimaries[index] << endl ;  
-               }
-           }
-       }   // endfor
-      SaveHistograms();
-    }       // endif
-}           // endfunction
-
+       }
+       
+        if(closestPrimary >=0 ){
+         
+          Int_t primaryType = ((TParticle *)primaryList->At(closestPrimary))->GetPdgCode() ;
+         
+          if(primaryType==22)
+            recPhot->Fill(recZ,recX,recParticle->Energy()) ;
+          else
+            if(primaryType==-2112)
+              recNbar->Fill(recZ,recX,recParticle->Energy()) ; 
+        }
+      }
+    }
+  
+  
+  digitOccupancy->Draw("box") ;
+  emcOccupancy->SetLineColor(2) ;
+  emcOccupancy->Draw("boxsame") ;
+  ppsdUp->SetLineColor(3) ;
+  ppsdUp->Draw("boxsame") ;
+  ppsdLow->SetLineColor(4) ;
+  ppsdLow->Draw("boxsame") ;
+  phot->SetLineColor(8) ;
+  phot->Draw("boxsame") ;
+  nbar->SetLineColor(6) ;
+  nbar->Draw("boxsame") ;
+  
+}
 //____________________________________________________________________________
- void AliPHOSAnalyze::Reconstruct(Int_t Nevents,Int_t FirstEvent )    
+ void AliPHOSAnalyze::Reconstruct(Int_t nevents,Int_t firstEvent )    
 {     
 
-  // Performs reconstruction of EMC and CPV (GPS2 or IHEP)
+  // Performs reconstruction of EMC and CPV (GPS2, IHEP or MIXT)
   // for events from FirstEvent to Nevents
 
   Int_t ievent ;   
-  for ( ievent=FirstEvent; ievent<Nevents; ievent++) {  
-    if (ievent==FirstEvent) {
+  for ( ievent=firstEvent; ievent<nevents; ievent++) {  
+    if (ievent==firstEvent) {
       cout << "Analyze > Starting Reconstructing " << endl ; 
       //========== Create the Clusterizer
       fClu = new AliPHOSClusterizerv1() ; 
       fClu->SetEmcEnergyThreshold(0.05) ; 
       fClu->SetEmcClusteringThreshold(0.20) ; 
       fClu->SetLocalMaxCut(0.03) ;
-      if      (strcmp(fGeom->GetName(),"GPS2") == 0) {
+      if      (strcmp(fGeom->GetName(),"GPS2") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0) {
        fClu->SetPpsdEnergyThreshold    (0.0000002) ; 
        fClu->SetPpsdClusteringThreshold(0.0000001) ; 
       }
-      else if (strcmp(fGeom->GetName(),"IHEP") == 0) {
+      else if (strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0) {
        fClu->SetLocalMaxCutCPV(0.03) ;
        fClu->SetLogWeightCutCPV(4.0) ;
-       fClu->SetPpsdEnergyThreshold    (0.09) ;
+       fClu->SetCpvEnergyThreshold(0.09) ;
       }
       fClu->SetCalibrationParameters(0., 0.00000001) ; 
       
@@ -335,7 +335,7 @@ void AliPHOSAnalyze::AnalyzeManyEvents(Int_t Nevents, Int_t module)
          //      fTrs->UnsetUnfoldFlag() ; 
      
       //========== Creates the particle identifier for GPS2 only
-      if      (strcmp(fGeom->GetName(),"GPS2") == 0) {
+      if      (strcmp(fGeom->GetName(),"GPS2") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0) {
        fPID = new AliPHOSPIDv1() ;
        fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ; 
       }          
@@ -350,15 +350,15 @@ void AliPHOSAnalyze::AnalyzeManyEvents(Int_t Nevents, Int_t module)
       cout <<  "======= Analyze ======> Event " << ievent+1 << endl ;
     
     //=========== Connects the various Tree's for evt
-    gAlice->GetEvent(ievent);
-    
-    //=========== Gets the Digit TTree
-    gAlice->TreeD()->GetEvent(0) ;
-    
+    Int_t tracks = gAlice->GetEvent(ievent);
+
+    fPHOS->Hit2Digit(tracks) ;
+        
     //=========== Do the reconstruction
-    fPHOS->Reconstruction(fRec);
-  }
+    fPHOS->Reconstruction(fRec);    
 
+  }
+  
   if(fClu)      {delete fClu      ; fClu     =0 ;}
   if(fPID)      {delete fPID      ; fPID     =0 ;}
   if(fRec)      {delete fRec      ; fRec     =0 ;}
@@ -487,7 +487,7 @@ void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
       
     Int_t nOfModules = fGeom->GetNModules();
     TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
-    Int_t iModule ; 
+    Int_t iModule = 0;         
     for (iModule=0; iModule < nOfModules; iModule++)
       hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100);
 
@@ -506,7 +506,7 @@ void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
       // Get the Hits Tree for the Primary track itrack
       gAlice->ResetHits();
       gAlice->TreeH()->GetEvent(itrack);
-      for (iModule=0; iModule < nOfModules; iModule++) {
+      for (Int_t iModule=0; iModule < nOfModules; iModule++) {
        cpvModule = fPHOS->GetCPVModule(iModule);
        cpvHits   = cpvModule.Hits();
        nCPVhits  = cpvHits->GetEntriesFast();
@@ -635,8 +635,9 @@ void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
  void AliPHOSAnalyze::InvariantMass(Int_t Nevents )    
 {
   // Calculates Real and Mixed invariant mass distributions
-  const Int_t NMixedEvents = 4 ; //# of events used for calculation of 'mixed' distribution 
-  Int_t MixedLoops = (Int_t )TMath::Ceil(Nevents/NMixedEvents) ;
+
+  Int_t nMixedEvents = 4 ; //# of events used for calculation of 'mixed' distribution 
+  Int_t mixedLoops = (Int_t )TMath::Ceil(Nevents/nMixedEvents) ;
   
   //========== Booking Histograms
   TH2D * hRealEM   = new TH2D("hRealEM",   "Real for EM particles",      250,0.,1.,40,0.,4.) ;
@@ -645,92 +646,88 @@ void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
   TH2D * hMixedPhot= new TH2D("hMixedPhot","Mixed for kPhoton particles",250,0.,1.,40,0.,4.) ;
   
   Int_t ievent;
-  Int_t EventInMixedLoop ;
+  Int_t eventInMixedLoop ;
   
-  Int_t NRecParticles[NMixedEvents] ;
+  Int_t nRecParticles[nMixedEvents] ;
   
-  AliPHOSRecParticle::RecParticlesList * AllRecParticleList  = new TClonesArray("AliPHOSRecParticle", NMixedEvents*1000) ;
+  AliPHOSRecParticle::RecParticlesList * allRecParticleList  = new TClonesArray("AliPHOSRecParticle", nMixedEvents*1000) ;
   
-  for(EventInMixedLoop = 0; EventInMixedLoop < MixedLoops; EventInMixedLoop++  ){
+  for(eventInMixedLoop = 0; eventInMixedLoop < mixedLoops; eventInMixedLoop++  ){
     Int_t iRecPhot = 0 ;
     
-    for ( ievent=0; ievent < NMixedEvents; ievent++){        
+    for ( ievent=0; ievent < nMixedEvents; ievent++){        
       
-      Int_t AbsEventNumber = EventInMixedLoop*NMixedEvents + ievent ;
+      Int_t absEventNumber = eventInMixedLoop*nMixedEvents + ievent ;
       
       //=========== Connects the various Tree's for evt
-      gAlice->GetEvent(AbsEventNumber);
-      
-      //=========== Get the Digit Tree
-      gAlice->TreeD()->GetEvent(0) ;
-      
+      gAlice->GetEvent(absEventNumber);
+
       //========== Creating branches ===================================       
+      fPHOS->SetTreeAddress() ;
       
-      AliPHOSRecParticle::RecParticlesList ** RecParticleList  = fPHOS->RecParticles() ;
-      if( (*RecParticleList) )
-       (*RecParticleList)->Clear() ;
-      gAlice->TreeR()->SetBranchAddress( "PHOSRP", RecParticleList ) ;
-      
-      //=========== Gets the Reconstraction TTree
+      gAlice->TreeD()->GetEvent(0) ;
       gAlice->TreeR()->GetEvent(0) ;
       
-      AliPHOSRecParticle * RecParticle ;
+      TClonesArray ** recParticleList  = fPHOS->RecParticles() ;
+      
+            
+      AliPHOSRecParticle * recParticle ;
       Int_t iRecParticle ;
-      for(iRecParticle = 0; iRecParticle < (*RecParticleList)->GetEntries() ;iRecParticle++ )
+      for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ )
        {
-         RecParticle = (AliPHOSRecParticle *) (*RecParticleList)->At(iRecParticle) ;
-         if((RecParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
-            (RecParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)){ 
-           new( (*AllRecParticleList)[iRecPhot] ) AliPHOSRecParticle(*RecParticle) ;
+         recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
+         if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
+            (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)){ 
+           new( (*allRecParticleList)[iRecPhot] ) AliPHOSRecParticle(*recParticle) ;
            iRecPhot++;
          }
        }
       
-       NRecParticles[ievent] = iRecPhot-1 ;  
+       nRecParticles[ievent] = iRecPhot-1 ;  
     }
     
     //Now calculate invariant mass:
     Int_t irp1,irp2 ;
-    Int_t NCurEvent = 0 ;
+    Int_t nCurEvent = 0 ;
 
-    for(irp1 = 0; irp1 < AllRecParticleList->GetEntries()-1; irp1++){
-      AliPHOSRecParticle * rp1 = (AliPHOSRecParticle *)AllRecParticleList->At(irp1) ;
+    for(irp1 = 0; irp1 < allRecParticleList->GetEntries()-1; irp1++){
+      AliPHOSRecParticle * rp1 = (AliPHOSRecParticle *)allRecParticleList->At(irp1) ;
 
-      for(irp2 = irp1+1; irp2 < AllRecParticleList->GetEntries(); irp2++){
-       AliPHOSRecParticle * rp2 = (AliPHOSRecParticle *)AllRecParticleList->At(irp2) ;
+      for(irp2 = irp1+1; irp2 < allRecParticleList->GetEntries(); irp2++){
+       AliPHOSRecParticle * rp2 = (AliPHOSRecParticle *)allRecParticleList->At(irp2) ;
            
-       Double_t InvMass ;
-       InvMass = (rp1->Energy()+rp2->Energy())*(rp1->Energy()+rp2->Energy())-
+       Double_t invMass ;
+       invMass = (rp1->Energy()+rp2->Energy())*(rp1->Energy()+rp2->Energy())-
          (rp1->Px()+rp2->Px())*(rp1->Px()+rp2->Px())-
          (rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py())-
          (rp1->Pz()+rp2->Pz())*(rp1->Pz()+rp2->Pz()) ;
        
-       if(InvMass> 0)
-         InvMass = TMath::Sqrt(InvMass);
+       if(invMass> 0)
+         invMass = TMath::Sqrt(invMass);
        
-       Double_t Pt ; 
-       Pt = TMath::Sqrt((rp1->Px()+rp2->Px() )*( rp1->Px()+rp2->Px() ) +(rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py()));
+       Double_t pt ; 
+       pt = TMath::Sqrt((rp1->Px()+rp2->Px() )*( rp1->Px()+rp2->Px() ) +(rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py()));
 
-       if(irp1 > NRecParticles[NCurEvent])
-         NCurEvent++;
+       if(irp1 > nRecParticles[nCurEvent])
+         nCurEvent++;
            
-       if(irp2 <= NRecParticles[NCurEvent]){ //'Real' event
-         hRealEM->Fill(InvMass,Pt);
+       if(irp2 <= nRecParticles[nCurEvent]){ //'Real' event
+         hRealEM->Fill(invMass,pt);
          if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
-           hRealPhot->Fill(InvMass,Pt);
+           hRealPhot->Fill(invMass,pt);
        }
        else{
-         hMixedEM->Fill(InvMass,Pt);
+         hMixedEM->Fill(invMass,pt);
          if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
-           hMixedPhot->Fill(InvMass,Pt);
+           hMixedPhot->Fill(invMass,pt);
        } //real-mixed
            
       } //loop over second rp
     }//loop over first rp
-    AllRecParticleList->Delete() ;
+    allRecParticleList->Delete() ;
   } //Loop over events
   
-  delete AllRecParticleList ;
+  delete allRecParticleList ;
   
   //writing output
   TFile output("invmass.root","RECREATE");
@@ -747,6 +744,263 @@ void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
 }
 
 //____________________________________________________________________________
+ void AliPHOSAnalyze::ReadAndPrintEMC(Int_t EvFirst=0, Int_t EvLast=0)    
+{
+  //
+  // 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 ;
+
+    //=========== Connects the various Tree's for evt
+    Int_t ntracks = gAlice->GetEvent(ievent);
+    fPHOS->SetTreeAddress() ;
+    
+    gAlice->TreeD()->GetEvent(0) ;
+    gAlice->TreeR()->GetEvent(0) ;
+
+    // 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);
+    }
+
+    // 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);
+       }
+      }
+    }
+
+//      // 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++) {  
+      
+    //========== 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);
+
+    fPHOS->SetTreeAddress() ;
+    
+    gAlice->TreeD()->GetEvent(0) ;
+    gAlice->TreeR()->GetEvent(0) ;
+
+    // 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();
+      }
+    }
+
+    // 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");
+
+}
+
+//____________________________________________________________________________
  void AliPHOSAnalyze::AnalyzeResolutions(Int_t Nevents )    
 {
   // analyzes Nevents events and calculate Energy and Position resolution as well as
@@ -756,15 +1010,15 @@ void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
   cout << "AnalyzeResolutions > " << "Booking Histograms" << endl ; 
   BookResolutionHistograms();
 
-  Int_t Counter[9][5] ;     
-  Int_t i1,i2,TotalInd = 0 ;
+  Int_t counter[9][5] ;     
+  Int_t i1,i2,totalInd = 0 ;
   for(i1 = 0; i1<9; i1++)
     for(i2 = 0; i2<5; i2++)
-      Counter[i1][i2] = 0 ;
+      counter[i1][i2] = 0 ;
   
-  Int_t TotalPrimary = 0 ;
-  Int_t TotalRecPart = 0 ;
-  Int_t TotalRPwithPrim = 0 ;
+  Int_t totalPrimary = 0 ;
+  Int_t totalRecPart = 0 ;
+  Int_t totalRPwithPrim = 0 ;
   Int_t ievent;
 
   cout << "Start Analysing"<< endl ;
@@ -782,275 +1036,268 @@ void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
       gAlice->TreeK()->GetEvent(0) ;
       
       //=========== Gets the list of Primari Particles
-      TClonesArray * PrimaryList  = gAlice->Particles();     
+      TClonesArray * primaryList  = gAlice->Particles();     
 
-      TParticle * Primary ;
+      TParticle * primary ;
       Int_t iPrimary ;
-      for ( iPrimary = 0 ; iPrimary < PrimaryList->GetEntries() ; iPrimary++)
+      for ( iPrimary = 0 ; iPrimary < primaryList->GetEntries() ; iPrimary++)
        {
-         Primary = (TParticle*)PrimaryList->UncheckedAt(iPrimary) ;
-         Int_t PrimaryType = Primary->GetPdgCode() ;
-         if( PrimaryType == 22 ) {
-           Int_t ModuleNumber ;
-           Double_t PrimX, PrimZ ;
-           fGeom->ImpactOnEmc(Primary->Theta(), Primary->Phi(), ModuleNumber, PrimX, PrimZ) ;
-           if(ModuleNumber){
-             fhPrimary->Fill(Primary->Energy()) ;
-             if(Primary->Energy() > 0.3)
-               TotalPrimary++ ;
+         primary = (TParticle*)primaryList->UncheckedAt(iPrimary) ;
+         Int_t primaryType = primary->GetPdgCode() ;
+         if( primaryType == 22 ) {
+           Int_t moduleNumber ;
+           Double_t primX, primZ ;
+           fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
+           if(moduleNumber){
+             fhPrimary->Fill(primary->Energy()) ;
+             if(primary->Energy() > 0.3)
+               totalPrimary++ ;
            }
          } 
        }
       
-      //=========== Get the Digit Tree
-      gAlice->TreeD()->GetEvent(0) ;
-      
-      //========== Creating branches ===================================       
-      AliPHOSRecPoint::RecPointsList ** EmcRecPoints =  fPHOS->EmcRecPoints() ;
-      gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP", EmcRecPoints ) ;
-      
-      AliPHOSRecPoint::RecPointsList ** PpsdRecPoints = fPHOS->PpsdRecPoints() ;
-      gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", PpsdRecPoints ) ;
-      
-      AliPHOSTrackSegment::TrackSegmentsList **  TrackSegmentsList = fPHOS->TrackSegments() ;
-      if( (*TrackSegmentsList) )
-       (*TrackSegmentsList)->Clear() ;
-      gAlice->TreeR()->SetBranchAddress( "PHOSTS", TrackSegmentsList ) ;
-      
-      AliPHOSRecParticle::RecParticlesList ** RecParticleList  = fPHOS->RecParticles() ;
-      if( (*RecParticleList) )
-       (*RecParticleList)->Clear() ;
-      gAlice->TreeR()->SetBranchAddress( "PHOSRP", RecParticleList ) ;
+      fPHOS->SetTreeAddress() ;
       
-      //=========== Gets the Reconstraction TTree
+      gAlice->TreeD()->GetEvent(0) ;
       gAlice->TreeR()->GetEvent(0) ;
       
-      AliPHOSRecParticle * RecParticle ;
+      TClonesArray ** recParticleList  = fPHOS->RecParticles() ;     
+      
+      AliPHOSRecParticle * recParticle ;
       Int_t iRecParticle ;
-      for(iRecParticle = 0; iRecParticle < (*RecParticleList)->GetEntries() ;iRecParticle++ )
+      for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ )
        {
-         RecParticle = (AliPHOSRecParticle *) (*RecParticleList)->At(iRecParticle) ;
-         fhAllRP->Fill(CorrectEnergy(RecParticle->Energy())) ;
+         recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
+         fhAllRP->Fill(CorrectEnergy(recParticle->Energy())) ;
          
-         Int_t ModuleNumberRec ;
-         Double_t RecX, RecZ ;
-         fGeom->ImpactOnEmc(RecParticle->Theta(), RecParticle->Phi(), ModuleNumberRec, RecX, RecZ) ;
+         Int_t moduleNumberRec ;
+         Double_t recX, recZ ;
+         fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
          
-         Double_t MinDistance = 2. ;
-         Int_t ClosestPrimary = -1 ;
+         Double_t minDistance  = 100. ;
+         Int_t closestPrimary = -1 ;
          
          Int_t numberofprimaries ;
-         Int_t * listofprimaries  = RecParticle->GetPrimaries(numberofprimaries)  ;
+         Int_t * listofprimaries  = recParticle->GetPrimaries(numberofprimaries)  ;
          Int_t index ;
-         TParticle * Primary ;
-         Double_t Distance = MinDistance ;
+         TParticle * primary ;
+         Double_t distance = minDistance ;
+         Double_t dX, dZ; 
+         Double dXmin = 0.; 
+         Double dZmin = 0. ;
          for ( index = 0 ; index < numberofprimaries ; index++){
-           Primary = (TParticle*)PrimaryList->UncheckedAt(listofprimaries[index]) ;
-           Int_t ModuleNumber ;
-           Double_t PrimX, PrimZ ;
-           fGeom->ImpactOnEmc(Primary->Theta(), Primary->Phi(), ModuleNumber, PrimX, PrimZ) ;
-           if(ModuleNumberRec == ModuleNumber)
-             Distance = TMath::Sqrt((RecX-PrimX)*(RecX-PrimX)+(RecZ-PrimZ)*(RecZ-PrimZ) ) ;
-           if(MinDistance > Distance)
-             {
-               MinDistance = Distance ;
-               ClosestPrimary = listofprimaries[index] ;
+           primary = (TParticle*)primaryList->UncheckedAt(listofprimaries[index]) ;
+           Int_t moduleNumber ;
+           Double_t primX, primZ ;
+           fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
+           if(moduleNumberRec == moduleNumber) {
+             dX = recX - primX;
+             dZ = recZ - primZ;
+             distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
+             if(minDistance > distance) {
+               minDistance = distance ;
+               dXmin = dX;
+               dZmin = dZ;
+               closestPrimary = listofprimaries[index] ;
              }
+           }
          }
-         TotalRecPart++ ;
+         totalRecPart++ ;
 
-         if(ClosestPrimary >=0 ){
-           TotalRPwithPrim++;
+         if(closestPrimary >=0 ){
+           totalRPwithPrim++;
            
-           Int_t PrimaryType = ((TParticle *)PrimaryList->At(ClosestPrimary))->GetPdgCode() ;
-//         TParticlePDG* PDGparticle = ((TParticle *)PrimaryList->At(ClosestPrimary))->GetPDG();
+           Int_t primaryType = ((TParticle *)primaryList->At(closestPrimary))->GetPdgCode() ;
+//         TParticlePDG* pDGparticle = ((TParticle *)primaryList->At(closestPrimary))->GetPDG();
 //         Double_t charge =  PDGparticle->Charge() ;
 //         if(charge)
-//           cout <<"Primary " <<PrimaryType << " E " << ((TParticle *)PrimaryList->At(ClosestPrimary))->Energy() << endl ;
-           Int_t PrimaryCode ;
-           switch(PrimaryType)
+//           cout <<"Primary " <<primaryType << " E " << ((TParticle *)primaryList->At(closestPrimary))->Energy() << endl ;
+           Int_t primaryCode ;
+           switch(primaryType)
              {
              case 22:
-               PrimaryCode = 0;  //Photon
-               fhAllEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy()) ;
-               fhAllPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),MinDistance) ;
+               primaryCode = 0;  //Photon
+               fhAllEnergy   ->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), recParticle->Energy()) ;
+               fhAllPosition ->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), minDistance) ;
+               fhAllPositionX->Fill(dXmin);
+               fhAllPositionZ->Fill(dZmin);
                break;
              case 11 :
-               PrimaryCode = 1;  //Electron
+               primaryCode = 1;  //Electron
                break;
              case -11 :
-               PrimaryCode = 1;  //positron
+               primaryCode = 1;  //positron
                break;
              case 321 :
-               PrimaryCode = 4;  //K+
+               primaryCode = 4;  //K+
                break;
              case -321 :
-               PrimaryCode = 4;  //K-
+               primaryCode = 4;  //K-
                break;
              case 310 :
-               PrimaryCode = 4;  //K0s
+               primaryCode = 4;  //K0s
                break;
              case 130 :
-               PrimaryCode = 4;  //K0l
+               primaryCode = 4;  //K0l
                break;
              case 211 :
-               PrimaryCode = 2;  //K0l
+               primaryCode = 2;  //K0l
                break;
              case -211 :
-               PrimaryCode = 2;  //K0l
+               primaryCode = 2;  //K0l
                break;
              case 2212 :
-               PrimaryCode = 2;  //K0l
+               primaryCode = 2;  //K0l
                break;
              case -2212 :
-               PrimaryCode = 2;  //K0l
+               primaryCode = 2;  //K0l
                break;
              default:
-               PrimaryCode = 3; //ELSE
+               primaryCode = 3; //ELSE
                break;
              }
            
-           switch(RecParticle->GetType())
+           switch(recParticle->GetType())
              {
              case AliPHOSFastRecParticle::kGAMMA:
-               if(PrimaryType == 22){
-                 fhPhotEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ; 
-                 fhEMEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ; 
-                 fhPPSDEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ; 
+               if(primaryType == 22){
+                 fhPhotEnergy->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), recParticle->Energy() ) ; 
+                 fhEMEnergy->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), recParticle->Energy() ) ; 
+                 fhPPSDEnergy->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), recParticle->Energy() ) ; 
 
-                 fhPhotPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),MinDistance) ;
-                 fhEMPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),MinDistance) ;
-                 fhPPSDPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),MinDistance) ;
+                 fhPhotPosition->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(),minDistance) ;
+                 fhEMPosition->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(),minDistance) ;
+                 fhPPSDPosition->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(),minDistance) ;
 
-                 fhPhotReg->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-                 fhPhotEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-                 fhPhotPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+                 fhPhotReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+                 fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+                 fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
 
-                 fhPhotPhot->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+                 fhPhotPhot->Fill(CorrectEnergy(recParticle->Energy()) ) ;
                }
-               if(PrimaryType == 2112){ //neutron
-                 fhNReg->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-                 fhNEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-                 fhNPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+               if(primaryType == 2112){ //neutron
+                 fhNReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+                 fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+                 fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
                }
                
-               if(PrimaryType == -2112){ //neutron ~
-                 fhNBarReg->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-                 fhNBarEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-                 fhNBarPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+               if(primaryType == -2112){ //neutron ~
+                 fhNBarReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+                 fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+                 fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
                  
                }
-               if(PrimaryCode == 2){
-                 fhChargedReg->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-                 fhChargedEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-                 fhChargedPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+               if(primaryCode == 2){
+                 fhChargedReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+                 fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+                 fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
                }
                
-               fhAllReg->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-               fhAllEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-               fhAllPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-               fhShape->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-               fhVeto->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-               fhPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-               Counter[0][PrimaryCode]++;
+               fhAllReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+               fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+               fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+               fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+               fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+               fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+               counter[0][primaryCode]++;
                break;
              case  AliPHOSFastRecParticle::kELECTRON:
-               if(PrimaryType == 22){ 
-                 fhPhotElec->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-                 fhEMEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ; 
-                 fhEMPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),MinDistance) ;
-                 fhPhotEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-                 fhPhotPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+               if(primaryType == 22){ 
+                 fhPhotElec->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+                 fhEMEnergy->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), recParticle->Energy() ) ; 
+                 fhEMPosition->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(),minDistance) ;
+                 fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+                 fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
                }         
-               if(PrimaryType == 2112){ //neutron
-                 fhNEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-                 fhNPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+               if(primaryType == 2112){ //neutron
+                 fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+                 fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
                }
                
-               if(PrimaryType == -2112){ //neutron ~
-                 fhNBarEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-                 fhNBarPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+               if(primaryType == -2112){ //neutron ~
+                 fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+                 fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
                  
                }
-               if(PrimaryCode == 2){
-                 fhChargedEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-                 fhChargedPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+               if(primaryCode == 2){
+                 fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+                 fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
                }
                
-               fhAllEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-               fhAllPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-               fhShape->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-               fhPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-               Counter[1][PrimaryCode]++;
+               fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+               fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+               fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+               fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+               counter[1][primaryCode]++;
                break;
              case  AliPHOSFastRecParticle::kNEUTRALHA:
-               if(PrimaryType == 22) 
-                 fhPhotNeuH->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+               if(primaryType == 22) 
+                 fhPhotNeuH->Fill(CorrectEnergy(recParticle->Energy()) ) ;
 
-               fhVeto->Fill(CorrectEnergy(RecParticle->Energy()) ) ;           
-               Counter[2][PrimaryCode]++;
+               fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;           
+               counter[2][primaryCode]++;
                break ;
              case  AliPHOSFastRecParticle::kNEUTRALEM:
-               if(PrimaryType == 22){
-                 fhEMEnergy->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),RecParticle->Energy() ) ; 
-                 fhEMPosition->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),MinDistance ) ;
+               if(primaryType == 22){
+                 fhEMEnergy->Fill(((TParticle *)primaryList->At(closestPrimary))->Energy(),recParticle->Energy() ) ; 
+                 fhEMPosition->Fill(((TParticle *)primaryList->At(closestPrimary))->Energy(),minDistance ) ;
                
-                 fhPhotNuEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-                 fhPhotEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+                 fhPhotNuEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+                 fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
                }
-               if(PrimaryType == 2112) //neutron
-                 fhNEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+               if(primaryType == 2112) //neutron
+                 fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
                
-               if(PrimaryType == -2112) //neutron ~
-                 fhNBarEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+               if(primaryType == -2112) //neutron ~
+                 fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
                
-               if(PrimaryCode == 2)
-                 fhChargedEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+               if(primaryCode == 2)
+                 fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
                
-               fhAllEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-               fhShape->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-               fhVeto->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+               fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+               fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+               fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
 
-               Counter[3][PrimaryCode]++;
+               counter[3][primaryCode]++;
                break ;
              case  AliPHOSFastRecParticle::kCHARGEDHA:
-               if(PrimaryType == 22) //photon
-                 fhPhotChHa->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+               if(primaryType == 22) //photon
+                 fhPhotChHa->Fill(CorrectEnergy(recParticle->Energy()) ) ;
                
-               Counter[4][PrimaryCode]++ ;
+               counter[4][primaryCode]++ ;
                break ;
              case  AliPHOSFastRecParticle::kGAMMAHA:
-                 if(PrimaryType == 22){ //photon
-                   fhPhotGaHa->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-                   fhPPSDEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ; 
-                   fhPPSDPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),MinDistance) ;
-                   fhPhotPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+                 if(primaryType == 22){ //photon
+                   fhPhotGaHa->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+                   fhPPSDEnergy->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), recParticle->Energy() ) ; 
+                   fhPPSDPosition->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(),minDistance) ;
+                   fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
                  }
-                 if(PrimaryType == 2112){ //neutron
-                   fhNPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+                 if(primaryType == 2112){ //neutron
+                   fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
                  }
                
-                 if(PrimaryType == -2112){ //neutron ~
-                   fhNBarPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ; 
+                 if(primaryType == -2112){ //neutron ~
+                   fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ; 
                  }
-                 if(PrimaryCode == 2){
-                   fhChargedPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+                 if(primaryCode == 2){
+                   fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
                  }
                
-                 fhAllPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-                 fhVeto->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-                 fhPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-                 Counter[5][PrimaryCode]++ ;
+                 fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+                 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+                 fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+                 counter[5][primaryCode]++ ;
                  break ;       
              case  AliPHOSFastRecParticle::kABSURDEM:        
-               Counter[6][PrimaryCode]++ ;
-               fhShape->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+               counter[6][primaryCode]++ ;
+               fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
                break;
              case  AliPHOSFastRecParticle::kABSURDHA:
-               Counter[7][PrimaryCode]++ ;
+               counter[7][primaryCode]++ ;
                break;
              default:
-               Counter[8][PrimaryCode]++ ;
+               counter[8][primaryCode]++ ;
                break;
              }
          }
@@ -1058,24 +1305,24 @@ void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
     }   // endfor
   SaveHistograms();
   cout << "Resolutions: Analyzed " << Nevents << " event(s)" << endl ;
-  cout << "Resolutions: Total primary       " << TotalPrimary << endl ;
-  cout << "Resoluitons: Total reconstracted " << TotalRecPart << endl ;
-  cout << "TotalReconstructed with Primarie " << TotalRPwithPrim << endl ;
+  cout << "Resolutions: Total primary       " << totalPrimary << endl ;
+  cout << "Resoluitons: Total reconstracted " << totalRecPart << endl ;
+  cout << "TotalReconstructed with Primarie " << totalRPwithPrim << endl ;
   cout << "                        Primary:   Photon   Electron   Ch. Hadr.  Neutr. Hadr  Kaons" << endl ; 
-  cout << "             Detected as photon       " << Counter[0][0] << "          " << Counter[0][1] << "          " << Counter[0][2] << "          " <<Counter[0][3] << "          " << Counter[0][4] << endl ;
-  cout << "           Detected as electron       " << Counter[1][0] << "          " << Counter[1][1] << "          " << Counter[1][2] << "          " <<Counter[1][3] << "          " << Counter[1][4] << endl ; 
-  cout << "     Detected as neutral hadron       " << Counter[2][0] << "          " << Counter[2][1] << "          " << Counter[2][2] << "          " <<Counter[2][3] << "          " << Counter[2][4] << endl ;
-  cout << "         Detected as neutral EM       " << Counter[3][0] << "          " << Counter[3][1] << "          " << Counter[3][2] << "          " <<Counter[3][3] << "          " << Counter[3][4] << endl ;
-  cout << "     Detected as charged hadron       " << Counter[4][0] << "          " << Counter[4][1] << "          " << Counter[4][2] << "          " <<Counter[4][3] << "          " << Counter[4][4] << endl ;
-  cout << "       Detected as gamma-hadron       " << Counter[5][0] << "          " << Counter[5][1] << "          " << Counter[5][2] << "          " <<Counter[5][3] << "          " << Counter[5][4] << endl ;
-  cout << "          Detected as Absurd EM       " << Counter[6][0] << "          " << Counter[6][1] << "          " << Counter[6][2] << "          " <<Counter[6][3] << "          " << Counter[6][4] << endl ;
-  cout << "      Detected as absurd hadron       " << Counter[7][0] << "          " << Counter[7][1] << "          " << Counter[7][2] << "          " <<Counter[7][3] << "          " << Counter[7][4] << endl ;
-  cout << "          Detected as undefined       " << Counter[8][0] << "          " << Counter[8][1] << "          " << Counter[8][2] << "          " <<Counter[8][3] << "          " << Counter[8][4] << endl ;
+  cout << "             Detected as photon       " << counter[0][0] << "          " << counter[0][1] << "          " << counter[0][2] << "          " <<counter[0][3] << "          " << counter[0][4] << endl ;
+  cout << "           Detected as electron       " << counter[1][0] << "          " << counter[1][1] << "          " << counter[1][2] << "          " <<counter[1][3] << "          " << counter[1][4] << endl ; 
+  cout << "     Detected as neutral hadron       " << counter[2][0] << "          " << counter[2][1] << "          " << counter[2][2] << "          " <<counter[2][3] << "          " << counter[2][4] << endl ;
+  cout << "         Detected as neutral EM       " << counter[3][0] << "          " << counter[3][1] << "          " << counter[3][2] << "          " <<counter[3][3] << "          " << counter[3][4] << endl ;
+  cout << "     Detected as charged hadron       " << counter[4][0] << "          " << counter[4][1] << "          " << counter[4][2] << "          " <<counter[4][3] << "          " << counter[4][4] << endl ;
+  cout << "       Detected as gamma-hadron       " << counter[5][0] << "          " << counter[5][1] << "          " << counter[5][2] << "          " <<counter[5][3] << "          " << counter[5][4] << endl ;
+  cout << "          Detected as Absurd EM       " << counter[6][0] << "          " << counter[6][1] << "          " << counter[6][2] << "          " <<counter[6][3] << "          " << counter[6][4] << endl ;
+  cout << "      Detected as absurd hadron       " << counter[7][0] << "          " << counter[7][1] << "          " << counter[7][2] << "          " <<counter[7][3] << "          " << counter[7][4] << endl ;
+  cout << "          Detected as undefined       " << counter[8][0] << "          " << counter[8][1] << "          " << counter[8][2] << "          " <<counter[8][3] << "          " << counter[8][4] << endl ;
       
       for(i1 = 0; i1<9; i1++)
        for(i2 = 0; i2<5; i2++)
-         TotalInd+=Counter[i1][i2] ;
-      cout << "Indentified particles            " << TotalInd << endl ;
+         totalInd+=counter[i1][i2] ;
+      cout << "Indentified particles            " << totalInd << endl ;
       
 }           // endfunction
 
@@ -1137,6 +1384,9 @@ void  AliPHOSAnalyze::BookResolutionHistograms()
   fhEMPosition   = new TH2F("hEMPosition",   "Position of EM with primary photon",    100, 0., 5., 100, 0., 5.);
   fhPPSDPosition = new TH2F("hPPSDPosition", "Position of PPSD with primary photon",  100, 0., 5., 100, 0., 5.);
 
+  fhAllPositionX = new TH1F("hAllPositionX", "#Delta X of any RP with primary photon",100, -2., 2.);
+  fhAllPositionZ = new TH1F("hAllPositionZ", "#Delta X of any RP with primary photon",100, -2., 2.);
+
 //   if(fhAllReg)
 //     delete fhAllReg ;
 //   if(fhPhotReg)
@@ -1226,500 +1476,8 @@ void  AliPHOSAnalyze::BookResolutionHistograms()
   fhPhotNuEM = new TH1F("hPhotNuEM","hPhotNuEM", 100, 0., 5.);   //Photon registered as Neutral EM
   fhPhotChHa = new TH1F("hPhotChHa","hPhotChHa", 100, 0., 5.);   //Photon registered as Charged Hadron
   fhPhotGaHa = new TH1F("hPhotGaHa","hPhotGaHa", 100, 0., 5.);   //Photon registered as Gamma-Hadron
-
-
 }
-//____________________________________________________________________________
-Bool_t AliPHOSAnalyze::Init(Int_t evt)
-{
-  // Do a few initializations: open the root file
-  //                           get the AliRun object
-  //                           defines the clusterizer, tracksegment maker and particle identifier
-  //                           sets the associated parameters
-
-  Bool_t ok = kTRUE ; 
-  
-   //========== Open galice root file  
-
-  if ( fRootFile == 0 ) {
-    Text_t * name  = new Text_t[80] ; 
-    cout << "AnalyzeOneEvent > Enter file root file name : " ;  
-    cin >> name ; 
-    Bool_t ok = OpenRootFile(name) ; 
-    if ( !ok )
-      cout << " AliPHOSAnalyze > Error opening " << name << endl ; 
-    else { 
-      //========== Get AliRun object from file 
-      
-      gAlice = (AliRun*) fRootFile->Get("gAlice") ;
-      
-      //=========== Get the PHOS object and associated geometry from the file 
-      
-      fPHOS  = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
-      fGeom = fPHOS->GetGeometry();
-      //      fGeom  = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
-
-    } // else !ok
-  } // if fRootFile
-  
-  if ( ok ) {
-    
-    //========== Create the Clusterizer
-
-    fClu =  new AliPHOSClusterizerv1() ; 
-    fClu->SetEmcEnergyThreshold(0.030) ; 
-    fClu->SetEmcClusteringThreshold(0.20) ; 
-    fClu->SetPpsdEnergyThreshold    (0.0000002) ; 
-    fClu->SetPpsdClusteringThreshold(0.0000001) ; 
-    fClu->SetLocalMaxCut(0.03) ;
-    fClu->SetCalibrationParameters(0., 0.00000001) ;  
-    cout <<  "AnalyzeOneEvent > using clusterizer " << fClu->GetName() << endl ; 
-    fClu->PrintParameters() ; 
-    
-    //========== Creates the track segment maker
-    
-    fTrs = new AliPHOSTrackSegmentMakerv1() ;
-    cout <<  "AnalyzeOneEvent > using tack segment maker " << fTrs->GetName() << endl ; 
-    //   fTrs->UnsetUnfoldFlag() ;
-    
-    //========== Creates the particle identifier
-    
-    fPID = new AliPHOSPIDv1() ;
-    cout <<  "AnalyzeOneEvent > using particle identifier " << fPID->GetName() << endl ; 
-    //fPID->SetShowerProfileCuts(Float_t l1m, Float_t l1M, Float_t l2m, Float_t l2M) ; 
-    fPID->SetShowerProfileCuts(0.7, 2.0 , 0.6 , 1.5) ; 
 
-    //========== Creates the Reconstructioner  
-    
-    fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
-//      fRec -> SetDebugReconstruction(kFALSE);     
-    fRec -> SetDebugReconstruction(kTRUE);     
-    
-    //=========== Connect the various Tree's for evt
-    
-    if ( evt == -999 ) {
-      cout <<  "AnalyzeOneEvent > Enter event number : " ; 
-      cin >> evt ;  
-      cout << evt << endl ; 
-    }
-    fEvt = evt ; 
-    gAlice->GetEvent(evt);
-    
-    //=========== Get the Digit TTree
-    
-    gAlice->TreeD()->GetEvent(0) ;     
-    
-  } // ok
-  
-  return ok ; 
-}
-
-
-//____________________________________________________________________________
-void AliPHOSAnalyze::DisplayKineEvent(Int_t evt)
-{
-  // Display particles from the Kine Tree in global Alice (theta, phi) coordinates. 
-  // One PHOS module at the time.
-  // The particle type can be selected.
-  
-  if (evt == -999) 
-    evt = fEvt ;
-
-  Int_t module ; 
-  cout <<  "DisplayKineEvent > which module (1-5,  -1: all) ? " ; 
-  cin >> module ; cout << module << endl ; 
-
-  Int_t testparticle ; 
-  cout << " 22      : PHOTON " << endl 
-       << " (-)11   : (POSITRON)ELECTRON " << endl 
-       << " (-)2112 : (ANTI)NEUTRON " << endl  
-       << " -999    : Everything else " << endl ; 
-  cout  <<  "DisplayKineEvent > enter PDG particle code to display " ; 
-  cin >> testparticle ; cout << testparticle << endl ; 
-
-  Text_t histoname[80] ;
-  sprintf(histoname,"Event %d: Incident particles in module %d", evt, module) ; 
-
-  Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module   
-  fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, AliPHOSGeometry::Degre() ) ;
-
-  Double_t theta, phi ; 
-  fGeom->EmcXtalCoverage(theta, phi, AliPHOSGeometry::Degre() ) ;
-
-  Int_t tdim = (Int_t)( (tM - tm) / theta ) ; 
-  Int_t pdim = (Int_t)( (pM - pm) / phi ) ; 
-
-  tm -= theta ; 
-  tM += theta ; 
-  pm -= phi ; 
-  pM += phi ; 
-
-  TH2F * histoparticle = new TH2F("histoparticle",  histoname, 
-                                         pdim, pm, pM, tdim, tm, tM) ; 
-  histoparticle->SetStats(kFALSE) ;
-
-  // Get pointers to Alice Particle TClonesArray
-
-  TParticle * particle;
-  TClonesArray * particlearray  = gAlice->Particles();    
-
-  Text_t canvasname[80];
-  sprintf(canvasname,"Particles incident in PHOS/EMC module # %d",module) ;
-  TCanvas * kinecanvas = new TCanvas("kinecanvas", canvasname, 650, 500) ; 
-
-  // get the KINE Tree
-
-  TTree * kine =  gAlice->TreeK() ; 
-  Stat_t nParticles =  kine->GetEntries() ; 
-  cout << "DisplayKineEvent > events in kine " << nParticles << endl ; 
-  
-  // loop over particles
-
-  Double_t kRADDEG = 180. / TMath::Pi() ; 
-  Int_t index1 ; 
-  Int_t nparticlein = 0 ; 
-  for (index1 = 0 ; index1 < nParticles ; index1++){
-    Int_t nparticle = particlearray->GetEntriesFast() ;
-    Int_t index2 ; 
-    for( index2 = 0 ; index2 < nparticle ; index2++) {         
-      particle            = (TParticle*)particlearray->UncheckedAt(index2) ;
-      Int_t  particletype = particle->GetPdgCode() ;
-      if (testparticle == -999 || testparticle == particletype) { 
-       Double_t phi        = particle->Phi() ;
-       Double_t theta      = particle->Theta() ;
-       Int_t mod ; 
-       Double_t x, z ; 
-       fGeom->ImpactOnEmc(theta, phi, mod, z, x) ;
-       if ( mod == module ) {
-         nparticlein++ ; 
-         if (particle->Energy() >  fClu->GetEmcClusteringThreshold()  )
-           histoparticle->Fill(phi*kRADDEG, theta*kRADDEG, particle->Energy() ) ; 
-       } 
-      } 
-    }
-  }
-  kinecanvas->Draw() ; 
-  histoparticle->Draw("color") ; 
-  TPaveText *  pavetext = new TPaveText(294, 100, 300, 101); 
-  Text_t text[40] ; 
-  sprintf(text, "Particles: %d ", nparticlein) ;
-  pavetext->AddText(text) ; 
-  pavetext->Draw() ; 
-  kinecanvas->Update(); 
-
-}
-//____________________________________________________________________________
-void AliPHOSAnalyze::DisplayRecParticles()
-{
-  // Display reconstructed particles in global Alice(theta, phi) coordinates. 
-  // One PHOS module at the time.
-  // Click on symbols indicate the reconstructed particle type. 
-
-  if (fEvt == -999) {
-    cout << "DisplayRecParticles > Analyze an event first ... (y/n) " ; 
-    Text_t answer[1] ; 
-    cin >> answer ; cout << answer ; 
-//     if ( answer == "y" ) 
-//       AnalyzeOneEvent() ;
-  } 
-    if (fEvt != -999) {
-      
-      Int_t module ; 
-      cout <<  "DisplayRecParticles > which module (1-5,  -1: all) ? " ; 
-      cin >> module ; cout << module << endl ;
-      Text_t histoname[80] ; 
-      sprintf(histoname,"Event %d: Reconstructed particles in module %d", fEvt, module) ; 
-      Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module   
-      fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, AliPHOSGeometry::Degre() ) ;
-      Double_t theta, phi ; 
-      fGeom->EmcXtalCoverage(theta, phi, AliPHOSGeometry::Degre() ) ;
-      Int_t tdim = (Int_t)( (tM - tm) / theta ) ; 
-      Int_t pdim = (Int_t)( (pM - pm) / phi ) ; 
-      tm -= theta ; 
-      tM += theta ; 
-      pm -= phi ; 
-      TH2F * histoRparticle = new TH2F("histoRparticle",  histoname, 
-                                      pdim, pm, pM, tdim, tm, tM) ; 
-      histoRparticle->SetStats(kFALSE) ;
-      Text_t canvasname[80] ; 
-      sprintf(canvasname, "Reconstructed particles in PHOSmodule # %d", module) ;
-      TCanvas * rparticlecanvas = new TCanvas("RparticleCanvas", canvasname, 650, 500) ; 
-      AliPHOSRecParticle::RecParticlesList * rpl = *fPHOS->RecParticles() ; 
-      Int_t nRecParticles = rpl->GetEntries() ; 
-      Int_t nRecParticlesInModule = 0 ; 
-      TIter nextRecPart(rpl) ; 
-      AliPHOSRecParticle * rp ; 
-      cout << "DisplayRecParticles > " << nRecParticles << " reconstructed particles " << endl ; 
-      Double_t kRADDEG = 180. / TMath::Pi() ; 
-      while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
-       AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ; 
-       if ( ts->GetPHOSMod() == module ) {
-         Int_t numberofprimaries = 0 ;
-         Int_t * listofprimaries = 0;
-         rp->GetPrimaries(numberofprimaries) ;
-         cout << "Number of primaries = " << numberofprimaries << endl ; 
-         Int_t index ;
-         for ( index = 0 ; index < numberofprimaries ; index++)
-           cout << "    primary # " << index << " =  " << listofprimaries[index] << endl ;  
-         
-         nRecParticlesInModule++ ; 
-         Double_t theta = rp->Theta() * kRADDEG ;
-         Double_t phi   = rp->Phi() * kRADDEG ;
-         Double_t energy = rp->Energy() ; 
-         histoRparticle->Fill(phi, theta, energy) ;
-       }
-      }
-      histoRparticle->Draw("color") ; 
-
-      nextRecPart.Reset() ; 
-      while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
-       AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ; 
-       if ( ts->GetPHOSMod() == module )  
-         rp->Draw("P") ; 
-      }
-
-      Text_t text[80] ; 
-      sprintf(text, "reconstructed particles: %d", nRecParticlesInModule) ;
-      TPaveText *  pavetext = new TPaveText(292, 100, 300, 101); 
-      pavetext->AddText(text) ; 
-      pavetext->Draw() ; 
-      rparticlecanvas->Update() ; 
-    }
-}
-
-//____________________________________________________________________________
-void AliPHOSAnalyze::DisplayRecPoints()
-{
-  // Display reconstructed points in local PHOS-module (x, z) coordinates. 
-  // One PHOS module at the time.
-  // Click on symbols displays the EMC cluster, or PPSD information.
-
-  if (fEvt == -999) {
-    cout << "DisplayRecPoints > Analyze an event first ... (y/n) " ; 
-    Text_t answer[1] ; 
-    cin >> answer ; cout << answer ; 
-//     if ( answer == "y" ) 
-//       AnalyzeOneEvent() ;
-  } 
-    if (fEvt != -999) {
-      
-      Int_t module ; 
-      cout <<  "DisplayRecPoints > which module (1-5,  -1: all) ? " ; 
-      cin >> module ; cout << module << endl ; 
-
-      Text_t canvasname[80];
-      sprintf(canvasname,"Digits in PHOS/EMC module # %d",module) ;
-      TCanvas * modulecanvas = new TCanvas("module", canvasname, 650, 500) ; 
-      modulecanvas->Draw() ;
-
-      //=========== Creating 2d-histogram of the PHOS module
-      // a little bit junkie but is used to test Geom functinalities
-
-      Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module   
-      
-      fGeom->EmcModuleCoverage(module, tm, tM, pm, pM); 
-      // convert angles into coordinates local to the EMC module of interest
-
-      Int_t emcModuleNumber ;
-      Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
-      Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
-      fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
-      fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
-      Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;  
-      Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
-      Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ; 
-      Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ; 
-      Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ; 
-      Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;     
-      Text_t histoname[80];
-      sprintf(histoname,"Event %d: Digits and RecPoints in module %d", fEvt, module) ;
-      TH2F * hModule = new TH2F("HistoReconstructed", histoname,
-                               xdim, xmin, xMax, zdim, zmin, zMax) ;  
-      hModule->SetMaximum(2.0);
-      hModule->SetMinimum(0.0);
-      hModule->SetStats(kFALSE); 
-
-      TIter next(fPHOS->Digits()) ;
-      Float_t energy, y, z;
-      Float_t etot=0.;
-      Int_t relid[4]; Int_t nDigits = 0 ;
-      AliPHOSDigit * digit ; 
-
-      // Making 2D histogram of the EMC module
-      while((digit = (AliPHOSDigit *)next())) 
-       {  
-         fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
-         if (relid[0] == module && relid[1] == 0)  
-           {  
-             energy = fClu->Calibrate(digit->GetAmp()) ;
-              cout << "Energy is " << energy << " and threshold is " << fClu->GetEmcEnergyThreshold() << endl; 
-             if (energy >  fClu->GetEmcEnergyThreshold()  ){
-               nDigits++ ;
-               etot += energy ; 
-               fGeom->RelPosInModule(relid,y,z) ;   
-               hModule->Fill(y, z, energy) ;
-             }
-           } 
-       }
-      cout <<"DrawRecPoints >  Found in module " 
-          << module << " " << nDigits << "  digits with total energy " << etot << endl ;
-      hModule->Draw("col2") ;
-
-      //=========== Cluster in module
-
-      //      TClonesArray * emcRP = fPHOS->EmcClusters() ; 
-      TObjArray * emcRP = *(fPHOS->EmcRecPoints()) ; 
-      
-      etot = 0.; 
-      Int_t totalnClusters = 0 ; 
-      Int_t nClusters = 0 ;
-      TIter nextemc(emcRP) ;
-      AliPHOSEmcRecPoint * emc ;
-      while((emc = (AliPHOSEmcRecPoint *)nextemc())) 
-       {
-         //      Int_t numberofprimaries ;
-         //      Int_t * primariesarray = new Int_t[10] ;
-         //      emc->GetPrimaries(numberofprimaries, primariesarray) ;
-         totalnClusters++ ;
-         if ( emc->GetPHOSMod() == module )
-           { 
-             nClusters++ ; 
-             energy = emc->GetTotalEnergy() ;   
-             etot+= energy ;  
-             emc->Draw("M") ;
-           }
-       }
-      cout << "DrawRecPoints > Found " << totalnClusters << " EMC Clusters in PHOS" << endl ; 
-      cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " EMC Clusters " << endl ;
-      cout << "DrawRecPoints > total energy  " << etot << endl ; 
-
-      TPaveText *  pavetext = new TPaveText(22, 80, 83, 90); 
-      Text_t text[40] ; 
-      sprintf(text, "digits: %d;  clusters: %d", nDigits, nClusters) ;
-      pavetext->AddText(text) ; 
-      pavetext->Draw() ; 
-      modulecanvas->Update(); 
-      //=========== Cluster in module PPSD Down
-
-      //      TClonesArray * ppsdRP = fPHOS->PpsdClusters() ;
-      TObjArray * ppsdRP = *(fPHOS->PpsdRecPoints() );
-      etot = 0.; 
-      TIter nextPpsd(ppsdRP) ;
-      AliPHOSPpsdRecPoint * ppsd ;
-      while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd())) 
-       {
-         totalnClusters++ ;
-         if ( ppsd->GetPHOSMod() == module )
-           { 
-             nClusters++ ; 
-             energy = ppsd->GetEnergy() ;   
-             etot+=energy ;  
-             if (!ppsd->GetUp()) ppsd->Draw("P") ;
-           }
-       }
-      cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Down Clusters in PHOS" << endl ; 
-      cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " Ppsd Down Clusters " << endl ;
-      cout << "DrawRecPoints > total energy  " << etot << endl ; 
-
-      //=========== Cluster in module PPSD Up
-  
-      ppsdRP = *(fPHOS->PpsdRecPoints()) ;
-     
-      etot = 0.; 
-      TIter nextPpsdUp(ppsdRP) ;
-      while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsdUp())) 
-       {
-         totalnClusters++ ;
-         if ( ppsd->GetPHOSMod() == module )
-           { 
-             nClusters++ ; 
-             energy = ppsd->GetEnergy() ;   
-             etot+=energy ;  
-             if (ppsd->GetUp()) ppsd->Draw("P") ;
-           }
-       }
-  cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Up Clusters in PHOS" << endl ; 
-  cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " Ppsd Up Clusters " << endl ;
-  cout << "DrawRecPoints > total energy  " << etot << endl ; 
-    
-    } // if !-999
-}
-
-//____________________________________________________________________________
-void AliPHOSAnalyze::DisplayTrackSegments()
-{
-  // Display track segments in local PHOS-module (x, z) coordinates. 
-  // One PHOS module at the time.
-  // One symbol per PHOS subsystem: EMC, upper PPSD, lower PPSD.
-
-  if (fEvt == -999) {
-    cout << "DisplayTrackSegments > Analyze an event first ... (y/n) " ; 
-    Text_t answer[1] ; 
-    cin >> answer ; cout << answer ; 
-//     if ( answer == "y" ) 
-//       AnalyzeOneEvent() ;
-  } 
-    if (fEvt != -999) {
-
-      Int_t module ; 
-      cout <<  "DisplayTrackSegments > which module (1-5,  -1: all) ? " ; 
-      cin >> module ; cout << module << endl ; 
-      //=========== Creating 2d-histogram of the PHOS module
-      // a little bit junkie but is used to test Geom functinalities
-      
-      Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module   
-      
-      fGeom->EmcModuleCoverage(module, tm, tM, pm, pM); 
-      // convert angles into coordinates local to the EMC module of interest
-      
-      Int_t emcModuleNumber ;
-      Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
-      Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
-      fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
-      fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
-      Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;  
-      Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
-      Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ; 
-      Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ; 
-      Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ; 
-      Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;     
-      Text_t histoname[80];
-      sprintf(histoname,"Event %d: Track Segments in module %d", fEvt, module) ; 
-      TH2F * histotrack = new TH2F("histotrack",  histoname, 
-                                  xdim, xmin, xMax, zdim, zmin, zMax) ;  
-      histotrack->SetStats(kFALSE); 
-      Text_t canvasname[80];
-      sprintf(canvasname,"Track segments in PHOS/EMC-PPSD module # %d", module) ;
-      TCanvas * trackcanvas = new TCanvas("TrackSegmentCanvas", canvasname, 650, 500) ; 
-      histotrack->Draw() ; 
-
-      AliPHOSTrackSegment::TrackSegmentsList * trsegl = *(fPHOS->TrackSegments()) ;
-      AliPHOSTrackSegment * trseg ;
-      Int_t nTrackSegments = trsegl->GetEntries() ;
-      Int_t index ;
-      Float_t etot = 0 ;
-      Int_t nTrackSegmentsInModule = 0 ; 
-      for(index = 0; index < nTrackSegments ; index++){
-       trseg = (AliPHOSTrackSegment * )trsegl->At(index) ;
-       etot+= trseg->GetEnergy() ;
-       if ( trseg->GetPHOSMod() == module ) { 
-         nTrackSegmentsInModule++ ; 
-         trseg->Draw("P");
-       }
-      } 
-      Text_t text[80] ; 
-      sprintf(text, "track segments: %d", nTrackSegmentsInModule) ;
-      TPaveText *  pavetext = new TPaveText(22, 80, 83, 90); 
-      pavetext->AddText(text) ; 
-      pavetext->Draw() ; 
-      trackcanvas->Update() ; 
-      cout << "DisplayTrackSegments > Found " << trsegl->GetEntries() << " Track segments with total energy "<< etot << endl ;
-    
-   }
-}
 //____________________________________________________________________________
 Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
 {
@@ -1728,6 +1486,7 @@ Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
   fRootFile   = new TFile(name, "update") ;
   return  fRootFile->IsOpen() ; 
 }
+
 //____________________________________________________________________________
 void AliPHOSAnalyze::SaveHistograms()
 {
@@ -1748,6 +1507,10 @@ void AliPHOSAnalyze::SaveHistograms()
     fhPPSDEnergy->Write() ;
   if(fhAllPosition)
     fhAllPosition->Write() ;
+  if(fhAllPositionX)
+    fhAllPositionX->Write() ;
+  if(fhAllPositionZ)
+    fhAllPositionZ->Write() ;
   if(fhPhotPosition)
     fhPhotPosition->Write() ;
   if(fhEMPosition)
@@ -1838,6 +1601,8 @@ void AliPHOSAnalyze::ResetHistograms()
    fhEMEnergy = 0 ;         // Spectrum of detected electrons with electron primary
    fhPPSDEnergy = 0 ;
    fhAllPosition = 0 ; 
+   fhAllPositionX = 0 ; 
+   fhAllPositionZ = 0 ; 
    fhPhotPosition = 0 ; 
    fhEMPosition = 0 ; 
    fhPPSDPosition = 0 ; 
@@ -1867,5 +1632,4 @@ void AliPHOSAnalyze::ResetHistograms()
    fhPhotChHa = 0 ;
    fhPhotGaHa = 0 ;
 
-
 }
index 07217a8..db247d3 100644 (file)
@@ -36,22 +36,19 @@ public:
   AliPHOSAnalyze(const AliPHOSAnalyze & ana) ; // cpy ctor                   
   virtual ~AliPHOSAnalyze() ;     // dtor
 
-  void ActivePPSD(Int_t Nevents) ;
-  void AnalyzeManyEvents(Int_t Nevtents = 100, Int_t Module=0) ;  // analyzes many events   ;
-  void AnalyzeResolutions(Int_t Nevtents) ; // analyzes Energy and Position resolutions   ;
+  void DrawRecon(Int_t Nevent= 0,Int_t Nmod = 1) ;  // draws positions of entering and 
+                                                    //reconstructed points in PHOS plain
+  void AnalyzeResolutions (Int_t Nevents) ; // analyzes Energy and Position resolutions   ;
+  void ReadAndPrintEMC(Int_t EvFirst=0, Int_t EvLast=0); // Read & print generated and reconstructed hits in EMC
   void ReadAndPrintCPV(Int_t EvFirst=0, Int_t EvLast=0); // Read & print generated and reconstructed hits in CPV
   void AnalyzeCPV(Int_t Nevents);           // analyzes various CPV characteristics
+  void AnalyzeEMC(Int_t Nevents);           // analyzes EMC resolution
   void InvariantMass(Int_t Nevents = 100) ; 
   void Reconstruct(Int_t Nevtents = 100,Int_t FirstEvent = 0) ;
   void BookingHistograms() ;                // booking histograms for the ManyEvent analysis ;
   void BookResolutionHistograms() ;         // booking histograms for the Resoluion analysis ;
   void Copy(TObject & obj) ;                // copies an analysis into an other one   
   Float_t CorrectEnergy(Float_t ERecPart) ;   //Corrects reconstracted energy
-  Bool_t Init(Int_t evt) ;                  // does various initialisations
-  void DisplayKineEvent(Int_t evt = -999) ; // displays the Kine events in ALICE coordinate 
-  void DisplayRecParticles() ;              // displays RecParticles in ALICE coordinate  
-  void DisplayRecPoints() ;                 // displays RecPoints in module coordinate  
-  void DisplayTrackSegments() ;             // displays TrackSegments in module coordinate  
   Bool_t OpenRootFile(Text_t * name) ;      // opens the root file
   void SaveHistograms() ;                   // Save histograms in a root file
   void ResetHistograms() ;                  // 
@@ -88,15 +85,18 @@ public:
   TH1F * fhConvertorCluster ;       // Histo of Cluster energies in Convertor
   TH2F * fhConvertorEmc ;           // 2d Convertor versus Emc energies
 
-  TH2F * fhAllEnergy ;           // Spectrum of detected photons with photon primary
-  TH2F * fhPhotEnergy ;        // Total spectrum of detected photons
+  TH2F * fhAllEnergy ;       // Spectrum of detected photons with photon primary
+  TH2F * fhPhotEnergy ;      // Total spectrum of detected photons
   TH2F * fhEMEnergy ;        // Spectrum of detected neutral EM with EM primary
-  TH2F * fhPPSDEnergy ;
+  TH2F * fhPPSDEnergy ;      // 
 
-  TH2F * fhAllPosition ;        // Position Resolution of  photons with photon primary
-  TH2F * fhPhotPosition ;     // Position Resolution of  photons
-  TH2F * fhEMPosition ;     // Position Resolution of neutral EM with EM primary
-  TH2F * fhPPSDPosition ;  // Position Resolution of neutral EM
+  TH2F * fhAllPosition ;     // Position Resolution of  photons with photon primary
+  TH2F * fhPhotPosition ;    // Position Resolution of  photons
+  TH2F * fhEMPosition ;      // Position Resolution of neutral EM with EM primary
+  TH2F * fhPPSDPosition ;    // Position Resolution of neutral EM
+
+  TH1F * fhAllPositionX ;    // X-Position Resolution of  photons with photon primary
+  TH1F * fhAllPositionZ ;    // Z-Position Resolution of  photons with photon primary
 
   TH1F * fhPhotonPositionY ;        // Y distribution of detected photons
   TH1F * fhElectronPositionY ;      // Y distribution of detected electrons
index 623584f..fa73553 100644 (file)
@@ -48,6 +48,7 @@ public:
   virtual Float_t GetPhiDisplacement(void)          = 0 ;
   virtual Float_t GetPPSDModuleSize(Int_t index)    = 0 ;
   virtual Float_t GetZDisplacement(void)            = 0 ;
+  virtual void    SetLeadConverterThickness(Float_t x) = 0 ;
 
   // Return CPV geometrical parameters
 
index 1d3d7a9..01bdeb8 100644 (file)
@@ -64,6 +64,7 @@ public:
   virtual Int_t   GetNumberOfPadsZ(void)           { AssertPPSD("GetNumberOfPadsZ");           return 0; }
   virtual Int_t   GetNumberOfModulesPhi(void)      { AssertPPSD("GetNumberOfModulesPhi");      return 0; }
   virtual Int_t   GetNumberOfModulesZ(void)        { AssertPPSD("GetNumberOfModulesZ");        return 0; }
+  virtual void    SetLeadConverterThickness(Float_t x)  { AssertPPSD("SetLeadConverterThickness");       }
  
 private:
 
index 4e73945..23df88e 100644 (file)
@@ -15,6 +15,9 @@
 
 /*
   $Log$
+  Revision 1.4  2000/11/20 09:50:03  schutz
+  AliPHOSCPVHit inherits from AliHit
+
   Revision 1.3  2000/11/13 11:45:36  schutz
   DEC and HP compilers sufisfied
 
@@ -144,29 +147,29 @@ void AliPHOSCPVModule::Print(Option_t *opt)
 
 //______________________________________________________________________________
 
-void AliPHOSCPVModule::MakeBranch(Int_t i)
+void AliPHOSCPVModule::MakeBranch(char *title,Int_t i)
 {
   //
-  // Create a new branch for a CPV module #i in the current Root Tree
+  // Create a new branch for a EMC or CPV module #i in the current Root Tree
   //
 
   char branchname[10];
-  sprintf(branchname,"CPV%d",i);
+  sprintf(branchname,"%s%d",title,i);
   gAlice->TreeH()->Branch(branchname,&fHits, 1000);
 }
 
 //_____________________________________________________________________________
-void AliPHOSCPVModule::SetTreeAddress(Int_t i)
+void AliPHOSCPVModule::SetTreeAddress(char *title, Int_t i)
 {
   //
-  // Set branch address for the Hits Tree for a CPV module #i
+  // Set branch address for the Hits Tree for a CPV or EMC module #i
   //
 
   TBranch *branch;
   char branchname[20];
   TTree *treeH = gAlice->TreeH();
   if (treeH){
-    sprintf(branchname,"CPV%d",i);
+    sprintf(branchname,"%s%d",title,i);
     branch = treeH->GetBranch(branchname);
     if (branch) branch->SetAddress(&fHits);
   }
index 82f44f0..a1506d5 100644 (file)
@@ -35,8 +35,8 @@ public:
   void     Clear(Option_t *opt="");
   void     Print(Option_t *opt="");
   void     AddHit(Int_t shunt, Int_t track, TLorentzVector p, Float_t *xy, Int_t ipart);
-  void     MakeBranch(Int_t i);
-  void     SetTreeAddress(Int_t i);
+  void     MakeBranch    (char *title, Int_t i);
+  void     SetTreeAddress(char *title, Int_t i);
   
   TClonesArray *Hits         (void) {return fHits;}
 
index d6c6bdc..0df3eee 100644 (file)
@@ -52,6 +52,8 @@ public:
   virtual void SetLogWeightCutCPV(Float_t w) = 0 ; 
   virtual void SetPpsdClusteringThreshold(Float_t cluth) = 0 ; 
   virtual void SetPpsdEnergyThreshold(Float_t enth) = 0 ; 
+  virtual void SetCpvClusteringThreshold(Float_t cluth) = 0 ; 
+  virtual void SetCpvEnergyThreshold(Float_t enth) = 0 ;  
  
   ClassDef(AliPHOSClusterizer,1)  // Clusterization algorithm class (abstract base class)
 
index 980f2e4..549a07a 100644 (file)
@@ -34,6 +34,8 @@
 //*-- Author: Yves Schutz (SUBATECH) 
 //////////////////////////////////////////////////////////////////////////////
 
+#include <assert.h>
+
 // --- ROOT system ---
 
 #include "TMath.h" 
@@ -68,6 +70,8 @@ AliPHOSClusterizerv1::AliPHOSClusterizerv1()
   fEmcEnergyThreshold      = 0.01;    
   fPpsdClusteringThreshold = 0.0;
   fPpsdEnergyThreshold     = 0.1;  
+  fCpvClusteringThreshold  = 0.0;
+  fCpvEnergyThreshold      = 0.1;  
   fW0                      = 4.5 ;
   fLocMaxCut               = 0.06 ;
   fW0CPV                   = 4.5 ;
@@ -113,8 +117,8 @@ Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)
   }
 
   //Do NOT clusterize upper PPSD  
-  if(fGeom->GetNumberOfCPVLayers() > 1 &&
-     relid1[1] > 0                     &&
+  if( IsInPpsd(d1) && IsInPpsd(d2) &&
+     relid1[1] > 0                 &&
      relid1[1] < fGeom->GetNumberOfPadsPhi()*fGeom->GetNumberOfPadsPhi() ) rv = 2 ;
 
   return rv ; 
@@ -176,8 +180,37 @@ Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit)
   Int_t relid[4] ; 
   fGeom->AbsToRelNumbering(digit->GetId(), relid) ; 
 
-  if ( relid[1] == 0  )
-    rv = kTRUE; 
+  if ( relid[1] == 0  ) rv = kTRUE; 
+
+  return rv ; 
+}
+
+//____________________________________________________________________________
+Bool_t AliPHOSClusterizerv1::IsInPpsd(AliPHOSDigit * digit) 
+{
+  // Tells if (true) or not (false) the digit is in a PHOS-EMC module
+  Bool_t rv = kFALSE ; 
+
+  Int_t relid[4] ; 
+  fGeom->AbsToRelNumbering(digit->GetId(), relid) ; 
+
+  if ( relid[1] > 0 && relid[0] > fGeom->GetNCPVModules() ) rv = kTRUE; 
+
+  return rv ; 
+}
+
+//____________________________________________________________________________
+Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) 
+{
+  // Tells if (true) or not (false) the digit is in a PHOS-EMC module
+  Bool_t rv = kFALSE ; 
+
+  Int_t relid[4] ; 
+  fGeom->AbsToRelNumbering(digit->GetId(), relid) ; 
+
+  if ( relid[1] > 0 && relid[0] <= fGeom->GetNCPVModules() ) rv = kTRUE; 
 
   return rv ; 
 }
@@ -207,8 +240,9 @@ void AliPHOSClusterizerv1::MakeClusters(const DigitsList * dl,
 
     AliPHOSDigit ** clusterdigitslist = new AliPHOSDigit*[dl->GetEntries()] ;   
     Int_t index ;
-    if (( (  IsInEmc(digit) ) && ( Calibrate(digit->GetAmp() ) > fEmcClusteringThreshold  ) ) || 
-        ( ( !IsInEmc(digit) ) && ( Calibrate(digit->GetAmp() ) > fPpsdClusteringThreshold ) ) ) {
+    if (( (  IsInEmc (digit) ) && ( Calibrate(digit->GetAmp() ) > fEmcClusteringThreshold  ) ) || 
+        ( (  IsInPpsd(digit) ) && ( Calibrate(digit->GetAmp() ) > fPpsdClusteringThreshold ) ) ||
+        ( (  IsInCpv (digit) ) && ( Calibrate(digit->GetAmp() ) > fCpvClusteringThreshold  ) ) ) {
   
       Int_t iDigitInCluster = 0 ; 
 
@@ -227,14 +261,18 @@ void AliPHOSClusterizerv1::MakeClusters(const DigitsList * dl,
        
        // start a new PPSD cluster
        if(fNumberOfPpsdClusters >= ppsdl->GetSize()) ppsdl->Expand(2*fNumberOfPpsdClusters+1);
-       if      (strcmp(fGeom->GetName(),"GPS2") == 0) {
+       if      (IsInPpsd(digit)) {
          (*ppsdl)[fNumberOfPpsdClusters] = new AliPHOSPpsdRecPoint() ;
          clu =  (AliPHOSPpsdRecPoint *) ppsdl->At(fNumberOfPpsdClusters)  ;  
        }
-       else if (strcmp(fGeom->GetName(),"IHEP") == 0) {
+       else if (IsInCpv(digit) ) {
          (*ppsdl)[fNumberOfPpsdClusters] = new AliPHOSCpvRecPoint(fW0CPV, fLocMaxCutCPV) ;
          clu =  (AliPHOSCpvRecPoint *) ppsdl->At(fNumberOfPpsdClusters)  ;  
        }
+       else {
+         cout << "AliPHOSClusterizerv1::MakeClusters: unknown configuration " << fGeom->GetName() << endl;
+         assert(0==1);
+       }
        fNumberOfPpsdClusters++ ; 
        clu->AddDigit(*digit, Calibrate(digit->GetAmp()) ) ;    
        clusterdigitslist[iDigitInCluster] = digit  ;   
@@ -307,5 +345,7 @@ void AliPHOSClusterizerv1::PrintParameters()
        << "                       EMC  Clustering threshold = " << fEmcClusteringThreshold << endl
        << "                       EMC  Energy threshold     = " << fEmcEnergyThreshold << endl                  
        << "                      PPSD  Clustering threshold = " << fPpsdClusteringThreshold << endl
-       << "                      PPSD  Energy threshold     = " << fPpsdEnergyThreshold << endl ;                
+       << "                      PPSD  Energy threshold     = " << fPpsdEnergyThreshold << endl
+       << "                       CPV  Clustering threshold = " << fCpvClusteringThreshold << endl
+       << "                       CPV  Energy threshold     = " << fCpvEnergyThreshold << endl ;                
 }
index 7b3e45a..9a02b69 100644 (file)
@@ -43,8 +43,12 @@ public:
   virtual Float_t GetLogWeightCutCPV()         { return fW0CPV;}  
   virtual Float_t GetPpsdClusteringThreshold() { return fPpsdClusteringThreshold;  } 
   virtual Float_t GetPpsdEnergyThreshold()     { return fPpsdEnergyThreshold;  }
+  virtual Float_t GetCpvClusteringThreshold()  { return fCpvClusteringThreshold;  } 
+  virtual Float_t GetCpvEnergyThreshold()      { return fCpvEnergyThreshold;  }
 
-  virtual Bool_t IsInEmc(AliPHOSDigit * digit) ;                      // Tells if id digit is in EMC
+  virtual Bool_t IsInEmc (AliPHOSDigit * digit) ;                      // Tells if id digit is in EMC
+  virtual Bool_t IsInPpsd(AliPHOSDigit * digit) ;                      // Tells if id digit is in PPSD
+  virtual Bool_t IsInCpv (AliPHOSDigit * digit) ;                      // Tells if id digit is in CPV
   virtual void MakeClusters(const DigitsList * dl, 
                            AliPHOSRecPoint::RecPointsList * emcl, 
                            AliPHOSRecPoint::RecPointsList * ppsdl) ; // does the job 
@@ -58,25 +62,29 @@ public:
   virtual void SetLogWeightCutCPV(Float_t w)             { fW0CPV = w ; }
   virtual void SetPpsdClusteringThreshold(Float_t cluth) { fPpsdClusteringThreshold = cluth ; }
   virtual void SetPpsdEnergyThreshold(Float_t enth)      { fPpsdEnergyThreshold = enth ; } 
+  virtual void SetCpvClusteringThreshold(Float_t cluth)  { fCpvClusteringThreshold = cluth ; }
+  virtual void SetCpvEnergyThreshold(Float_t enth)       { fCpvEnergyThreshold = enth ; } 
   
 private:
   
   Float_t fA ;                       // offset of the energy calibration
   Float_t fB ;                       // gain of the energy calibration
   AliPHOSGeometry * fGeom ;          // pointer to geometry
+  Int_t   fNumberOfEmcClusters ;     // number of EMC clusters found 
   Float_t fEmcClusteringThreshold ;  // minimum energy to include a EMC digit in a cluster
   Float_t fEmcEnergyThreshold ;      // minimum energy of EMC digit to be considered
-  Int_t   fNumberOfEmcClusters ;     // number of EMC clusters found 
   Int_t   fNumberOfPpsdClusters ;    // number of PPSD clusters found
   Float_t fPpsdClusteringThreshold ; // minimum energy to include a PPSD digit in a cluster
   Float_t fPpsdEnergyThreshold ;     // minimum energy of PPSD digit to be considered
+  Int_t   fNumberOfCpvClusters ;     // number of CPV clusters found
+  Float_t fCpvClusteringThreshold ;  // minimum energy to include a CPV digit in a cluster
+  Float_t fCpvEnergyThreshold ;      // minimum energy of CPV digit to be considered
   Float_t fLocMaxCut ;               // minimum energy difference to distinguish local maxima in a cluster
   Float_t fW0 ;                      // logarithmic weight for the cluster center of gravity calculation
-
   Float_t fLocMaxCutCPV ;            // minimum energy difference to distinguish local maxima in a CPV cluster
   Float_t fW0CPV ;                   // logarithmic weight for the CPV cluster center of gravity calculation
     
-  ClassDef(AliPHOSClusterizerv1,1)  // Clusterizer implementation version 1
+  ClassDef(AliPHOSClusterizerv1,1)   // Clusterizer implementation version 1
 
 };
 
index efbd3b1..c20c740 100644 (file)
@@ -75,7 +75,7 @@ public:
   void        GetClusterLengths(Int_t &lengX, Int_t &lengZ); // cluster lengths along x and z
   Bool_t      IsCPV(void) {
     // true if the recpoint is in CPV
-    return kTRUE ; } 
+    return (fPHOSMod <= fGeom->GetNCPVModules()) ; } 
   Bool_t      IsSortable() const { 
     // says that emcrecpoints are sortable objects 
     return kTRUE ; } 
index 14abc0b..ad831f2 100644 (file)
@@ -39,12 +39,12 @@ AliPHOSEMCAGeometry::AliPHOSEMCAGeometry()
   fNZ         = 64 ; 
 
   fXtlSize[0] =  2.2 ;
-  fXtlSize[1] = 18.0 ;
+  fXtlSize[1] = 22.0 ;
   fXtlSize[2] =  2.2 ;
 
   // all these numbers coming next are subject to changes
 
-  fOuterBoxThickness[0] = 2.8 ;
+  fOuterBoxThickness[0] = 2.5 ;
   fOuterBoxThickness[1] = 5.0 ;      
   fOuterBoxThickness[2] = 5.0 ;
   
@@ -73,7 +73,7 @@ AliPHOSEMCAGeometry::AliPHOSEMCAGeometry()
   fTextolitBoxThickness[1] = 0.0 ;   
   fTextolitBoxThickness[2] = 3.0 ; 
   
-  fAirThickness[0] =  1.56   ;
+  fAirThickness[0] =  0.4   ;
   fAirThickness[1] = 20.5175 ;  
   fAirThickness[2] =  2.48   ;  
   
index c631f83..cf3c15d 100644 (file)
@@ -54,8 +54,6 @@ AliPHOSGeometry::~AliPHOSGeometry(void)
   if (fRotMatrixArray) fRotMatrixArray->Delete() ; 
   if (fRotMatrixArray) delete fRotMatrixArray ; 
   if (fPHOSAngle     ) delete fPHOSAngle ; 
-//    if (fGeometryEMCA  ) detete fGeometryEMCA;
-//    if (fGeometryCPV   ) detete fGeometryCPV ;
 }
 
 //____________________________________________________________________________
@@ -66,13 +64,32 @@ void AliPHOSGeometry::Init(void)
 
   if ( ((strcmp( fName, "default" )) == 0) || 
        ((strcmp( fName, "GPS2" ))    == 0) ||
-       ((strcmp( fName, "IHEP" ))    == 0) ) {
+       ((strcmp( fName, "IHEP" ))    == 0) ||
+       ((strcmp( fName, "MIXT" ))    == 0) ) {
     fgInit     = kTRUE ; 
-                                             fGeometryEMCA = new AliPHOSEMCAGeometry();
-    if ( ((strcmp( fName, "GPS2" ))  == 0) ) fGeometryCPV  = new AliPHOSPPSDGeometry();
-    if ( ((strcmp( fName, "IHEP" ))  == 0) ) fGeometryCPV  = new AliPHOSCPVGeometry ();
-                                             fGeometrySUPP = new AliPHOSSupportGeometry();
-    fNModules = 5;
+
+    fNModules     = 5;
+    fNPPSDModules = 0;
+    fAngle        = 20;
+
+      fGeometryEMCA = new AliPHOSEMCAGeometry();
+    if      ( ((strcmp( fName, "GPS2" ))  == 0) ) {
+      fGeometryPPSD = new AliPHOSPPSDGeometry();
+      fGeometryCPV  = 0;
+      fNPPSDModules = fNModules;
+    }
+    else if ( ((strcmp( fName, "IHEP" ))  == 0) ) {
+      fGeometryCPV  = new AliPHOSCPVGeometry ();
+      fGeometryPPSD = 0;
+      fNPPSDModules = 0;
+    }
+    else if ( ((strcmp( fName, "MIXT" ))  == 0) ) {
+      fGeometryCPV  = new AliPHOSCPVGeometry ();
+      fGeometryPPSD = new AliPHOSPPSDGeometry();
+      fNPPSDModules = 1;
+    }
+      fGeometrySUPP = new AliPHOSSupportGeometry();
+
     fPHOSAngle = new Float_t[fNModules] ;
     Int_t index ;
     for ( index = 0; index < fNModules; index++ )
@@ -81,10 +98,10 @@ void AliPHOSGeometry::Init(void)
     this->SetPHOSAngles() ; 
     fRotMatrixArray = new TObjArray(fNModules) ; 
   }
- else {
-   fgInit = kFALSE ; 
-   cout << "PHOS Geometry setup: option not defined " << fName << endl ; 
- }
+  else {
+    fgInit = kFALSE ; 
+    cout << "PHOS Geometry setup: option not defined " << fName << endl ; 
+  }
 }
 
 //____________________________________________________________________________
@@ -130,11 +147,13 @@ void AliPHOSGeometry::SetPHOSAngles()
   // Calculates the position in ALICE of the PHOS modules
   
   Double_t const kRADDEG = 180.0 / kPI ;
-  Float_t pphi =  TMath::ATan( GetOuterBoxSize(0)  / ( 2.0 * GetIPtoOuterCoverDistance() ) ) ;
+  Float_t pphi =  2 * TMath::ATan( GetOuterBoxSize(0)  / ( 2.0 * GetIPtoOuterCoverDistance() ) ) ;
   pphi *= kRADDEG ;
+  if (pphi > fAngle) cout << "AliPHOSGeometry: PHOS modules overlap!\n";
+  pphi = fAngle;
   
   for( Int_t i = 1; i <= fNModules ; i++ ) {
-    Float_t angle = pphi * 2 * ( i - fNModules / 2.0 - 0.5 ) ;
+    Float_t angle = pphi * ( i - fNModules / 2.0 - 0.5 ) ;
     fPHOSAngle[i-1] = -  angle ;
   } 
 }
@@ -168,11 +187,34 @@ Bool_t AliPHOSGeometry::AbsToRelNumbering(const Int_t AbsId, Int_t * relid)
     }
     else if ( strcmp(fName,"IHEP") == 0 ) {
       id -=  GetNPhi() * GetNZ() *  GetNModules() ; 
-      relid[0] = (Int_t) TMath::Ceil( id / ( GetNumberOfPadsPhi() * GetNumberOfPadsZ() ) ) ;
+      Float_t nCPV  = GetNumberOfCPVPadsPhi() * GetNumberOfCPVPadsZ() ;
+      relid[0] = (Int_t) TMath::Ceil( id / nCPV ) ;
       relid[1] = 1 ;
-      id -= ( relid[0] - 1 ) * GetNumberOfPadsPhi() * GetNumberOfPadsZ()  ; 
-      relid[2] = (Int_t) TMath::Ceil( id / GetNumberOfPadsZ() ) ;
-      relid[3] = (Int_t) ( id - ( relid[2] - 1 ) * GetNumberOfPadsZ() ) ; 
+      id -= ( relid[0] - 1 ) * nCPV ; 
+      relid[2] = (Int_t) TMath::Ceil( id / GetNumberOfCPVPadsZ() ) ;
+      relid[3] = (Int_t) ( id - ( relid[2] - 1 ) * GetNumberOfCPVPadsZ() ) ; 
+    }
+    else if ( strcmp(fName,"MIXT") == 0 ) {
+      id -=  GetNPhi() * GetNZ() *  GetNModules() ; 
+      Float_t nPPSD = 2 * GetNumberOfModulesPhi() * GetNumberOfModulesZ() *  GetNumberOfPadsPhi() * GetNumberOfPadsZ() ; 
+      Float_t nCPV  = GetNumberOfCPVPadsPhi() * GetNumberOfCPVPadsZ() ;
+      if (id <= nCPV*GetNCPVModules()) { // this pad belons to CPV
+       relid[0] = (Int_t) TMath::Ceil( id / nCPV ) ;
+       relid[1] = 1 ;
+       id -= ( relid[0] - 1 ) * nCPV ; 
+       relid[2] = (Int_t) TMath::Ceil( id / GetNumberOfCPVPadsZ() ) ;
+       relid[3] = (Int_t) ( id - ( relid[2] - 1 ) * GetNumberOfCPVPadsZ() ) ; 
+      }
+      else {                             // this pad belons to PPSD
+       id -= nCPV*GetNCPVModules();
+       relid[0] = (Int_t)TMath::Ceil( id / nPPSD ); 
+       id -= ( relid[0] - 1 ) * nPPSD ;
+       relid[0] += GetNCPVModules();
+       relid[1] = (Int_t)TMath::Ceil( id / ( GetNumberOfPadsPhi() * GetNumberOfPadsZ() ) ) ; 
+       id -= ( relid[1] - 1 ) * GetNumberOfPadsPhi() * GetNumberOfPadsZ() ;
+       relid[2] = (Int_t)TMath::Ceil( id / GetNumberOfPadsPhi() ) ;
+       relid[3] = (Int_t) ( id - ( relid[2] - 1 )  * GetNumberOfPadsPhi() ) ; 
+      }
     }
   } 
   else { // its a PW04 crystal
@@ -352,21 +394,36 @@ void AliPHOSGeometry::ImpactOnEmc(const Double_t theta, const Double_t phi, Int_
 Bool_t AliPHOSGeometry::RelToAbsNumbering(const Int_t * relid, Int_t &  AbsId)
 {
   // Converts the relative numbering into the absolute numbering
-  //  AbsId = 1 to fNModules * fNPhi * fNZ                                                                       -> PbWO4
-  //  AbsId = N(total PHOS crystals) +
-  //          1 to fNModules * 2 * (fNumberOfModulesPhi * fNumberOfModulesZ) * fNumberOfPadsPhi * fNumberOfPadsZ -> PPSD
-  //  AbsId = N(total PHOS crystals) +
-  //          1:fNModules * fNumberOfCPVPadsPhi * fNumberOfCPVPadsZ                                              -> CPV
+  // EMCA crystals:
+  //  AbsId = from 1 to fNModules * fNPhi * fNZ
+  // PPSD gas cell:
+  //  AbsId = from N(total EMCA crystals) + 1
+  //          to NCPVModules * fNumberOfCPVPadsPhi * fNumberOfCPVPadsZ +
+  //          fNModules * 2 * (fNumberOfModulesPhi * fNumberOfModulesZ) * fNumberOfPadsPhi * fNumberOfPadsZ
+  // CPV pad:
+  //  AbsId = from N(total PHOS crystals) + 1
+  //          to NCPVModules * fNumberOfCPVPadsPhi * fNumberOfCPVPadsZ
 
   Bool_t rv = kTRUE ; 
  
-  if      ( relid[1]  >  0 ) { // it is a PPSD pad
-    AbsId =    GetNPhi() * GetNZ() *  GetNModules()                          // the offset to separate EMCA crystals from PPSD pads
-      + ( relid[0] - 1 ) * GetNumberOfModulesPhi() * GetNumberOfModulesZ()   // the pads offset of PHOS modules 
+  if      ( relid[1] > 0 && strcmp(fName,"GPS2")==0) { // it is a PPSD pad
+    AbsId =    GetNPhi() * GetNZ() * GetNModules()                         // the offset to separate EMCA crystals from PPSD pads
+      + ( relid[0] - 1 ) * GetNumberOfModulesPhi() * GetNumberOfModulesZ() // the pads offset of PPSD modules 
                          * GetNumberOfPadsPhi() * GetNumberOfPadsZ() * 2
-      + ( relid[1] - 1 ) * GetNumberOfPadsPhi() * GetNumberOfPadsZ()         // the pads offset of PPSD modules 
-      + ( relid[2] - 1 ) * GetNumberOfPadsPhi()                              // the pads offset of a PPSD row
-      +   relid[3] ;                                                         // the column number
+      + ( relid[1] - 1 ) * GetNumberOfPadsPhi() * GetNumberOfPadsZ()       // the pads offset of PPSD modules 
+      + ( relid[2] - 1 ) * GetNumberOfPadsPhi()                            // the pads offset of a PPSD row
+      +   relid[3] ;                                                       // the column number
+  } 
+
+  else if ( relid[1] > 0 && strcmp(fName,"MIXT")==0) { // it is a PPSD pad
+    AbsId =    GetNPhi() * GetNZ() * GetNModules()                         // the offset to separate EMCA crystals from PPSD pads
+      + GetNCPVModules() * GetNumberOfCPVPadsPhi() * GetNumberOfCPVPadsZ() // the pads offset of CPV modules if any
+      + ( relid[0] - 1 - GetNCPVModules())
+                         * GetNumberOfModulesPhi() * GetNumberOfModulesZ() // the pads offset of PPSD modules 
+                         * GetNumberOfPadsPhi() * GetNumberOfPadsZ() * 2
+      + ( relid[1] - 1 ) * GetNumberOfPadsPhi() * GetNumberOfPadsZ()       // the pads offset of PPSD modules 
+      + ( relid[2] - 1 ) * GetNumberOfPadsPhi()                            // the pads offset of a PPSD row
+      +   relid[3] ;                                                       // the column number
   } 
 
   else if ( relid[1] ==  0 ) { // it is a Phos crystal
@@ -378,8 +435,8 @@ Bool_t AliPHOSGeometry::RelToAbsNumbering(const Int_t * relid, Int_t &  AbsId)
 
   else if ( relid[1] == -1 ) { // it is a CPV pad
     AbsId =    GetNPhi() * GetNZ() *  GetNModules()                          // the offset to separate EMCA crystals from CPV pads
-      + ( relid[0] - 1 ) * GetNumberOfPadsPhi() * GetNumberOfPadsZ()         // the pads offset of PHOS modules 
-      + ( relid[2] - 1 ) * GetNumberOfPadsZ()                                // the pads offset of a CPV row
+      + ( relid[0] - 1 ) * GetNumberOfCPVPadsPhi() * GetNumberOfCPVPadsZ()         // the pads offset of PHOS modules 
+      + ( relid[2] - 1 ) * GetNumberOfCPVPadsZ()                                // the pads offset of a CPV row
       +   relid[3] ;                                                         // the column number
   }
   
@@ -391,56 +448,56 @@ Bool_t AliPHOSGeometry::RelToAbsNumbering(const Int_t * relid, Int_t &  AbsId)
 void AliPHOSGeometry::RelPosInAlice(const Int_t id, TVector3 & pos ) 
 {
   // Converts the absolute numbering into the global ALICE coordinates
+  // It works only for the GPS2 geometry
   
-   if (id > 0) { 
-
-  Int_t relid[4] ;
-  AbsToRelNumbering(id , relid) ;
-
-  Int_t phosmodule = relid[0] ; 
-
-  Float_t y0 = 0 ; 
-
-  if ( relid[1] == 0 ) // it is a PbW04 crystal 
-  {  y0 =  -(GetIPtoOuterCoverDistance() + GetUpperPlateThickness()
-      + GetSecondUpperPlateThickness() + GetUpperCoolingPlateThickness())  ;  
-  }
-  if ( relid[1] > 0 ) { // its a PPSD pad
-    if ( relid[1] >  GetNumberOfModulesPhi() *  GetNumberOfModulesZ() ) // its an bottom module
-     {
-       y0 = -( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() / 2.0)  ;
-     } 
-    else // its an upper module
-      y0 = -( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() - GetLeadToMicro2Gap()
-       -  GetLeadConverterThickness() -  GetMicro1ToLeadGap() - GetMicromegas1Thickness() / 2.0) ; 
-  }
-
-  Float_t x, z ; 
-  RelPosInModule(relid, x, z) ; 
-
-  pos.SetX(x) ;
-  pos.SetZ(z) ;
-  pos.SetY( TMath::Sqrt(x*x + z*z + y0*y0) ) ; 
-
-
-
-   Float_t phi           = GetPHOSAngle( phosmodule) ; 
-   Double_t const kRADDEG = 180.0 / kPI ;
-   Float_t rphi          = phi / kRADDEG ; 
-
-   TRotation rot ;
-   rot.RotateZ(-rphi) ; // a rotation around Z by angle  
-  
-   TRotation dummy = rot.Invert() ;  // to transform from original frame to rotate frame
-  
-   pos.Transform(rot) ; // rotate the baby 
+  if (id > 0 && strcmp(fName,"GPS2")==0) { 
+    
+    Int_t relid[4] ;
+    
+    AbsToRelNumbering(id , relid) ;
+    
+    Int_t phosmodule = relid[0] ; 
+    
+    Float_t y0 = 0 ; 
+    
+    if ( relid[1] == 0 ) { // it is a PbW04 crystal 
+      y0 =  -(GetIPtoOuterCoverDistance() + GetUpperPlateThickness()
+             + GetSecondUpperPlateThickness() + GetUpperCoolingPlateThickness())  ;  
+    }
+    if ( relid[1] > 0 ) { // its a PPSD pad
+      if ( relid[1] >  GetNumberOfModulesPhi() *  GetNumberOfModulesZ() ) { // its an bottom module
+       y0 = -( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() / 2.0)  ;
+      } 
+      else // its an upper module
+       y0 = -( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() - GetLeadToMicro2Gap()
+               -  GetLeadConverterThickness() -  GetMicro1ToLeadGap() - GetMicromegas1Thickness() / 2.0) ; 
+    }
+    
+    Float_t x, z ; 
+    RelPosInModule(relid, x, z) ; 
+    
+    pos.SetX(x) ;
+    pos.SetZ(z) ;
+    pos.SetY( TMath::Sqrt(x*x + z*z + y0*y0) ) ; 
+    
+    
+    
+    Float_t phi           = GetPHOSAngle( phosmodule) ; 
+    Double_t const kRADDEG = 180.0 / kPI ;
+    Float_t rphi          = phi / kRADDEG ; 
+    
+    TRotation rot ;
+    rot.RotateZ(-rphi) ; // a rotation around Z by angle  
+    
+    TRotation dummy = rot.Invert() ;  // to transform from original frame to rotate frame
+    
+    pos.Transform(rot) ; // rotate the baby 
   }
   else {
- pos.SetX(0.);
- pos.SetY(0.);
- pos.SetZ(0.);
-       }
+    pos.SetX(0.);
+    pos.SetY(0.);
+    pos.SetZ(0.);
+  }
 } 
 
 //____________________________________________________________________________
@@ -448,34 +505,53 @@ void AliPHOSGeometry::RelPosInModule(const Int_t * relid, Float_t & x, Float_t &
 {
   // Converts the relative numbering into the local PHOS-module (x, z) coordinates
   // Note: sign of z differs from that in the previous version (Yu.Kharlov, 12 Oct 2000)
+
+  Bool_t padOfCPV  = (strcmp(fName,"IHEP")==0) ||
+                    ((strcmp(fName,"MIXT")==0) && relid[0]<=GetNCPVModules()) ;
+  Bool_t padOfPPSD = (strcmp(fName,"GPS2")==0) ||
+                    ((strcmp(fName,"MIXT")==0) && relid[0]> GetNCPVModules()) ;
   
   Int_t ppsdmodule  ; 
   Float_t x0,z0;
   Int_t row        = relid[2] ; //offset along x axiz
   Int_t column     = relid[3] ; //offset along z axiz
 
-  Float_t padsizeZ = GetPadSizeZ();
-  Float_t padsizeX = GetPadSizePhi();
-
+  Float_t padsizeZ = 0;
+  Float_t padsizeX = 0;
+  Int_t   nOfPadsPhi = 0;
+  Int_t   nOfPadsZ   = 0;
+  if      ( padOfPPSD ) {
+    padsizeZ   = GetPPSDModuleSize(2) / GetNumberOfPadsZ();
+    padsizeX   = GetPPSDModuleSize(0) / GetNumberOfPadsPhi();
+    nOfPadsPhi = GetNumberOfPadsPhi();
+    nOfPadsZ   = GetNumberOfPadsZ();
+  }
+  else if ( padOfCPV  ) {
+    padsizeZ   = GetPadSizeZ();
+    padsizeX   = GetPadSizePhi();
+    nOfPadsPhi = GetNumberOfCPVPadsPhi();
+    nOfPadsZ   = GetNumberOfCPVPadsZ();
+  }
+  
   if ( relid[1] == 0 ) { // its a PbW04 crystal 
     x = - ( GetNPhi()/2. - row    + 0.5 ) *  GetCrystalSize(0) ; // position ox Xtal with respect
     z =   ( GetNZ()  /2. - column + 0.5 ) *  GetCrystalSize(2) ; // of center of PHOS module  
   }  
   else  {    
-    if ( relid[1] >  GetNumberOfModulesPhi() *  GetNumberOfModulesZ() )
-      ppsdmodule =  relid[1]-GetNumberOfModulesPhi() *  GetNumberOfModulesZ(); 
-    else
-      ppsdmodule =  relid[1] ;
-    Int_t modrow = 1+(Int_t)TMath::Ceil( (Float_t)ppsdmodule / GetNumberOfModulesPhi()-1. ) ; 
-    Int_t modcol = ppsdmodule -  ( modrow - 1 ) * GetNumberOfModulesPhi() ;     
-    if ( ((strcmp( fName, "GPS2" ))  == 0) ) {
+    if ( padOfPPSD ) {
+      if ( relid[1] >  GetNumberOfModulesPhi() *  GetNumberOfModulesZ() )
+       ppsdmodule =  relid[1]-GetNumberOfModulesPhi() *  GetNumberOfModulesZ(); 
+      else
+       ppsdmodule =  relid[1] ;
+      Int_t modrow = 1+(Int_t)TMath::Ceil( (Float_t)ppsdmodule / GetNumberOfModulesPhi()-1. ) ; 
+      Int_t modcol = ppsdmodule -  ( modrow - 1 ) * GetNumberOfModulesPhi() ;     
       x0 = (  GetNumberOfModulesPhi() / 2.  - modrow  + 0.5 ) * GetPPSDModuleSize(0) ;
       z0 = (  GetNumberOfModulesZ()   / 2.  - modcol  + 0.5 ) * GetPPSDModuleSize(2)  ;     
     } else {
       x0 = 0;
       z0 = 0;
     }
-    x = - ( GetNumberOfPadsPhi()/2. - row    - 0.5 ) * padsizeX + x0 ; // position of pad  with respect
-    z =   ( GetNumberOfPadsZ()  /2. - column - 0.5 ) * padsizeZ - z0 ; // of center of PHOS module  
+    x = - ( nOfPadsPhi/2. - row    - 0.5 ) * padsizeX + x0 ; // position of pad  with respect
+    z =   ( nOfPadsZ  /2. - column - 0.5 ) * padsizeZ - z0 ; // of center of PHOS module  
   }
 }
index d62a8d3..0224831 100644 (file)
@@ -97,6 +97,8 @@ public:
   // Return general PHOS parameters
 
   Int_t   GetNModules(void)                    const { return fNModules ; }
+  Int_t   GetNPPSDModules(void)                const { return fNPPSDModules ; }
+  Int_t   GetNCPVModules(void)                 const { return fNModules - fNPPSDModules ; }
   Float_t GetPHOSAngle(Int_t index)            const { return fPHOSAngle[index-1] ; }
 
   // Return EMCA geometrical parameters
@@ -126,83 +128,53 @@ public:
 
   // Return PPSD geometrical parameters
 
-  Float_t GetAnodeThickness(void)              const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetAnodeThickness();         }
-  Float_t GetAvalancheGap(void)                const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetAvalancheGap();           }
-  Float_t GetCathodeThickness(void)            const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetCathodeThickness();       }
-  Float_t GetCompositeThickness(void)          const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetCompositeThickness();     }
-  Float_t GetConversionGap(void)               const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetConversionGap();          }
-  Float_t GetLeadConverterThickness(void)      const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetLeadConverterThickness(); }
-  Float_t GetLeadToMicro2Gap(void)             const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetLeadToMicro2Gap();        }
-  Float_t GetLidThickness(void)                const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetLidThickness();           }
-  Float_t GetMicromegas1Thickness(void)        const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetMicromegas1Thickness();   }
-  Float_t GetMicromegas2Thickness(void)        const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetMicromegas2Thickness();   }
-  Float_t GetMicromegasWallThickness(void)     const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetMicromegasWallThickness();}
-  Float_t GetMicro1ToLeadGap(void)             const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetMicro1ToLeadGap();        }
-  Float_t GetPCThickness(void)                 const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetPCThickness();            }
-  Float_t GetPhiDisplacement(void)             const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetPhiDisplacement();        }
-  Float_t GetPPSDModuleSize(Int_t index)       const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetPPSDModuleSize(index);    }
-  Float_t GetZDisplacement(void)               const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetZDisplacement();          }
+  Float_t GetAnodeThickness(void)              const { return fGeometryPPSD->GetAnodeThickness();         }
+  Float_t GetAvalancheGap(void)                const { return fGeometryPPSD->GetAvalancheGap();           }
+  Float_t GetCathodeThickness(void)            const { return fGeometryPPSD->GetCathodeThickness();       }
+  Float_t GetCompositeThickness(void)          const { return fGeometryPPSD->GetCompositeThickness();     }
+  Float_t GetConversionGap(void)               const { return fGeometryPPSD->GetConversionGap();          }
+  Float_t GetLeadConverterThickness(void)      const { return fGeometryPPSD->GetLeadConverterThickness(); }
+  Float_t GetLeadToMicro2Gap(void)             const { return fGeometryPPSD->GetLeadToMicro2Gap();        }
+  Float_t GetLidThickness(void)                const { return fGeometryPPSD->GetLidThickness();           }
+  Float_t GetMicromegas1Thickness(void)        const { return fGeometryPPSD->GetMicromegas1Thickness();   }
+  Float_t GetMicromegas2Thickness(void)        const { return fGeometryPPSD->GetMicromegas2Thickness();   }
+  Float_t GetMicromegasWallThickness(void)     const { return fGeometryPPSD->GetMicromegasWallThickness();}
+  Float_t GetMicro1ToLeadGap(void)             const { return fGeometryPPSD->GetMicro1ToLeadGap();        }
+  Int_t   GetNumberOfModulesPhi(void)          const { return fGeometryPPSD->GetNumberOfModulesPhi();     }
+  Int_t   GetNumberOfModulesZ(void)            const { return fGeometryPPSD->GetNumberOfModulesZ();       }
+  Int_t   GetNumberOfPadsPhi(void)             const { return fGeometryPPSD->GetNumberOfPadsPhi();        }
+  Int_t   GetNumberOfPadsZ(void)               const { return fGeometryPPSD->GetNumberOfPadsZ();          }
+  Float_t GetPCThickness(void)                 const { return fGeometryPPSD->GetPCThickness();            }
+  Float_t GetPhiDisplacement(void)             const { return fGeometryPPSD->GetPhiDisplacement();        }
+  Float_t GetPPSDModuleSize(Int_t index)       const { return fGeometryPPSD->GetPPSDModuleSize(index);    }
+  Float_t GetZDisplacement(void)               const { return fGeometryPPSD->GetZDisplacement();          }
+  void    SetLeadConverterThickness(Float_t x) const {        fGeometryPPSD->SetLeadConverterThickness(x);}
 
   // Return CPV geometrical parameters
 
-  Bool_t  IsLeadConverterExists(void)          const { return ((AliPHOSCPVGeometry*) fGeometryCPV)->IsLeadConverterExists();      }
-  Float_t GetCPVActiveSize(Int_t index)        const { return ((AliPHOSCPVGeometry*) fGeometryCPV)->GetCPVActiveSize(index);         }
-  Int_t   GetNumberOfCPVChipsPhi(void)         const { return ((AliPHOSCPVGeometry*) fGeometryCPV)->GetNumberOfCPVChipsPhi();     }
-  Int_t   GetNumberOfCPVChipsZ(void)           const { return ((AliPHOSCPVGeometry*) fGeometryCPV)->GetNumberOfCPVChipsZ();       }
-  Float_t GetGassiplexChipSize(Int_t index)    const { return ((AliPHOSCPVGeometry*) fGeometryCPV)->GetGassiplexChipSize(index);  }
-  Float_t GetCPVGasThickness(void)             const { return ((AliPHOSCPVGeometry*) fGeometryCPV)->GetCPVGasThickness();         }
-  Float_t GetCPVTextoliteThickness(void)       const { return ((AliPHOSCPVGeometry*) fGeometryCPV)->GetCPVTextoliteThickness();   }
-  Float_t GetCPVCuNiFoilThickness(void)        const { return ((AliPHOSCPVGeometry*) fGeometryCPV)->GetCPVCuNiFoilThickness();    }
-  Float_t GetFTPosition(Int_t index)           const { return ((AliPHOSCPVGeometry*) fGeometryCPV)->GetFTPosition(index);         }
-  Float_t GetCPVFrameSize(Int_t index)         const { return ((AliPHOSCPVGeometry*) fGeometryCPV)->GetCPVFrameSize(index);       }
+  Int_t   GetNumberOfCPVLayers(void)           const { return fGeometryCPV ->GetNumberOfCPVLayers();      }
+  Bool_t  IsLeadConverterExists(void)          const { return fGeometryCPV->IsLeadConverterExists();      }
+  Float_t GetCPVActiveSize(Int_t index)        const { return fGeometryCPV->GetCPVActiveSize(index);      }
+  Int_t   GetNumberOfCPVChipsPhi(void)         const { return fGeometryCPV->GetNumberOfCPVChipsPhi();     }
+  Int_t   GetNumberOfCPVChipsZ(void)           const { return fGeometryCPV->GetNumberOfCPVChipsZ();       }
+  Int_t   GetNumberOfCPVPadsPhi(void)          const { return fGeometryCPV->GetNumberOfCPVPadsPhi();      }
+  Int_t   GetNumberOfCPVPadsZ(void)            const { return fGeometryCPV->GetNumberOfCPVPadsZ();        }
+  Float_t GetPadSizePhi(void)                  const { return fGeometryCPV->GetCPVPadSizePhi();           }
+  Float_t GetPadSizeZ(void)                    const { return fGeometryCPV->GetCPVPadSizeZ();             }
+  Float_t GetGassiplexChipSize(Int_t index)    const { return fGeometryCPV->GetGassiplexChipSize(index);  }
+  Float_t GetCPVGasThickness(void)             const { return fGeometryCPV->GetCPVGasThickness();         }
+  Float_t GetCPVTextoliteThickness(void)       const { return fGeometryCPV->GetCPVTextoliteThickness();   }
+  Float_t GetCPVCuNiFoilThickness(void)        const { return fGeometryCPV->GetCPVCuNiFoilThickness();    }
+  Float_t GetFTPosition(Int_t index)           const { return fGeometryCPV->GetFTPosition(index);         }
+  Float_t GetCPVFrameSize(Int_t index)         const { return fGeometryCPV->GetCPVFrameSize(index);       }
 
   // Common PPSD and CPV parameters
 
-  Int_t   GetNumberOfCPVLayers(void)  const {
-    if      (strcmp(fName,"GPS2")==0) return 2;
-    else if (strcmp(fName,"IHEP")==0) return ((AliPHOSCPVGeometry*) fGeometryCPV)->GetNumberOfCPVLayers();
-    else                              return 0;
-  }
-
   Float_t GetCPVBoxSize(Int_t index)  const { 
-    if      (strcmp(fName,"GPS2")==0) return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetCPVBoxSize(index);
-    else if (strcmp(fName,"IHEP")==0) return ((AliPHOSCPVGeometry* ) fGeometryCPV)->GetCPVBoxSize(index);
-    else                              return 0;
-  }
-
-  Int_t   GetNumberOfModulesPhi(void) const {
-    if      (strcmp(fName,"GPS2")==0) return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetNumberOfModulesPhi();
-    else if (strcmp(fName,"IHEP")==0) return 1;
-    else                              return 0;
-  }
-
-  Int_t   GetNumberOfModulesZ(void)   const {
-    if      (strcmp(fName,"GPS2")==0) return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetNumberOfModulesZ();
-    else if (strcmp(fName,"IHEP")==0) return 1;
-    else                              return 0;
-  }
-
-  Int_t   GetNumberOfPadsPhi(void)    const { 
-    if      (strcmp(fName,"GPS2")==0) return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetNumberOfPadsPhi();
-    else if (strcmp(fName,"IHEP")==0) return ((AliPHOSCPVGeometry* ) fGeometryCPV)->GetNumberOfCPVPadsPhi();
-    else                              return 0;
-  }
-
-  Int_t   GetNumberOfPadsZ(void)      const { 
-    if      (strcmp(fName,"GPS2")==0) return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetNumberOfPadsZ();
-    else if (strcmp(fName,"IHEP")==0) return ((AliPHOSCPVGeometry* ) fGeometryCPV)->GetNumberOfCPVPadsZ();
-    else                              return 0;
-  }
-
-  Float_t GetPadSizePhi(void)         const {
-    if      (strcmp(fName,"GPS2")==0) return GetPPSDModuleSize(0) / GetNumberOfPadsPhi();
-    else if (strcmp(fName,"IHEP")==0) return ((AliPHOSCPVGeometry*) fGeometryCPV)->GetCPVPadSizePhi();
-    else                              return 0;
-  }
-
-  Float_t GetPadSizeZ(void)           const {
-    if      (strcmp(fName,"GPS2")==0) return GetPPSDModuleSize(2) / GetNumberOfPadsZ();
-    else if (strcmp(fName,"IHEP")==0) return ((AliPHOSCPVGeometry*) fGeometryCPV)->GetCPVPadSizeZ();
+    if      (strcmp(fName,"GPS2")==0) return fGeometryPPSD->GetCPVBoxSize(index);
+    else if (strcmp(fName,"IHEP")==0) return fGeometryCPV ->GetCPVBoxSize(index);
+    else if (strcmp(fName,"MIXT")==0) return TMath::Max(fGeometryCPV ->GetCPVBoxSize(index),
+                                                       fGeometryPPSD->GetCPVBoxSize(index));
     else                              return 0;
   }
 
@@ -246,13 +218,16 @@ protected:
 private:
 
   Int_t                    fNModules ;       // Number of modules constituing PHOS
+  Int_t                    fNPPSDModules ;   // Number of PPSD modules
+  Float_t                  fAngle ;          // Position angles between modules
   Float_t                 *fPHOSAngle ;      //[fNModules] Position angles of modules
   TObjArray               *fRotMatrixArray ; // Liste of rotation matrices (one per phos module)
   AliPHOSEMCAGeometry     *fGeometryEMCA ;   // Geometry object for Electromagnetic calorimeter
-  AliPHOSCPVBaseGeometry  *fGeometryCPV ;    // Geometry object for CPV (either GPS2 or IHEP)
+  AliPHOSCPVGeometry      *fGeometryCPV ;    // Geometry object for CPV  (IHEP)
+  AliPHOSPPSDGeometry     *fGeometryPPSD ;   // Geometry object for PPSD (GPS2)
   AliPHOSSupportGeometry  *fGeometrySUPP ;   // Geometry object for PHOS support
 
-  void                 SetPHOSAngles(); // calculates the PHOS modules PHI angle
+  void                     SetPHOSAngles();  // calculates the PHOS modules PHI angle
 
   static AliPHOSGeometry * fgGeom ; // pointer to the unique instance of the singleton 
   static Bool_t fgInit ;            // Tells if geometry has been succesfully set up 
index 866eff9..febee33 100644 (file)
@@ -75,7 +75,7 @@ Bool_t AliPHOSHit::operator==(AliPHOSHit const &rValue) const
 
   Bool_t rv = kFALSE ; 
 
-  if ( fId == rValue.GetId() && fPrimary == rValue.GetPrimary() ) 
+  if ( fId == rValue.GetId() && fPrimary == rValue.GetPrimary() )
     rv = kTRUE;
   
   return rv;
index d653289..936da55 100644 (file)
@@ -43,6 +43,10 @@ class AliPHOSHit : public AliHit {
     // return the identificator of this his
     return fId ; 
   }
+  Int_t   GetPid(void)      const { 
+    // return the particle PDG code which initiates this hit
+    return fPid ; 
+  }
   Int_t   GetPrimary(void)  const { 
     // returns the primary particle id at the origine of this hit 
     return fPrimary ; 
index d58a4ae..72c2c91 100644 (file)
@@ -22,6 +22,8 @@
 //*-- Author: Yves Schutz (SUBATECH)
 //////////////////////////////////////////////////////////////////////////////
 
+#include <iostream.h>
+
 // --- ROOT system ---
 
 #include "TTree.h"
index da5a23c..96904f7 100644 (file)
@@ -49,7 +49,7 @@ Float_t  AliPHOSPIDv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcclu,AliPHO
 {
   // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
  
-  Float_t r;
+  Float_t r = 1.e+10;
   TVector3 vecEmc ;
   TVector3 vecPpsd ;
   
index e581a2c..e3c4652 100644 (file)
@@ -47,6 +47,7 @@ public:
   virtual Float_t GetCPVBoxSize(Int_t index)       { return  fPPSDBoxSize[index] ;      }
   virtual Float_t GetPPSDModuleSize(Int_t index)   { return  fPPSDModuleSize[index] ;   }
   virtual Float_t GetZDisplacement(void)           { return  fZDisplacement ;           }
+  virtual void    SetLeadConverterThickness(Float_t x)  { fLeadConverterThickness = x;  }
   
   // CPV functions cannot be used for PPSD
   
index f32259e..d83fea1 100644 (file)
@@ -26,6 +26,8 @@
 
 // --- Standard library ---
 
+#include <assert.h>
+
 // --- AliRoot header files ---
 
 #include "AliPHOSRecParticle.h"
index 8d3c9f1..bc5716c 100644 (file)
@@ -8,6 +8,8 @@
 //                  
 //*-- Author: Gines Martinez (SUBATECH)
 
+#include <assert.h>
+
 // --- ROOT system ---
 
 #include "TMarker.h"
@@ -20,6 +22,8 @@
 
 #include "AliRecPoint.h"
 #include "AliPHOSDigit.h"
+#include "AliPHOSGeometry.h"
+class AliPHOSGeometry;
 
 class AliPHOSRecPoint : public AliRecPoint {
 
@@ -71,11 +75,11 @@ class AliPHOSRecPoint : public AliRecPoint {
     return *this ; 
   }
 
- protected:
+protected:
   
   Int_t fPHOSMod ;      // PHOS Module number in which the RecPoint is found
-
-
+  AliPHOSGeometry * fGeom ; // pointer to the PHOS geometry class
+  
   ClassDef(AliPHOSRecPoint,1) // RecPoint for PHOS (Base Class)
  
 };
index 772abdb..8e02248 100644 (file)
@@ -68,11 +68,11 @@ AliPHOSReconstructioner::AliPHOSReconstructioner(AliPHOSClusterizer * Clusterize
 } 
 
 //____________________________________________________________________________
- void AliPHOSReconstructioner::Make(DigitsList * dl, 
-                                   AliPHOSRecPoint::RecPointsList * emccl, 
-                                   AliPHOSRecPoint::RecPointsList * ppsdl, 
-                                   AliPHOSTrackSegment::TrackSegmentsList * trsl, 
-                                   AliPHOSRecParticle::RecParticlesList * rpl)
+ void AliPHOSReconstructioner::MakePPSD(DigitsList * dl, 
+                                       AliPHOSRecPoint::RecPointsList * emccl, 
+                                       AliPHOSRecPoint::RecPointsList * ppsdl, 
+                                       AliPHOSTrackSegment::TrackSegmentsList * trsl, 
+                                       AliPHOSRecParticle::RecParticlesList * rpl)
 {
   // Launches the Reconstruction process in the sequence: Make the reconstructed poins (clusterize)
   //                                                      Make the track segments 
@@ -118,8 +118,6 @@ AliPHOSReconstructioner::AliPHOSReconstructioner(AliPHOSClusterizer * Clusterize
     
   }
 
-
-
   // Making Clusters
   if  (fDebugReconstruction)  cout << "DebugReconstruction>>> Start making reconstructed points (clusterizing)" << endl;
 
@@ -337,9 +335,9 @@ AliPHOSReconstructioner::AliPHOSReconstructioner(AliPHOSClusterizer * Clusterize
 }
 
 //____________________________________________________________________________
- void AliPHOSReconstructioner::Make(DigitsList * dl, 
-                                   AliPHOSRecPoint::RecPointsList * emccl,
-                                   AliPHOSRecPoint::RecPointsList * cpvcl)
+ void AliPHOSReconstructioner::MakeCPV(DigitsList * dl, 
+                                      AliPHOSRecPoint::RecPointsList * emccl,
+                                      AliPHOSRecPoint::RecPointsList * cpvcl)
 {
 
   // Launches the Reconstruction process of EMC and CPV in the sequence:
@@ -358,7 +356,7 @@ AliPHOSReconstructioner::AliPHOSReconstructioner(AliPHOSClusterizer * Clusterize
   fClusterizer->MakeClusters(dl, emccl, cpvcl);
   
   if  (fDebugReconstruction){
-  // Digit Debuging
+    // Digit Debuging
     cout << "AliPHOSReconstructioner: Digit list entries are " << dl->GetEntries()    << endl ;
     cout << "AliPHOSReconstructioner: EMC   list entries are " << emccl->GetEntries() << endl ;
     cout << "AliPHOSReconstructioner: CPV   list entries are " << cpvcl->GetEntries() << endl ;
@@ -401,6 +399,7 @@ AliPHOSReconstructioner::AliPHOSReconstructioner(AliPHOSClusterizer * Clusterize
   AliPHOSCpvRecPoint * cpvrp ; 
   for (index = 0 ; index < cpvcl->GetEntries() ; index++) {
     cpvrp = (AliPHOSCpvRecPoint * )cpvcl->At(index) ; 
+    if (cpvrp->IsCPV()) break;
     cpvrp ->SetIndexInList(index) ; 
     TVector3  locpos;   cpvrp->GetLocalPosition(locpos);
     Int_t lengX,lengZ;  cpvrp->GetClusterLengths(lengX,lengZ);
index bd77a58..6f467d0 100644 (file)
@@ -42,14 +42,14 @@ public:
   AliPHOSClusterizer * GetClusterizer() { return fClusterizer ; }
   void Init(AliPHOSClusterizer * Clusterizer, AliPHOSTrackSegmentMaker * Tracker, 
                          AliPHOSPID * Identifier) ;  
-  void Make(TClonesArray * DL, 
-           AliPHOSRecPoint::RecPointsList * emccl, 
-           AliPHOSRecPoint::RecPointsList * ppsdl, 
-           AliPHOSTrackSegment::TrackSegmentsList * trsl, 
-           AliPHOSRecParticle::RecParticlesList * rpl)    ; // does the job for EMC+PPSD
-  void Make(TClonesArray * DL, 
-           AliPHOSRecPoint::RecPointsList * emccl, 
-           AliPHOSRecPoint::RecPointsList * ppsdl) ;        // does the job for EMC+CPV
+  void MakePPSD(TClonesArray * DL, 
+               AliPHOSRecPoint::RecPointsList * emccl, 
+               AliPHOSRecPoint::RecPointsList * ppsdl, 
+               AliPHOSTrackSegment::TrackSegmentsList * trsl, 
+               AliPHOSRecParticle::RecParticlesList * rpl)    ; // does the job for EMC+PPSD
+  void MakeCPV (TClonesArray * DL, 
+               AliPHOSRecPoint::RecPointsList * emccl, 
+               AliPHOSRecPoint::RecPointsList * ppsdl) ;        // does the job for EMC+CPV
 
   void SetDebugReconstruction(Bool_t deb) { fDebugReconstruction = deb; }
 
index a5d0d54..9d2a180 100644 (file)
@@ -382,7 +382,7 @@ void  AliPHOSTrackSegmentMakerv1::MakeTrackSegments(DigitsList * dl,
 {
   // Makes the track segments out of the list of EMC and PPSD Recpoints and stores them in a list
   
-  Int_t phosmod      = 1 ;
+  Int_t phosmod      = fGeom->GetNCPVModules() + 1 ;
   Int_t emcStopedAt  = 0 ; 
   Int_t ppsdStopedAt = 0 ; 
 
index 4149719..28a6c53 100644 (file)
@@ -54,7 +54,9 @@ AliPHOSv0::AliPHOSv0(const char *name, const char *title):
   AliPHOS(name,title)
 {
   // ctor : title is used to identify the layout
-  //        GPS2 = 5 modules (EMC + PPSD)   
+  //        GPS2 = 5 modules (EMC + PPSD)
+  //        IHEP = 5 modules (EMC + CPV)
+  //        MIXT = 4 modules (EMC + CPV) and 1 module (EMC + PPSD)
  
   // gets an instance of the geometry parameters class  
 
@@ -98,12 +100,17 @@ void AliPHOSv0::BuildGeometry()
   //END_HTML  
 
   this->BuildGeometryforPHOS() ; 
-  if      ( ( strcmp(fGeom->GetName(), "GPS2" ) == 0 ) ) 
+  if      (strcmp(fGeom->GetName(),"GPS2") == 0)
     this->BuildGeometryforPPSD() ;
-  else if ( ( strcmp(fGeom->GetName(), "IHEP" ) == 0 ) ) 
+  else if (strcmp(fGeom->GetName(),"IHEP") == 0)
     this->BuildGeometryforCPV() ;
+  else if (strcmp(fGeom->GetName(),"MIXT") == 0) {
+    this->BuildGeometryforPPSD() ;
+    this->BuildGeometryforCPV() ;
+  }
   else
-    cout << "AliPHOSv0::BuildGeometry : no charged particle identification system installed" << endl; 
+    cout << "AliPHOSv0::BuildGeometry : no charged particle identification system installed: "
+        << "Geometry name = " << fGeom->GetName() << endl; 
 
 }
 
@@ -310,9 +317,15 @@ void AliPHOSv0:: BuildGeometryforPPSD(void)
   Int_t number = 988 ; 
   TNode * top = gAlice->GetGeometry()->GetNode("alice") ;
  
-  for( Int_t i = 1; i <= fGeom->GetNModules(); i++ ) { // the number of PHOS modules
+  Int_t firstModule = 0 ; 
+  if      (strcmp(fGeom->GetName(),"GPS2") == 0) 
+    firstModule = 1;
+  else if (strcmp(fGeom->GetName(),"MIXT") == 0) 
+    firstModule = fGeom->GetNModules() - fGeom->GetNPPSDModules() + 1;
+  
+  for( Int_t i = firstModule; i <= fGeom->GetNModules(); i++ ) { // the number of PHOS modules
     Float_t angle = fGeom->GetPHOSAngle(i) ;
-    sprintf(rotname, "%s%d", "rotg", number++) ;
+    sprintf(rotname, "%s%d", "rotg", number+i) ;
     new TRotMatrix(rotname, rotname, 90, angle, 90, 90 + angle, 0, 0);
     top->cd();
     sprintf(nodename, "%s%d", "Moduleg", i) ;    
@@ -535,12 +548,19 @@ void AliPHOSv0:: BuildGeometryforCPV(void)
   Float_t r = fGeom->GetIPtoCPVDistance() + fGeom->GetCPVBoxSize(1) / 2.0 ;
   Int_t number = 988 ; 
   TNode * top = gAlice->GetGeometry()->GetNode("alice") ;
-  for( Int_t i = 1; i <= fGeom->GetNModules(); i++ ) { // the number of PHOS modules
+
+  Int_t lastModule = 0 ;
+  if      (strcmp(fGeom->GetName(),"IHEP") == 0) 
+    lastModule = fGeom->GetNModules();
+  else if (strcmp(fGeom->GetName(),"MIXT") == 0) 
+    lastModule = fGeom->GetNModules() - fGeom->GetNPPSDModules();
+  
+  for( Int_t i = 1; i <= lastModule; i++ ) { // the number of PHOS modules
 
     // One CPV module
 
     Float_t angle = fGeom->GetPHOSAngle(i) ;
-    sprintf(rotname, "%s%d", "rotg", number++) ;
+    sprintf(rotname, "%s%d", "rotg", number+i) ;
     new TRotMatrix(rotname, rotname, 90, angle, 90, 90 + angle, 0, 0);
     top->cd();
     sprintf(nodename, "%s%d", "CPVModule", i) ;    
@@ -617,41 +637,73 @@ void AliPHOSv0::CreateGeometry()
   // Get pointer to the array containing media indeces
   Int_t *idtmed = fIdtmed->GetArray() - 699 ;
 
+  // Create a box a PHOS module.
+  // In case of MIXT geometry 2 different boxes are needed
+
   Float_t bigbox[3] ; 
   bigbox[0] =   fGeom->GetOuterBoxSize(0) / 2.0 ;
   bigbox[1] = ( fGeom->GetOuterBoxSize(1) + fGeom->GetCPVBoxSize(1) ) / 2.0 ;
   bigbox[2] =   fGeom->GetOuterBoxSize(2) / 2.0 ;
   
-  gMC->Gsvolu("PHOS", "BOX ", idtmed[798], bigbox, 3) ;
+    gMC->Gsvolu("PHOS", "BOX ", idtmed[798], bigbox, 3) ;
+
+  if ( strcmp( fGeom->GetName(),"MIXT") == 0 && fGeom->GetNPPSDModules() > 0) 
+    gMC->Gsvolu("PHO1", "BOX ", idtmed[798], bigbox, 3) ;
   
-  this->CreateGeometryforPHOS() ; 
+    this->CreateGeometryforPHOS() ; 
   if      ( strcmp( fGeom->GetName(), "GPS2") == 0  ) 
     this->CreateGeometryforPPSD() ;
   else if ( strcmp( fGeom->GetName(), "IHEP") == 0  ) 
     this->CreateGeometryforCPV() ;
+  else if ( strcmp( fGeom->GetName(), "MIXT") == 0  ) {
+    this->CreateGeometryforPPSD() ;
+    this->CreateGeometryforCPV() ;
+  }
   else
     cout << "AliPHOSv0::CreateGeometry : no charged particle identification system installed" << endl; 
-  this->CreateGeometryforSupport() ; 
+
+//      this->CreateGeometryforSupport() ; 
   
   // --- Position  PHOS mdules in ALICE setup ---
   
   Int_t idrotm[99] ;
   Double_t const kRADDEG = 180.0 / kPI ;
   
-  for( Int_t i = 1; i <= fGeom->GetNModules(); i++ ) {
+  Int_t lastModule;
+  if (strcmp(fGeom->GetName(),"MIXT") == 0) 
+    lastModule = fGeom->GetNModules() - fGeom->GetNPPSDModules();
+  else
+    lastModule = fGeom->GetNModules();
+
+  Int_t i;
+  for( i = 1; i <= lastModule ; i++ ) {
     
     Float_t angle = fGeom->GetPHOSAngle(i) ;
     AliMatrix(idrotm[i-1], 90.0, angle, 90.0, 90.0+angle, 0.0, 0.0) ;
  
     Float_t r = fGeom->GetIPtoOuterCoverDistance() + ( fGeom->GetOuterBoxSize(1) + fGeom->GetCPVBoxSize(1) ) / 2.0 ;
 
-    Float_t xP1 = r * TMath::Sin( angle / kRADDEG ) ;
+    Float_t xP1 =  r * TMath::Sin( angle / kRADDEG ) ;
     Float_t yP1 = -r * TMath::Cos( angle / kRADDEG ) ;
 
     gMC->Gspos("PHOS", i, "ALIC", xP1, yP1, 0.0, idrotm[i-1], "ONLY") ;
  
   } // for GetNModules
 
+  for( i = lastModule+1; i <= fGeom->GetNModules(); i++ ) {
+    
+    Float_t angle = fGeom->GetPHOSAngle(i) ;
+    AliMatrix(idrotm[i-1], 90.0, angle, 90.0, 90.0+angle, 0.0, 0.0) ;
+    Float_t r = fGeom->GetIPtoOuterCoverDistance() + ( fGeom->GetOuterBoxSize(1) + fGeom->GetCPVBoxSize(1) ) / 2.0 ;
+
+    Float_t xP1 =  r * TMath::Sin( angle / kRADDEG ) ;
+    Float_t yP1 = -r * TMath::Cos( angle / kRADDEG ) ;
+
+    gMC->Gspos("PHO1", i-lastModule, "ALIC", xP1, yP1, 0.0, idrotm[i-1], "ONLY") ;
+  } // for GetNModules
+
 }
 
 //____________________________________________________________________________
@@ -686,7 +738,9 @@ void AliPHOSv0::CreateGeometryforPHOS()
 
   Float_t yO =  - fGeom->GetCPVBoxSize(1)  / 2.0 ;
 
-  gMC->Gspos("EMCA", 1, "PHOS", 0.0, yO, 0.0, 0, "ONLY") ; 
+    gMC->Gspos("EMCA", 1, "PHOS", 0.0, yO, 0.0, 0, "ONLY") ; 
+  if ( strcmp( fGeom->GetName(),"MIXT") == 0 && fGeom->GetNPPSDModules() > 0) 
+    gMC->Gspos("EMCA", 1, "PHO1", 0.0, yO, 0.0, 0, "ONLY") ; 
 
   // ---
   // --- Define Textolit Wall box, position inside EMCA ---
@@ -950,7 +1004,10 @@ void AliPHOSv0::CreateGeometryforPPSD()
 
   Float_t yO =  fGeom->GetOuterBoxSize(1) / 2.0 ;
 
-  gMC->Gspos("PPSD", 1, "PHOS", 0.0, yO, 0.0, 0, "ONLY") ; 
+  if ( strcmp( fGeom->GetName(),"MIXT") == 0 && fGeom->GetNPPSDModules() > 0) 
+    gMC->Gspos("PPSD", 1, "PHO1", 0.0, yO, 0.0, 0, "ONLY") ; 
+  else
+    gMC->Gspos("PPSD", 1, "PHOS", 0.0, yO, 0.0, 0, "ONLY") ; 
 
   // Now we build a micromegas module
   // The box containing the whole module filled with epoxy (FR4)
@@ -1013,7 +1070,7 @@ void AliPHOSv0::CreateGeometryforPPSD()
   
   // --- Divide GGPP in X (phi) and Z directions --
   gMC->Gsdvn("GROW", "GGPS", fGeom->GetNumberOfPadsPhi(), 1) ;
-  gMC->Gsdvn("GCEL", "GROW", fGeom->GetNumberOfPadsZ() , 3) ;
+  gMC->Gsdvn("GCEL", "GROW", fGeom->GetNumberOfPadsZ() ,  3) ;
 
   y0 = y0 - fGeom->GetAnodeThickness() / 2.  - ( fGeom->GetConversionGap() +  fGeom->GetAvalancheGap() ) / 2. ; 
 
@@ -1388,6 +1445,22 @@ void AliPHOSv0::CreateGeometryforSupport()
 }
 
 //____________________________________________________________________________
+Float_t AliPHOSv0::ZMin(void) const
+{
+  // Overall dimension of the PHOS (min)
+  // Take it twice more than the PHOS module size
+  return -fGeom->GetOuterBoxSize(2);
+}
+
+//____________________________________________________________________________
+Float_t AliPHOSv0::ZMax(void) const
+{
+  // Overall dimension of the PHOS (max)
+  // Take it twice more than the PHOS module size
+  return  fGeom->GetOuterBoxSize(2);
+}
+
+//____________________________________________________________________________
 void AliPHOSv0::Init(void)
 {
   // Just prints an information message
index 5a382ca..0ff4936 100644 (file)
@@ -42,15 +42,18 @@ class AliPHOSv0 : public AliPHOS {
     // useless since there are no hits
     assert(0==1) ; 
   }
-  virtual void   BuildGeometry(void) ;                              // creates the geometry for the ROOT display
-  void           BuildGeometryforPHOS(void) ;                       // creates the PHOS geometry for the ROOT display
-  void           BuildGeometryforPPSD(void) ;                       // creates the PPSD geometry for the ROOT display
-  void           BuildGeometryforCPV(void) ;                        // creates the CPV  geometry for the ROOT display
-  virtual void   CreateGeometry(void) ;                             // creates the geometry for GEANT
-  void           CreateGeometryforPHOS(void) ;                      // creates the PHOS geometry for GEANT
-  void           CreateGeometryforPPSD(void) ;                      // creates the PPSD geometry for GEANT
-  void           CreateGeometryforCPV(void) ;                       // creates the CPV  geometry for GEANT
-  void           CreateGeometryforSupport(void) ;                   // creates the Support geometry for GEANT
+  virtual void   BuildGeometry(void) ;             // creates the geometry for the ROOT display
+  void           BuildGeometryforPHOS(void) ;      // creates the PHOS geometry for the ROOT display
+  void           BuildGeometryforPPSD(void) ;      // creates the PPSD geometry for the ROOT display
+  void           BuildGeometryforCPV(void) ;       // creates the CPV  geometry for the ROOT display
+  virtual void   CreateGeometry(void) ;            // creates the geometry for GEANT
+  void           CreateGeometryforPHOS(void) ;     // creates the PHOS geometry for GEANT
+  void           CreateGeometryforPPSD(void) ;     // creates the PPSD geometry for GEANT
+  void           CreateGeometryforCPV(void) ;      // creates the CPV  geometry for GEANT
+  void           CreateGeometryforSupport(void) ;  // creates the Support geometry for GEANT
+  virtual Float_t ZMin() const;                    // overall dimension of the module (min)
+  virtual Float_t ZMax() const;                    // overall dimension of the module (max)
+
   virtual AliPHOSGeometry * GetGeometry() {
     // gets the pointer to the AliPHOSGeometry unique instance  
     return fGeom ; 
index da91a13..dcf9bb0 100644 (file)
 // Implementation version v1 of PHOS Manager class 
 //---
 // Layout EMC + PPSD has name GPS2:
-// Produces cumulated hits (no hits) and digits                  
+// Produces cumulated hits
 //---
 // Layout EMC + CPV  has name IHEP:
-// Produces hits for CPV, cumulated hits and digits                  
+// Produces hits for CPV, cumulated hits
+//---
+// Layout EMC + CPV + PPSD has name GPS:
+// Produces hits for CPV, cumulated hits
+//---
 //*-- Author: Yves Schutz (SUBATECH)
 
 
@@ -57,13 +61,19 @@ ClassImp(AliPHOSv1)
 AliPHOSv1::AliPHOSv1()
 {
   // ctor
-  fNTmpHits = 0 ; 
-  fTmpHits  = 0 ;
 
   // Create an empty array of AliPHOSCPVModule to satisfy
   // AliPHOSv1::Streamer when reading root file
+  fReconstructioner  = 0;
+  fTrackSegmentMaker = 0;
+  if ( 0==(fEMCModules=new TClonesArray("AliPHOSCPVModule",0)) ) {
+    Error("AliPHOSv1","Can not create array of EMC modules");
+    exit(1);
+  }
 
-  if ( NULL==(fCPVModules=new TClonesArray("AliPHOSCPVModule",0)) ) {
+  if ( 0==(fCPVModules=new TClonesArray("AliPHOSCPVModule",0)) ) {
     Error("AliPHOSv1","Can not create array of CPV modules");
     exit(1);
   }
@@ -77,14 +87,13 @@ AliPHOSv0(name,title)
   // ctor : title is used to identify the layout
   //        GPS2 = 5 modules (EMC + PPSD)
   //        IHEP = 5 modules (EMC + CPV )
-  // We use 2 arrays of hits :
+  //        MIXT = 4 modules (EMC + CPV ) and 1 module (EMC + PPSD)
   //
+  // We store hits :
   //   - fHits (the "normal" one), which retains the hits associated with
   //     the current primary particle being tracked
   //     (this array is reset after each primary has been tracked).
   //
-  //   - fTmpHits, which retains all the hits of the current event. It 
-  //     is used for the digitization part.
 
   fPinElectronicNoise = 0.010 ;
   fDigitThreshold      = 0.1 ;   // 1 GeV 
@@ -94,25 +103,35 @@ AliPHOSv0(name,title)
   // It is put in the Digit Tree because the TreeH is filled after each primary
   // and the TreeD at the end of the event (branch is set in FinishEvent() ).
   
-  fTmpHits= new TClonesArray("AliPHOSHit",1000) ;
+  fHits= new TClonesArray("AliPHOSHit",1000) ;
 
-  fNTmpHits = fNhits = 0 ;
+  fNhits = 0 ;
 
-  fDigits = new TClonesArray("AliPHOSDigit",1000) ;
+  fReconstructioner  = 0;
+  fTrackSegmentMaker = 0;
 
   fIshunt     =  1 ; // All hits are associated with primary particles
   
+  // 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();
+
   // Create array of CPV modules for the IHEP's version of CPV
 
-  if ( strcmp(fGeom->GetName(),"IHEP") == 0 ) {
-    // Create array of AliPHOSCPVmodule 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 ( NULL==(fCPVModules=new TClonesArray("AliPHOSCPVModule",fGeom->GetNModules())) ) {
+    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->GetNModules(); i++) new(lcpvmodule[i]) AliPHOSCPVModule();
+    for (Int_t i=0; i<fGeom->GetNCPVModules(); i++) new(lcpvmodule[i]) AliPHOSCPVModule();
   }
   else {
     // Create an empty array of AliPHOSCPVModule to satisfy
@@ -143,10 +162,10 @@ AliPHOSv1::AliPHOSv1(AliPHOSReconstructioner * Reconstructioner, const char *nam
   // We do not want to save in TreeH the raw hits
   //fHits   = new TClonesArray("AliPHOSHit",100) ;
 
-  fDigits = new TClonesArray("AliPHOSDigit",1000) ;
-  fTmpHits= new TClonesArray("AliPHOSHit",1000) ;
+  fDigits = 0 ;
+  fHits= new TClonesArray("AliPHOSHit",1000) ;
 
-  fNTmpHits = fNhits = 0 ;
+  fNhits = 0 ;
 
   fIshunt     =  1 ; // All hits are associated with primary particles
  
@@ -169,30 +188,36 @@ AliPHOSv1::~AliPHOSv1()
 {
   // dtor
 
-  if ( fTmpHits) {
-    fTmpHits->Delete() ; 
-    delete fTmpHits ;
-    fTmpHits = 0 ; 
+  if ( fHits) {
+    fHits->Delete() ; 
+    delete fHits ;
+    fHits = 0 ; 
   }
 
-//   if ( fEmcRecPoints ) {
-//     fEmcRecPoints->Delete() ; 
-//     delete fEmcRecPoints ; 
-//     fEmcRecPoints = 0 ; 
-//   }
-
-//   if ( fPpsdRecPoints ) { 
-//     fPpsdRecPoints->Delete() ;
-//     delete fPpsdRecPoints ;
-//     fPpsdRecPoints = 0 ; 
-//   }
-  
-//   if ( fTrackSegments ) {
-//     fTrackSegments->Delete() ; 
-//     delete fTrackSegments ;
-//     fTrackSegments = 0 ; 
-//   }
+  if ( fDigits) {
+    fDigits->Delete() ; 
+    delete fDigits ;
+    fDigits = 0 ; 
+  }
 
+  if ( fEmcRecPoints ) {
+    fEmcRecPoints->Delete() ; 
+    delete fEmcRecPoints ; 
+    fEmcRecPoints = 0 ; 
+  }
+  
+  if ( fPpsdRecPoints ) { 
+    fPpsdRecPoints->Delete() ;
+    delete fPpsdRecPoints ;
+    fPpsdRecPoints = 0 ; 
+  }
+  
+  if ( fTrackSegments ) {
+    fTrackSegments->Delete() ; 
+    delete fTrackSegments ;
+    fTrackSegments = 0 ; 
+  }
+  
 }
 
 //____________________________________________________________________________
@@ -203,41 +228,26 @@ void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id,
   //   or in a single PPSD gas cell
 
   Int_t hitCounter ;
-  TClonesArray &ltmphits = *fTmpHits ;
   AliPHOSHit *newHit ;
   AliPHOSHit *curHit ;
   Bool_t deja = kFALSE ;
 
-  // In any case, fills the fTmpHit TClonesArray (with "accumulated hits")
-
   newHit = new AliPHOSHit(shunt, primary, tracknumber, Id, hits, trackpid) ;
 
-  // We do not want to save in TreeH the raw hits 
-  //  TClonesArray &lhits = *fHits;
-
-  for ( hitCounter = 0 ; hitCounter < fNTmpHits && !deja ; hitCounter++ ) {
-    curHit = (AliPHOSHit*) ltmphits[hitCounter] ;
-  if( *curHit == *newHit ) {
-    *curHit = *curHit + *newHit ;
-    deja = kTRUE ;
+  for ( hitCounter = 0 ; hitCounter < fNhits && !deja ; hitCounter++ ) {
+    curHit = (AliPHOSHit*) (*fHits)[hitCounter] ;
+    if( *curHit == *newHit ) {
+      *curHit = *curHit + *newHit ;
+      deja = kTRUE ;
     }
   }
          
   if ( !deja ) {
-    new(ltmphits[fNTmpHits]) AliPHOSHit(*newHit) ;
-    fNTmpHits++ ;
+    new((*fHits)[fNhits]) AliPHOSHit(*newHit) ;
+    fNhits++ ;
   }
 
-  // We do not want to save in TreeH the raw hits 
-  //   new(lhits[fNhits]) AliPHOSHit(*newHit) ;    
-  //   fNhits++ ;
-
-  // Please note that the fTmpHits array must survive up to the
-  // end of the events, so it does not appear e.g. in ResetHits() (
-  // which is called at the end of each primary).  
-
   delete newHit;
-
 }
 
 //___________________________________________________________________________
@@ -251,114 +261,117 @@ Int_t AliPHOSv1::Digitize(Float_t Energy)
   return chan ;
 }
 
-//___________________________________________________________________________
-void AliPHOSv1::FinishEvent()
-{
-  // Makes the digits from the sum of summed hit in a single crystal or PPSD gas cell
-  // Adds to the energy the electronic noise
-  // Keeps digits with energy above fDigitThreshold
-
-  // Save the cumulated hits instead of raw hits (need to create the branch myself)
-  // It is put in the Digit Tree because the TreeH is filled after each primary
-  // and the TreeD at the end of the event.
+//____________________________________________________________________________
+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);
+
   
   Int_t i ;
   Int_t relid[4];
   Int_t j ; 
-  TClonesArray &lDigits = *fDigits ;
   AliPHOSHit   * hit ;
   AliPHOSDigit * newdigit ;
   AliPHOSDigit * curdigit ;
   Bool_t deja = kFALSE ; 
   
-  for ( i = 0 ; i < fNTmpHits ; i++ ) {
-    hit = (AliPHOSHit*)fTmpHits->At(i) ;
-
-    // Assign primary number only if contribution is significant
-    if( hit->GetEnergy() > fDigitThreshold)
-      newdigit = new AliPHOSDigit( hit->GetPrimary(), hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
-    else
-      newdigit = new AliPHOSDigit( -1               , hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
-    deja =kFALSE ;
-    for ( j = 0 ; j < fNdigits ;  j++) { 
-      curdigit = (AliPHOSDigit*) lDigits[j] ;
-      if ( *curdigit == *newdigit) {
-       *curdigit = *curdigit + *newdigit ; 
-       deja = kTRUE ; 
+  Int_t itrack ;
+  for (itrack=0; itrack<ntracks; 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) ;
+    
+      // Assign primary number only if contribution is significant
+      if( hit->GetEnergy() > fDigitThreshold)
+       newdigit = new AliPHOSDigit( hit->GetPrimary(), hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
+      else
+       newdigit = new AliPHOSDigit( -1               , hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
+      deja =kFALSE ;
+
+
+      for ( j = 0 ; j < fNdigits ;  j++) { 
+       curdigit = (AliPHOSDigit*) fDigits->At(j) ;
+       if ( *curdigit == *newdigit) {
+         *curdigit = *curdigit + *newdigit ; 
+         deja = kTRUE ; 
+       }
       }
-    }
-    if ( !deja ) {
-      new(lDigits[fNdigits]) AliPHOSDigit(* newdigit) ;
-      fNdigits++ ;  
-    }
-    delete newdigit ;    
-  } 
-  
-  // Noise induced by the PIN diode of the PbWO crystals
 
+      if ( !deja ) {
+       new((*fDigits)[fNdigits]) AliPHOSDigit(* newdigit) ;
+       fNdigits++ ;  
+      }
+      
+      delete newdigit ;    
+    } 
+
+  } // loop over tracks
+    
+  // Noise induced by the PIN diode of the PbWO crystals
+  
   Float_t energyandnoise ;
   for ( i = 0 ; i < fNdigits ; i++ ) {
     newdigit =  (AliPHOSDigit * ) fDigits->At(i) ;
-
+    
     fGeom->AbsToRelNumbering(newdigit->GetId(), relid) ;
-
+    
     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) ; 
     }
   }
-  
+    
   fDigits->Compress() ;  
+  
+  fNdigits = fDigits->GetEntries() ;
+  fDigits->Expand(fNdigits) ;
 
-  fNdigits =  fDigits->GetEntries() ; 
   for (i = 0 ; i < fNdigits ; i++) { 
     newdigit = (AliPHOSDigit *) fDigits->At(i) ; 
-    newdigit->SetIndexInList(i) ; 
-
-//      fGeom->AbsToRelNumbering(newdigit->GetId(), relid) ;
-//      printf("FinishEvent(): relid=(%d,%d,%d,%d) Amp=%d\n",
-//        relid[0],relid[1],relid[2],relid[3], newdigit->GetAmp());
+    newdigit->SetIndexInList(i) ;     
   }
 
-}
+  gAlice->TreeD()->Fill() ;
 
+  gAlice->TreeD()->Write(0,TObject::kOverwrite) ;
+}
 //___________________________________________________________________________
 void AliPHOSv1::MakeBranch(Option_t* opt)
 {  
   // Create new branche in the current Root Tree in the digit Tree
   AliDetector::MakeBranch(opt) ;
   
-  char branchname[10];
-
-  // Create new branche PHOSCH in the current Root Tree in the Hit Tree for accumulated Hits
-  if ( ! (gAlice->IsLegoRun()) ) { // only when not in lego plot mode 
-    char *cdH = strstr(opt,"H");
-    if ( fTmpHits && gAlice->TreeH()  && cdH) {
-      sprintf(branchname, "%sCH", GetName()) ;
-      gAlice->TreeH()->Branch(branchname, &fTmpHits, fBufferSize) ;
-    }   
-  }
-  sprintf(branchname,"%s",GetName());
-  char *cdD = strstr(opt,"D");
-  if (fDigits && gAlice->TreeD() && cdD) {
-    gAlice->TreeD()->Branch(branchname, &fDigits, fBufferSize);
-  }
-
-
+  // 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);
+  
   // Create new branches CPV<i> for hits in CPV modules for IHEP geometry
-  // Yuri Kharlov, 28 September 2000.
-
-  if ( strcmp(fGeom->GetName(),"IHEP") == 0 ) {
-    for( Int_t i=0; i<fGeom->GetNModules(); i++ ) GetCPVModule(i).MakeBranch(i+1);
+  
+  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);
   }
-
+  
 }
 
 //_____________________________________________________________________________
@@ -402,7 +415,7 @@ void AliPHOSv1::Reconstruction(AliPHOSReconstructioner * Reconstructioner)
 
   fRecParticles->Delete() ; 
 
-  if      (strcmp(fGeom->GetName(),"GPS2") == 0) {
+  if      (strcmp(fGeom->GetName(),"GPS2") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0) {
     if ( fRecParticles && gAlice->TreeR() ) { 
       sprintf(branchname,"%sRP",GetName()) ;
       gAlice->TreeR()->Branch(branchname, &fRecParticles, fBufferSize) ;
@@ -410,10 +423,10 @@ void AliPHOSv1::Reconstruction(AliPHOSReconstructioner * Reconstructioner)
   }
   
   // 3.
-  if      (strcmp(fGeom->GetName(),"GPS2") == 0)
-    fReconstructioner->Make(fDigits, fEmcRecPoints, fPpsdRecPoints, fTrackSegments, fRecParticles);
-  else if (strcmp(fGeom->GetName(),"IHEP") == 0)
-    fReconstructioner->Make(fDigits, fEmcRecPoints, fPpsdRecPoints);
+  if      (strcmp(fGeom->GetName(),"GPS2") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0)
+    fReconstructioner->MakePPSD(fDigits, fEmcRecPoints, fPpsdRecPoints, fTrackSegments, fRecParticles);
+  if (strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0)
+    fReconstructioner->MakeCPV (fDigits, fEmcRecPoints, fPpsdRecPoints);
 
   // 4. Expand or Shrink the arrays to the proper size
   
@@ -421,7 +434,7 @@ void AliPHOSv1::Reconstruction(AliPHOSReconstructioner * Reconstructioner)
   
   size = fEmcRecPoints->GetEntries() ;
   fEmcRecPoints->Expand(size) ;
+
   size = fPpsdRecPoints->GetEntries() ;
   fPpsdRecPoints->Expand(size) ;
 
@@ -444,29 +457,15 @@ void AliPHOSv1::Reconstruction(AliPHOSReconstructioner * Reconstructioner)
 //____________________________________________________________________________
 void AliPHOSv1::ResetHits() 
 {              
-  // Reset hit tree for CPV in IHEP geometry
+  // Reset hit tree for EMC and CPV
   // Yuri Kharlov, 28 September 2000
 
-
-  if ( fTmpHits ) {
-    fTmpHits->Clear() ;
-    fNTmpHits = 0;
-  }
-
   AliDetector::ResetHits();
-  if ( strcmp(fGeom->GetName(),"IHEP") == 0 ) {
-    for (Int_t i=0; i<fGeom->GetNModules(); i++) ((AliPHOSCPVModule*)(*fCPVModules)[i]) -> Clear();
-  }
-}  
-//____________________________________________________________________________
-void AliPHOSv1::ResetDigits() 
-{ 
-  // May sound strange, but cumulative hits are store in digits Tree
-  AliDetector::ResetDigits();
-  if(  fTmpHits ) {
-    fTmpHits->Delete();
-    fNTmpHits = 0 ;
+  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();
   }
+  
 }  
 //____________________________________________________________________________
 void AliPHOSv1::ResetReconstruction() 
@@ -486,18 +485,23 @@ void AliPHOSv1::SetTreeAddress()
   //  TBranch *branch;
   AliPHOS::SetTreeAddress();
 
- //  //Branch address for TreeR: RecPpsdRecPoint
+//  //Branch address for TreeR: RecPpsdRecPoint
 //   TTree *treeR = gAlice->TreeR();
 //   if ( treeR && fPpsdRecPoints ) {
 //     branch = treeR->GetBranch("PHOSPpsdRP");
 //     if (branch) branch->SetAddress(&fPpsdRecPoints) ;
 //  }
 
+  // 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);
+
   // 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 ) {
-    for( Int_t i=0; i<fGeom->GetNModules(); i++ ) GetCPVModule(i).SetTreeAddress(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);
   }
 
 }
@@ -508,17 +512,23 @@ void AliPHOSv1::StepManager(void)
 {
   // Accumulates hits as long as the track stays in a single crystal or PPSD gas Cell
 
+//    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
   Int_t          absid    ;      // absolute cell ID number
   Float_t        xyze[4]  ;      // position wrt MRS and energy deposited
   TLorentzVector pos      ;      // Lorentz vector of the track current position
   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() ; 
-  if ( name == "GPS2" ) {                                       // ======> CPV is a GPS' PPSD
+  if ( name == "GPS2" || name == "MIXT" ) {            // ======> CPV is a GPS' PPSD
 
     if( gMC->CurrentVolID(copy) == gMC->VolId("GCEL") ) // We are inside a gas cell 
     {
@@ -530,9 +540,11 @@ void AliPHOSv1::StepManager(void)
 
       if ( xyze[3] != 0 ) { // there is deposited energy 
                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 
 
@@ -547,7 +559,7 @@ void AliPHOSv1::StepManager(void)
     } // We are inside the gas of the CPV  
   } // GPS2 configuration
 
-  else if ( name == "IHEP" ) {                                  // ======> CPV is a IHEP's one
+  if ( name == "IHEP" || name == "MIXT" ) {       // ======> CPV is a IHEP's one
 
     // Yuri Kharlov, 28 September 2000
 
@@ -590,13 +602,13 @@ void AliPHOSv1::StepManager(void)
 
       // Add the current particle in the list of the CPV hits.
 
-      phos.GetCPVModule(moduleNumber).AddHit(fIshunt,tracknumber,pmom,xyd,ipart);
+      phos.GetCPVModule(moduleNumber).AddHit(fIshunt,primary,pmom,xyd,ipart);
 
-      if (fDebugLevel == 1) {
+      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\n",
-               xyd[0],xyd[1],ipart);
+       printf( "                            xy = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
+               xyd[0],xyd[1],ipart,primary);
       }
 
       // Digitize the current CPV hit:
@@ -668,9 +680,64 @@ void AliPHOSv1::StepManager(void)
     xyze[1] = pos[1] ;
     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[2];
+       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 ) {
       gMC->CurrentVolOffID(10, relid[0]) ; // get the PHOS module number ;
+      if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(5),"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
@@ -678,11 +745,11 @@ void AliPHOSv1::StepManager(void)
       // get the absolute Id number
       
       fGeom->RelToAbsNumbering(relid, absid) ; 
-      
+
       // add current hit to the hit list
       
       AddHit(fIshunt, primary,tracknumber, absid, xyze, trackpid);
-      
+
     } // there is deposited energy
   } // we are inside a PHOS Xtal
 }
@@ -789,8 +856,8 @@ void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, Int_t moduleNumbe
 
   // Finite size of ionization region
 
-  Int_t nCellZ  = fGeom->GetNumberOfPadsZ();
-  Int_t nCellX  = fGeom->GetNumberOfPadsPhi();
+  Int_t nCellZ  = fGeom->GetNumberOfCPVPadsZ();
+  Int_t nCellX  = fGeom->GetNumberOfCPVPadsPhi();
   Int_t nz3     = (kNgamz+1)/2;
   Int_t nx3     = (kNgamx+1)/2;
   cpvDigits->Expand(nIter*kNgamx*kNgamz);
index 9283ff6..31cc779 100644 (file)
@@ -39,7 +39,7 @@ public:
 
   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   FinishEvent(void) ;                               
+  virtual void   Hit2Digit(Int_t event) ;
   virtual Int_t  IsVersion(void) const {
     // Gives the version number 
     return 1 ; 
@@ -47,7 +47,6 @@ public:
   virtual void   MakeBranch(Option_t* opt) ;
   void           Reconstruction(AliPHOSReconstructioner * Reconstructioner) ;
   void           ResetClusters(){} ;
-  virtual void   ResetDigits() ; 
   virtual void   ResetHits() ; 
   virtual void   ResetReconstruction() ; // Reset reconstructed objects
   void           SetReconstructioner(AliPHOSReconstructioner& Reconstructioner) {
@@ -68,9 +67,9 @@ public:
     return *this ; 
   }
 
-  // IHEP's CPV specific functions
-
+  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) ;
   Double_t   CPVCumulPadResponse(Double_t x, Double_t y) ;
@@ -78,13 +77,11 @@ public:
 protected:
 
   Float_t fDigitThreshold ;                       // Threshold for the digit registration 
-  Int_t fNTmpHits ;                               //!  Used internally for digitalization
   Float_t fPinElectronicNoise  ;                  // Electronic Noise in the PIN
-  AliPHOSReconstructioner * fReconstructioner ;   // Reconstrutioner of the PHOS event: Clusterization and subtracking procedures
-  TClonesArray * fTmpHits ;                       //!  Used internally for digitalization 
+  AliPHOSReconstructioner  * fReconstructioner ;  // Reconstrutioner of the PHOS event: Clusterization and subtracking procedures
   AliPHOSTrackSegmentMaker * fTrackSegmentMaker ; // Reconstructioner of the PHOS track segment: 2 x PPSD + 1 x EMC
-
-  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 fa30186..482510b 100644 (file)
 // AddHit, StepManager,and FinishEvent are redefined 
 //                  
 //*-- Author: Gines MARTINEZ (SUBATECH)
-
+//*-- Modified Nov. 22 2000 by Dmitri Peressounko
+// All hits are stored.
+// Note, that primaries will not be assigned to digits:
+// becouse of tiny energy deposition at each step.
+//  
 
 // --- ROOT system ---
 
@@ -49,8 +53,7 @@ ClassImp(AliPHOSv2)
 AliPHOSv2::AliPHOSv2()
 {
   // default ctor
-  fNTmpHits = 0 ; 
-  fTmpHits  = 0 ; 
+
 }
 
 //____________________________________________________________________________
@@ -58,80 +61,26 @@ AliPHOSv2::AliPHOSv2(const char *name, const char *title):
 AliPHOSv1(name,title)
 {
   // ctor
-  
-   fHits= new TClonesArray("AliPHOSHit",1000) ;
 }
 
 //____________________________________________________________________________
 AliPHOSv2::~AliPHOSv2()
 {
   // dtor
-  if ( fTmpHits) {
-    fTmpHits->Delete() ; 
-    delete fTmpHits ;
-    fTmpHits = 0 ; 
-  }
-  
-  if ( fEmcRecPoints ) { 
-    fEmcRecPoints->Delete() ; 
-    delete fEmcRecPoints ; 
-    fEmcRecPoints = 0 ; 
-  }
-  
-  if ( fPpsdRecPoints ) { 
-    fPpsdRecPoints->Delete() ;
-    delete fPpsdRecPoints ;
-    fPpsdRecPoints = 0 ; 
-  }
-  if ( fTrackSegments ) {
-    fTrackSegments->Delete() ; 
-    delete fTrackSegments ;
-    fTrackSegments = 0 ; 
-  }
 }
 
 //____________________________________________________________________________
 void AliPHOSv2::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id, Float_t * hits, Int_t pid)
 {
   // Add a hit to the hit list.
-  // A PHOS hit is the sum of all hits in a single crystal
-  //   or in a single PPSD gas cell
 
-  Int_t hitCounter ;
-  TClonesArray &ltmphits = *fTmpHits ;
   AliPHOSHit *newHit ;
-  AliPHOSHit *curHit ;
-  Bool_t deja = kFALSE ;
-
-  // In any case, fills the fTmpHit TClonesArray (with "accumulated hits")
 
   newHit = new AliPHOSHit(shunt, primary, tracknumber, Id, hits, pid) ;
 
-  // We do want to save in TreeH the raw hits 
-  TClonesArray &lhits = *fHits;
-
-  for ( hitCounter = 0 ; hitCounter < fNTmpHits && !deja ; hitCounter++ ) {
-    curHit = (AliPHOSHit*) ltmphits[hitCounter] ;
-  if( *curHit == *newHit ) {
-    *curHit = *curHit + *newHit ;
-    deja = kTRUE ;
-    }
-  }
-         
-  if ( !deja ) {
-    new(ltmphits[fNTmpHits]) AliPHOSHit(*newHit) ;
-    fNTmpHits++ ;
-  }
-
-  // We do want to save in TreeH the raw hits 
-  new(lhits[fNhits]) AliPHOSHit(*newHit) ;    
+  new((*fHits)[fNhits]) AliPHOSHit(*newHit) ;    
   fNhits++ ;
 
-  // Please note that the fTmpHits array must survive up to the
-  // end of the events, so it does not appear e.g. in ResetHits() (
-  // which is called at the end of each primary).  
-
   delete newHit;
 
 }