]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOSAnalyze.cxx
correct a typo SetTraskSegment to SetTrackSegment
[u/mrichter/AliRoot.git] / PHOS / AliPHOSAnalyze.cxx
index 9c74233fd22713622da6bfc1cfc5c629296bdf09..853ce1b29c0b0bf74f1111caa3f8142abcc056fe 100644 (file)
@@ -89,7 +89,6 @@
 #include "AliPHOSTrackSegment.h"
 #include "AliPHOSRecParticle.h"
 #include "AliPHOSCpvRecPoint.h"
-#include "AliPHOSPpsdRecPoint.h"
 #include "AliPHOSGetter.h"
 
 
@@ -128,24 +127,58 @@ void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod,const char * branchName,c
   //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 * sdigitOccupancy = new TH2F("sdigitOccupancy","EMC sdigits", 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 ObjectGetter
   AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ;
+  if(Nevent >= gAlice->TreeE()->GetEntries() ) {
+    cout << "There is no event " << Nevent << ", only " << gAlice->TreeE()->GetEntries() << "events available " <<endl ;
+    return ;
+  }
   const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ; 
   gime->Event(Nevent);
+
+  Int_t nx = phosgeom->GetNPhi() ;
+  Int_t nz = phosgeom->GetNZ() ;
+  Float_t * cri= phosgeom->GetEMCAGeometry()->GetCrystalHalfSize() ;
+  Float_t x = nx*cri[0] ;
+  Float_t z = nz*cri[2] ;
+  Int_t nxCPV = (Int_t) (nx*phosgeom->GetPadSizePhi()/(2.*cri[0])) ;
+  Int_t nzCPV = (Int_t) (nz*phosgeom->GetPadSizeZ()/(2.*cri[2])) ;
+  
+  TH2F * emcDigits = (TH2F*) gROOT->FindObject("emcDigits") ;
+  if(emcDigits)
+    emcDigits->Delete() ;
+  emcDigits = new TH2F("emcDigits","EMC digits",  nx,-x,x,nz,-z,z);
+  TH2F * emcSdigits =(TH2F*) gROOT->FindObject("emcSdigits") ;  
+  if(emcSdigits)
+    emcSdigits->Delete() ;
+  emcSdigits = new TH2F("emcSdigits","EMC sdigits", nx,-x,x,nz,-z,z);
+  TH2F * emcRecPoints = (TH2F*)gROOT->FindObject("emcRecPoints") ; 
+  if(emcRecPoints)
+    emcRecPoints->Delete() ;
+  emcRecPoints = new TH2F("emcRecPoints","EMC RecPoints",nx,-x,x,nz,-z,z);
+  TH2F * cpvSdigits =(TH2F*) gROOT->FindObject("cpvSdigits") ;
+  if(cpvSdigits)
+    cpvSdigits->Delete() ;
+  cpvSdigits = new TH2F("cpvSdigits","CPV sdigits", nx,-x,x,nz,-z,z);
+  TH2F * cpvDigits = (TH2F*)gROOT->FindObject("cpvDigits") ;
+  if(cpvDigits)
+    cpvDigits->Delete() ;
+  cpvDigits = new TH2F("cpvDigits","CPV digits",   nxCPV,-x,x,nzCPV,-z,z) ;
+  TH2F * cpvRecPoints= (TH2F*)gROOT->FindObject("cpvRecPoints") ; 
+  if(cpvRecPoints)
+    cpvRecPoints->Delete() ;
+  cpvRecPoints = new TH2F("cpvRecPoints","CPV RecPoints",    nxCPV,-x,x,nzCPV,-z,z) ;
+
+  TH2F * phot = (TH2F*)gROOT->FindObject("phot") ;
+  if(phot)
+    phot->Delete() ;
+  phot = new TH2F("phot","Primary Photon",  nx,-x,x,nz,-z,z);
+  TH2F * recPhot = (TH2F*)gROOT->FindObject("recPhot") ; 
+  if(recPhot)
+    recPhot->Delete() ;
+  recPhot = new TH2F("recPhot","RecParticles with primary Photon",nx,-x,x,nz,-z,z);
+  
   
   //Plot Primary Particles
   const TParticle * primary ;
@@ -154,13 +187,14 @@ void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod,const char * branchName,c
     {
       primary = gime->Primary(iPrimary) ;
       Int_t primaryType = primary->GetPdgCode() ;
-      if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212) ) {
-        Int_t moduleNumber ;
-        Double_t primX, primZ ;
-        phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
-        if(moduleNumber==Nmod)
-          charg->Fill(primZ,primX,primary->Energy()) ;
-      }
+//       if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212)
+//       ||(primaryType == 11)||(primaryType == -11) ) {
+//         Int_t moduleNumber ;
+//         Double_t primX, primZ ;
+//         phosgeom->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 ;
@@ -168,124 +202,114 @@ void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod,const char * branchName,c
         if(moduleNumber==Nmod) 
           phot->Fill(primZ,primX,primary->Energy()) ;
       }
-      else{
-        if( primaryType == -2112 ) {
-          Int_t moduleNumber ;
-          Double_t primX, primZ ;
-          phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
-          if(moduleNumber==Nmod)
-            nbar->Fill(primZ,primX,primary->Energy()) ;
-        }
-      }
+//       else{
+//         if( primaryType == -2112 ) {
+//           Int_t moduleNumber ;
+//           Double_t primX, primZ ;
+//           phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
+//           if(moduleNumber==Nmod)
+//             nbar->Fill(primZ,primX,primary->Energy()) ;
+//         }
+//       }
     }  
+
   
   Int_t iSDigit ;
   AliPHOSDigit * sdigit ;
-  
-  TClonesArray * sdigits = gime->SDigits(branchTitle,ffileName) ;
-  if(sdigits)
-    for(iSDigit = 0; iSDigit < sdigits->GetEntriesFast(); iSDigit++)
+  const TClonesArray * sdigits = gime->SDigits() ;
+  Int_t nsdig[5] = {0,0,0,0,0} ;
+  if(sdigits){
+    for(iSDigit = 0; iSDigit < sdigits->GetEntriesFast() ; iSDigit++)
       {
-       sdigit = (AliPHOSDigit*)sdigits->At(iSDigit) ;
+       sdigit = (AliPHOSDigit *) sdigits->At(iSDigit) ;
        Int_t relid[4];
        phosgeom->AbsToRelNumbering(sdigit->GetId(), relid) ;
        Float_t x,z ;
        phosgeom->RelPosInModule(relid,x,z) ;
-       Float_t e = gime->SDigitizer(branchTitle)->Calibrate(sdigit->GetAmp()) ;
+       Float_t e = gime->SDigitizer()->Calibrate(sdigit->GetAmp()) ;
+       nsdig[relid[0]-1]++ ;
        if(relid[0]==Nmod){
          if(relid[1]==0)  //EMC
-           sdigitOccupancy->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) ;
+           emcSdigits->Fill(x,z,e) ;
+         if( relid[1]!=0 )
+           cpvSdigits->Fill(x,z,e) ;
        }
       }
-  
+  }
+  cout << "Number of EMC + CPV SDigits per module: " <<endl ;
+  cout << nsdig[0] << " " << nsdig[1] << " " << nsdig[2] << " " << nsdig[3]<< " " << nsdig[4] << endl ;
+  cout << endl ;
+
+
   //Plot digits
   Int_t iDigit ;
   AliPHOSDigit * digit ;
-  TClonesArray * digits = gime->Digits(branchTitle) ;
-  if(digits)
+  const TClonesArray * digits = gime->Digits(); 
+  if(digits) {
     for(iDigit = 0; iDigit < digits->GetEntriesFast(); iDigit++)
       {
-       digit = (AliPHOSDigit*) digits->At(iDigit) ;
+       digit = (AliPHOSDigit *) digits->At(iDigit) ;
        Int_t relid[4];
        phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
        Float_t x,z ;
        phosgeom->RelPosInModule(relid,x,z) ;
-       Float_t e = gime->SDigitizer(branchTitle)->Calibrate(digit->GetAmp()) ;
+       Float_t e = gime->SDigitizer()->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) ;
+           emcDigits->Fill(x,z,e) ;
+         if( relid[1]!=0 )
+           cpvDigits->Fill(x,z,e) ;
        }
       }
+  }
   
   
   //Plot RecPoints
   Int_t irecp ;
   TVector3 pos ;
-  TObjArray * emcrp = gime->EmcRecPoints(branchTitle) ;
-  TObjArray * cpvrp = gime->CpvRecPoints(branchTitle) ;
-  if(cpvrp)
+  TObjArray * emcrp = gime->EmcRecPoints() ;
+  if(emcrp) {
     for(irecp = 0; irecp < emcrp->GetEntriesFast() ; irecp ++){
-      const AliPHOSEmcRecPoint * emc= (AliPHOSEmcRecPoint *) emcrp->At(irecp) ;
+      AliPHOSEmcRecPoint * emc = (AliPHOSEmcRecPoint *) emcrp->At(irecp) ;
       if(emc->GetPHOSMod()==Nmod){
        emc->GetLocalPosition(pos) ;
-       emcOccupancy->Fill(pos.X(),pos.Z(),emc->GetEnergy());
+       emcRecPoints->Fill(pos.X(),pos.Z(),emc->GetEnergy());
       }
     }
+  }
   
-  if(cpvrp)
+  TObjArray * cpvrp = gime->CpvRecPoints() ;
+  if(cpvrp) {
     for(irecp = 0; irecp < cpvrp->GetEntriesFast() ; irecp ++){
-      const AliPHOSRecPoint * cpv = (AliPHOSRecPoint *) cpvrp->At(irecp) ;
-      if((strcmp(cpv->ClassName(),"AliPHOSPpsdRecPoint" )) == 0){  // PPSD Rec Point
-       AliPHOSPpsdRecPoint * ppsd = (AliPHOSPpsdRecPoint*) cpv ;
-       ppsd->GetLocalPosition(pos) ;
-       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());
-       }
-      }
-      else{
-       AliPHOSCpvRecPoint * cpv1 = (AliPHOSCpvRecPoint*) cpv ;
-       if(cpv1->GetPHOSMod()==Nmod){
-         cpv1->GetLocalPosition(pos) ;
-         ppsdUpCl->Fill(pos.X(),pos.Z(),cpv1->GetEnergy());
-       }
+      AliPHOSRecPoint * cpv = (AliPHOSCpvRecPoint *) cpvrp->At(irecp) ;
+      if(cpv->GetPHOSMod()==Nmod){
+       cpv->GetLocalPosition(pos) ;
+       cpvRecPoints->Fill(pos.X(),pos.Z(),cpv->GetEnergy());
       }
     }
-  
-  
+  }
+    
   //Plot RecParticles
-  const AliPHOSRecParticle * recParticle ;
+  AliPHOSRecParticle * recParticle ;
   Int_t iRecParticle ;
-  TClonesArray * rps = gime->RecParticles(branchTitle) ;
-  TClonesArray * tss = gime->TrackSegments(branchTitle) ;
-  if(rps && tss && emcrp)
-    for(iRecParticle = 0; iRecParticle < rps->GetEntriesFast() ;iRecParticle++ )
+  TClonesArray * rp = gime->RecParticles() ;
+  TClonesArray * ts = gime->TrackSegments() ;
+  if(rp && ts && emcrp) {
+    for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ; iRecParticle++ )
       {
-       recParticle =  (AliPHOSRecParticle *) rps->At(iRecParticle) ;
+       recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
        Int_t moduleNumberRec ;
        Double_t recX, recZ ;
        phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
        if(moduleNumberRec == Nmod){
          
          Double_t minDistance = 5. ;
-         Int_t closestPrimary = -1 ;
-         
-         
+         Int_t closestPrimary = -1 ;   
+
          //extract list of primaries: it is stored at EMC RecPoints
-         Int_t emcIndex = ((AliPHOSTrackSegment *) tss->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
+         Int_t emcIndex = ((AliPHOSTrackSegment *) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
          Int_t numberofprimaries ;
-         Int_t * listofprimaries  = ((AliPHOSEmcRecPoint *)emcrp->At(emcIndex))->GetPrimaries(numberofprimaries)  ;
+         Int_t * listofprimaries  = ((AliPHOSRecPoint*) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries)  ;
          Int_t index ;
          const TParticle * primary ;
          Double_t distance = minDistance ;
@@ -310,28 +334,23 @@ void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod,const char * branchName,c
            
            if(primaryType==22)
              recPhot->Fill(recZ,recX,recParticle->Energy()) ;
-           else
-             if(primaryType==-2112)
-               recNbar->Fill(recZ,recX,recParticle->Energy()) ; 
+//         else
+//           if(primaryType==-2112)
+//             recNbar->Fill(recZ,recX,recParticle->Energy()) ; 
          }
        }
       }
-  
+
+  }
   
   //Plot made histograms
-  digitOccupancy->Draw("box") ;
-  sdigitOccupancy->SetLineColor(5) ;
-  sdigitOccupancy->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") ;
+  emcSdigits->Draw("box") ;
+  emcDigits->SetLineColor(5) ;
+  emcDigits->Draw("boxsame") ;
+  emcRecPoints->SetLineColor(2) ;
+  emcRecPoints->Draw("boxsame") ;
+  cpvSdigits->SetLineColor(1) ;
+  cpvSdigits->Draw("boxsame") ;
   
 }
 //____________________________________________________________________________
@@ -422,29 +441,33 @@ void AliPHOSAnalyze::Ls(){
   Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; 
   //  for(event = 0; event < gime->MaxEvent(); event++  ){
   for(event = 0; event < maxevent; event++  ){
-    gime->Event(event);
+    gime->Event(event,"R");  //will read only TreeR 
     
     //copy EM RecParticles to the "total" list        
     const AliPHOSRecParticle * recParticle ;
     Int_t iRecParticle ;
-    TClonesArray * rp = gime->RecParticles(branchTitle) ;
-    if(rp) 
-      for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast()  ;iRecParticle++ )
-       {
-         recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
-         if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
-            (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM))
-           new( (*allRecParticleList)[iRecPhot++] ) AliPHOSRecParticle(*recParticle) ;
-       }
+    TClonesArray * rp = gime->RecParticles() ;
+    if(!rp){
+      cout << "AliPHOSAnalyze::InvariantMass --> Can't find RecParticles " << endl ;
+      return ;
+    }
+
+    for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ )
+      {
+       recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
+       if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)||
+          (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW))
+         new( (*allRecParticleList)[iRecPhot++] ) AliPHOSRecParticle(*recParticle) ;
+      }
     
     Int_t mevent = event%nMixedEvents ; //event number in the "mixed" cicle
     nRecParticles[mevent] = iRecPhot-1 ;  
-
+    
     //check, if it is time to calculate invariant mass?
     Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; 
     if((mevent == 0) && (event +1 == maxevent)){
-
-    //   if((mevent == 0) && (event +1 == gime->MaxEvent())){
+      
+      //   if((mevent == 0) && (event +1 == gime->MaxEvent())){
       
       //calculate invariant mass:
       Int_t irp1,irp2 ;
@@ -474,14 +497,14 @@ void AliPHOSAnalyze::Ls(){
          
          if(irp2 <= nRecParticles[nCurEvent]){ //'Real' event
            hRealEM->Fill(invMass,pt);
-           if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&
-              (rp2->GetType() == AliPHOSFastRecParticle::kGAMMA) )
+           if((rp1->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)&&
+              (rp2->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST) )
              hRealPhot->Fill(invMass,pt);
          }
          else{
            hMixedEM->Fill(invMass,pt);
-           if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&
-              (rp2->GetType() == AliPHOSFastRecParticle::kGAMMA) )
+           if((rp1->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)&&
+              (rp2->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST) )
              hMixedPhot->Fill(invMass,pt);
          } //real-mixed
          
@@ -543,7 +566,6 @@ void AliPHOSAnalyze::Ls(){
 
   Int_t ievent;
   Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; 
-  //  for ( ievent=0; ievent < gime->MaxEvent() ; ievent++){
   for ( ievent=0; ievent < maxevent ; ievent++){
 
     //read the current event
@@ -551,67 +573,82 @@ void AliPHOSAnalyze::Ls(){
 
     const AliPHOSRecParticle * recParticle ;
     Int_t iRecParticle ;
-    TClonesArray * rp = gime->RecParticles(branchTitle) ;
-    TClonesArray * tss = gime->TrackSegments(branchTitle) ;
-    TObjArray * emcrp = gime->EmcRecPoints(branchTitle) ;
-    if(rp && tss && emcrp)
-      for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ; iRecParticle++ ){
-       recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
-       
-       //find the closest primary
-       Int_t moduleNumberRec ;
-       Double_t recX, recZ ;
-       phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
-       
-       Double_t minDistance  = 100. ;
-       Int_t closestPrimary = -1 ;
-       
-       //extract list of primaries: it is stored at EMC RecPoints
-       Int_t emcIndex = ((AliPHOSTrackSegment *)tss->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
-       Int_t numberofprimaries ;
-       Int_t * listofprimaries  = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries)  ;
-       
-       Int_t index ;
-       const TParticle * primary ;
-       Double_t distance = minDistance ;
-       Double_t dX, dZ; 
-       Double_t dXmin = 0.; 
-       Double_t dZmin = 0. ;
-       for ( index = 0 ; index < numberofprimaries ; index++){
-         primary = gime->Primary(listofprimaries[index]) ;
-         Int_t moduleNumber ;
-         Double_t primX, primZ ;
-         phosgeom->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] ;
-           }
+    TClonesArray * rp = gime->RecParticles() ;
+    if(!rp) {
+      cout << "AliPHOSAnalyze::EnergyResolution --> Event " <<ievent 
+          << ",  Can't find RecParticles " << endl ;
+      return ;
+    }
+    TClonesArray * ts = gime->TrackSegments() ;
+    if(!ts) {
+      cout << "AliPHOSAnalyze::EnergyResolution --> Event " <<ievent 
+          << ",  Can't find TrackSegments " << endl ;
+      return ;
+    }
+    TObjArray * emcrp = gime->EmcRecPoints() ;
+    if(!emcrp){
+      cout << "AliPHOSAnalyze::EnergyResolution --> Event " <<ievent 
+          << ",  Can't find EmcRecPoints " << endl ;
+      return ;
+    }
+      
+    for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ;iRecParticle++ ){
+      recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
+      
+      //find the closest primary
+      Int_t moduleNumberRec ;
+      Double_t recX, recZ ;
+      phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
+      
+      Double_t minDistance  = 100. ;
+      Int_t closestPrimary = -1 ;
+      
+      //extract list of primaries: it is stored at EMC RecPoints
+      Int_t emcIndex = ((AliPHOSTrackSegment*) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
+      Int_t numberofprimaries ;
+      Int_t * listofprimaries  = ((AliPHOSEmcRecPoint*) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries)  ;
+      
+      Int_t index ;
+      const TParticle * primary ;
+      Double_t distance = minDistance ;
+      Double_t dX, dZ; 
+      Double_t dXmin = 0.; 
+      Double_t dZmin = 0. ;
+      for ( index = 0 ; index < numberofprimaries ; index++){
+       primary = gime->Primary(listofprimaries[index]) ;
+       Int_t moduleNumber ;
+       Double_t primX, primZ ;
+       phosgeom->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] ;
          }
        }
-       
-       //if found primary, fill histograms
-       if(closestPrimary >=0 ){
-         const TParticle * primary = gime->Primary(closestPrimary) ;
-         if(primary->GetPdgCode() == 22){
-           hAllEnergy->Fill(primary->Energy(), recParticle->Energy()) ;
-           if(recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA){
-             hPhotEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; 
-             hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; 
-           }
-           else
-             if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)
-               hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; 
+      }
+
+      //if found primary, fill histograms
+      if(closestPrimary >=0 ){
+       const TParticle * primary = gime->Primary(closestPrimary) ;
+       if(primary->GetPdgCode() == 22){
+         hAllEnergy->Fill(primary->Energy(), recParticle->Energy()) ;
+         if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){
+           hPhotEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; 
+           hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; 
          }
+         else
+           if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)
+             hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; 
        }
       }
+    }
   }
-  
+
   //write filled histograms
   efile->cd() ;
   hAllEnergy->Write(0,kOverwrite) ;
@@ -620,7 +657,7 @@ void AliPHOSAnalyze::Ls(){
   //  efile->Write() ;
   efile->Close() ;
   delete efile ;
-  
+
 }
 //____________________________________________________________________________
 void AliPHOSAnalyze::PositionResolution(const char * branchTitle)    
@@ -666,77 +703,91 @@ void AliPHOSAnalyze::PositionResolution(const char * branchTitle)
   const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ; 
 
   Int_t ievent;
- Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; 
-  //  for ( ievent=0; ievent < gime->MaxEvent() ; ievent++){
+  Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; 
   for ( ievent=0; ievent < maxevent ; ievent++){
-
+    
     //read the current event
     gime->Event(ievent) ;
+    TClonesArray * rp = gime->RecParticles() ;
+    if(!rp) {
+      cout << "AliPHOSAnalyze::PositionResolution --> Event " <<ievent 
+          << ",  Can't find RecParticles " << endl ;
+      return ;
+    }
+    TClonesArray * ts = gime->TrackSegments() ;
+    if(!ts) {
+      cout << "AliPHOSAnalyze::PositionResolution --> Event " <<ievent 
+          << ",  Can't find TrackSegments " << endl ;
+      return ;
+    }
+    TObjArray * emcrp = gime->EmcRecPoints() ;
+    if(!emcrp){
+      cout << "AliPHOSAnalyze::PositionResolution --> Event " <<ievent 
+          << ",  Can't find EmcRecPoints " << endl ;
+      return ;
+    }
     
     const AliPHOSRecParticle * recParticle ;
     Int_t iRecParticle ;
-    TClonesArray * rp = gime->RecParticles(branchTitle) ;
-    TClonesArray * tss= gime->TrackSegments(branchTitle) ;
-    TObjArray * emcrp = gime->EmcRecPoints(branchTitle) ;
-    if( rp && tss && emcrp )
-      for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ;iRecParticle++ ){
-       recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
-       
-       //find the closest primary
-       Int_t moduleNumberRec ;
-       Double_t recX, recZ ;
-       phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
-       
-       Double_t minDistance  = 100. ;
-       Int_t closestPrimary = -1 ;
-       
-       //extract list of primaries: it is stored at EMC RecPoints
-       Int_t emcIndex = ((AliPHOSTrackSegment *) tss->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
-       Int_t numberofprimaries ;
-       Int_t * listofprimaries  = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries)  ;
-       
-       Int_t index ;
-       const TParticle * primary ;
-       Double_t distance = minDistance ;
-       Double_t dX = 1000; // incredible number
-       Double_t dZ = 1000; // for the case if no primary will be found
-       Double_t dXmin = 0.; 
-       Double_t dZmin = 0. ;
-       for ( index = 0 ; index < numberofprimaries ; index++){
-         primary = gime->Primary(listofprimaries[index]) ;
-         Int_t moduleNumber ;
-         Double_t primX, primZ ;
-         phosgeom->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] ;
-           }
+    for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ ){
+      recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
+      
+      //find the closest primary
+      Int_t moduleNumberRec ;
+      Double_t recX, recZ ;
+      phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
+      
+      Double_t minDistance  = 100. ;
+      Int_t closestPrimary = -1 ;
+      
+      //extract list of primaries: it is stored at EMC RecPoints
+      Int_t emcIndex = ((AliPHOSTrackSegment*) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
+      Int_t numberofprimaries ;
+      Int_t * listofprimaries  = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries)  ;
+
+      Int_t index ;
+      const TParticle * primary ;
+      Double_t distance = minDistance ;
+      Double_t dX = 1000; // incredible number
+      Double_t dZ = 1000; // for the case if no primary will be found
+      Double_t dXmin = 0.; 
+      Double_t dZmin = 0. ;
+      for ( index = 0 ; index < numberofprimaries ; index++){
+       primary = gime->Primary(listofprimaries[index]) ;
+       Int_t moduleNumber ;
+       Double_t primX, primZ ;
+       phosgeom->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] ;
          }
        }
-       
-       //if found primary, fill histograms
-       if(closestPrimary >=0 ){
-         const TParticle * primary = gime->Primary(closestPrimary) ;
-         if(primary->GetPdgCode() == 22){
-           hAllPosition->Fill(primary->Energy(), minDistance) ;
-           hAllPositionX->Fill(primary->Energy(), dX) ;
-           hAllPositionZ->Fill(primary->Energy(), dZ) ;
-           if(recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA){
-             hPhotPosition->Fill(primary->Energy(), minDistance ) ; 
-             hEMPosition->Fill(primary->Energy(), minDistance ) ; 
-           }
-           else
-             if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)
-               hEMPosition->Fill(primary->Energy(), minDistance ) ; 
+      }
+      
+      //if found primary, fill histograms
+      if(closestPrimary >=0 ){
+       const TParticle * primary = gime->Primary(closestPrimary) ;
+       if(primary->GetPdgCode() == 22){
+         hAllPosition->Fill(primary->Energy(), minDistance) ;
+         hAllPositionX->Fill(primary->Energy(), dX) ;
+         hAllPositionZ->Fill(primary->Energy(), dZ) ;
+         if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){
+           hPhotPosition->Fill(primary->Energy(), minDistance ) ; 
+           hEMPosition->Fill(primary->Energy(), minDistance ) ; 
          }
+         else
+           if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)
+             hEMPosition->Fill(primary->Energy(), minDistance ) ; 
        }
       }
+    }
   }
   
   //Write output histgrams
@@ -749,8 +800,8 @@ void AliPHOSAnalyze::PositionResolution(const char * branchTitle)
   pfile->Write() ;
   pfile->Close() ;
   delete pfile ;
-  
-  
+
+
 }
 //____________________________________________________________________________
 void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){
@@ -762,7 +813,6 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){
   TH1F * hPrimary = 0;   //spectrum (P_t distribution) of primary photons         
   TH1F * hAllRP   = 0;   //spectrum of all RecParticles in PHOS
   TH1F * hPhot    = 0;   //spectrum of kGAMMA RecParticles
-  TH1F * hPPSD    = 0;   //spectrum of all RecParticles with PPSD signal
   TH1F * hShape   = 0;   //spectrum of all EM RecParticles
   TH1F * hVeto    = 0;   //spectrum of all neutral RecParticles
 
@@ -770,27 +820,22 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){
   //primary - photon
   TH1F * hPhotReg = 0;   //Registeres as photon
   TH1F * hPhotEM  = 0;   //Registered as EM       
-  TH1F * hPhotPPSD= 0;   //Registered as RecParticle with PPSD signal        
 
   //primary - n
   TH1F * hNReg = 0;   //Registeres as photon          
   TH1F * hNEM  = 0;   //Registered as EM            
-  TH1F * hNPPSD= 0;   //Registered as RecParticle with PPSD signal           
 
   //primary - nBar
   TH1F * hNBarReg = 0;   //Registeres as photon
   TH1F * hNBarEM  = 0;   //Registered as EM          
-  TH1F * hNBarPPSD= 0;   //Registered as RecParticle with PPSD signal             
 
   //primary - charged hadron (pBar excluded)
   TH1F * hChargedReg = 0;   //Registeres as photon  
   TH1F * hChargedEM  = 0;   //Registered as EM           
-  TH1F * hChargedPPSD= 0;   //Registered as RecParticle with PPSD signal           
 
   //primary - pBar
   TH1F * hPbarReg = 0;   //Registeres as photon  
   TH1F * hPbarEM  = 0;   //Registered as EM 
-  TH1F * hPbarPPSD= 0;   //Registered as RecParticle with PPSD signal   
 
 
   //Reading histograms from the file
@@ -806,9 +851,6 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){
   hPhot  = (TH1F*)cfile->Get("hPhot") ;
   if(hPhot == 0)
     hPhot = new TH1F("hPhot","All kGAMMA RecParticles",100, 0., 5.);
-  hPPSD = (TH1F*)cfile->Get("hPPSD") ;
-  if(hPPSD == 0)
-    hPPSD  = new TH1F("hPPSD", "All PPSD Recparticles",    100, 0., 5.);
   hShape = (TH1F*) cfile->Get("hShape") ;
   if(hShape == 0)
     hShape = new TH1F("hShape","All particles with EM shower",100, 0., 5.);
@@ -824,9 +866,6 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){
   hPhotEM  =(TH1F*)cfile->Get("hPhotEM");
   if(hPhotEM== 0)
     hPhotEM   = new TH1F("hPhotEM",  "Photon registered as EM", 100, 0., 5.);
-  hPhotPPSD= (TH1F*)cfile->Get("hPhotPPSD");
-  if(hPhotPPSD== 0)
-    hPhotPPSD   = new TH1F("hPhotPPSD","Photon registered as PPSD", 100, 0., 5.);
 
   //primary - n
   hNReg = (TH1F*)cfile->Get("hNReg");
@@ -835,9 +874,6 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){
   hNEM  = (TH1F*)cfile->Get("hNEM"); 
   if(hNEM== 0)
     hNEM      = new TH1F("hNEM",    "N registered as EM",      100, 0., 5.);
-  hNPPSD=(TH1F*)cfile->Get("hNPPSD"); 
-  if(hNPPSD== 0)
-    hNPPSD      = new TH1F("hNPPSD","N registered as PPSD",      100, 0., 5.);
 
   //primary - nBar
   hNBarReg =(TH1F*)cfile->Get("hNBarReg");
@@ -846,9 +882,6 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){
   hNBarEM  =(TH1F*)cfile->Get("hNBarEM"); 
   if(hNBarEM== 0)
     hNBarEM   = new TH1F("hNBarEM",  "NBar registered as EM",   100, 0., 5.);
-  hNBarPPSD=(TH1F*)cfile->Get("hNBarPPSD");
-  if(hNBarPPSD== 0)
-    hNBarPPSD   = new TH1F("hNBarPPSD","NBar registered as PPSD",   100, 0., 5.);
 
   //primary - charged hadron (pBar excluded)
   hChargedReg = (TH1F*)cfile->Get("hChargedReg");
@@ -857,10 +890,7 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){
   hChargedEM  = (TH1F*)cfile->Get("hChargedEM"); 
   if(hChargedEM== 0)
     hChargedEM= new TH1F("hChargedEM","Charged registered as EM",100, 0., 5.);
-  hChargedPPSD= (TH1F*)cfile->Get("hChargedPPSD");
-  if(hChargedPPSD== 0)
-    hChargedPPSD= new TH1F("hChargedPPSD","Charged registered as PPSD",100, 0., 5.);
-  
   //primary - pBar
   hPbarReg = (TH1F*)cfile->Get("hPbarReg");
   if(hPbarReg== 0)
@@ -868,9 +898,6 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){
   hPbarEM  = (TH1F*)cfile->Get("hPbarEM");
   if(hPbarEM== 0)
     hPbarEM= new TH1F("hPbarEM","Pbar registered as EM",100, 0., 5.);
-  hPbarPPSD= (TH1F*)cfile->Get("hPbarPPSD");
-  if(hPbarPPSD== 0)
-    hPbarPPSD= new TH1F("hPbarPPSD","Pbar as PPSD",100, 0., 5.);
   
 
   //Now make some initializations
@@ -888,11 +915,30 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){
   
   Int_t ievent;
   Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; 
-  //  for ( ievent=0; ievent < gime->MaxEvent() ; ievent++){
   for ( ievent=0; ievent < maxevent ; ievent++){
     
     gime->Event(ievent) ;
     
+    TClonesArray * rp = gime->RecParticles() ;
+    if(!rp) {
+      cout << "AliPHOSAnalyze::Contamination --> Event " <<ievent 
+          << ",  Can't find RecParticles " << endl ;
+      return ;
+    }
+    TClonesArray * ts = gime->TrackSegments() ;
+    if(!ts) {
+      cout << "AliPHOSAnalyze::Contamination --> Event " <<ievent 
+          << ",  Can't find TrackSegments " << endl ;
+      return ;
+    }
+    TObjArray * emcrp = gime->EmcRecPoints() ;
+    if(!emcrp){
+      cout << "AliPHOSAnalyze::Contamination --> Event " <<ievent 
+          << ",  Can't find EmcRecPoints " << endl ;
+      return ;
+    }
+    
+    
     //=========== Make spectrum of the primary photons
     const TParticle * primary ;
     Int_t iPrimary ;
@@ -910,162 +956,136 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){
       }
       
     }
+    
     //========== Now scan over RecParticles            
     const AliPHOSRecParticle * recParticle ;
     Int_t iRecParticle ;
-    TClonesArray * rp = gime->RecParticles(RecPointsTitle) ;
-    TClonesArray * tss= gime->TrackSegments(RecPointsTitle) ;
-    TObjArray * emcrp = gime->EmcRecPoints(RecPointsTitle) ;
-    if( rp && tss && emcrp )      
-      for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ ){
-       recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
-       //fill histo spectrum of all RecParticles
-       hAllRP->Fill(CorrectedEnergy(recParticle->Energy())) ;
-       
-       //==========find the closest primary    
-       Int_t moduleNumberRec ;
-       Double_t recX, recZ ;
-       phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
-       
-       Double_t minDistance  = 100. ;
-       Int_t closestPrimary = -1 ;
-       
-       //extract list of primaries: it is stored at EMC RecPoints
-       Int_t emcIndex = ((AliPHOSTrackSegment *) tss->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
-       Int_t numberofprimaries ;
-       Int_t * listofprimaries  = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries)  ;
-       Int_t index ;
-       const TParticle * primary ;
-       Double_t distance = minDistance ;
-       Double_t dX, dZ; 
-       Double_t dXmin = 0.; 
-       Double_t dZmin = 0. ;
-       for ( index = 0 ; index < numberofprimaries ; index++){
-         primary = gime->Primary(listofprimaries[index]) ;
-         Int_t moduleNumber ;
-         Double_t primX, primZ ;
-         phosgeom->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] ;
-           }
+    for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ ){
+      recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
+      //fill histo spectrum of all RecParticles
+      hAllRP->Fill(CorrectedEnergy(recParticle->Energy())) ;
+      
+      //==========find the closest primary     
+      Int_t moduleNumberRec ;
+      Double_t recX, recZ ;
+      phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
+      
+      Double_t minDistance  = 100. ;
+      Int_t closestPrimary = -1 ;
+      
+      //extract list of primaries: it is stored at EMC RecPoints
+      Int_t emcIndex = ((AliPHOSTrackSegment *) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
+      Int_t numberofprimaries ;
+      Int_t * listofprimaries  = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries)  ;
+      Int_t index ;
+      const TParticle * primary ;
+      Double_t distance = minDistance ;
+      Double_t dX, dZ; 
+      Double_t dXmin = 0.; 
+      Double_t dZmin = 0. ;
+      for ( index = 0 ; index < numberofprimaries ; index++){
+       primary = gime->Primary(listofprimaries[index]) ;
+       Int_t moduleNumber ;
+       Double_t primX, primZ ;
+       phosgeom->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] ;
          }
        }
-       
-       //===========define the "type" of closest primary
-       if(closestPrimary >=0 ){
-         Int_t primaryCode = -1;
-         const TParticle * primary = gime->Primary(closestPrimary) ;
-         Int_t primaryType = primary->GetPdgCode() ;
-         if(primaryType == 22) // photon ?
-           primaryCode = 0 ;
+      }
+      
+      //===========define the "type" of closest primary
+      if(closestPrimary >=0 ){
+       Int_t primaryCode = -1;
+       const TParticle * primary = gime->Primary(closestPrimary) ;
+       Int_t primaryType = primary->GetPdgCode() ;
+       if(primaryType == 22) // photon ?
+         primaryCode = 0 ;
+       else
+         if(primaryType == 2112) // neutron
+           primaryCode = 1 ; 
          else
-           if(primaryType == 2112) // neutron
-             primaryCode = 1 ; 
+           if(primaryType == -2112) // Anti neutron
+             primaryCode = 2 ;
            else
-             if(primaryType == -2112) // Anti neutron
-               primaryCode = 2 ;
-             else
-               if(primaryType == -2122) //Anti proton
-                 primaryCode = 4 ;
-               else {
-                 TParticle tempo(*primary) ; 
-                 if(tempo.GetPDG()->Charge())
-                   primaryCode = 3 ;
-               }
-         
-         //==========Now look at the type of RecParticle
-         Float_t energy = CorrectedEnergy(recParticle->Energy()) ;
-         if(recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA){
-           hPhot->Fill(energy ) ;      
-           switch(primaryCode){
-           case 0:
-             hPhotReg->Fill(energy ) ; 
-             break ;
-           case 1:
-             hNReg->Fill(energy ) ; 
-             break ;
-           case 2:
-             hNBarReg->Fill(energy ) ; 
-             break ;
-           case 3:
-             hChargedReg->Fill(energy ) ;
-             break ;
-           case 4:
-             hPbarReg->Fill(energy ) ;
-             break ;
-           default:
-             break ;
-           }
-         }
-         if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
-            (recParticle->GetType() == AliPHOSFastRecParticle::kGAMMAHA)){ //with PPSD signal
-           hPPSD->Fill(energy ) ; 
-           switch(primaryCode){
-           case 0:
-             hPhotPPSD->Fill(energy ) ; 
-             break ;
-           case 1:
-             hNPPSD->Fill(energy ) ; 
-             break ;
-           case 2:
-             hNBarPPSD->Fill(energy ) ;
-             break ;
-           case 3:
-             hChargedPPSD->Fill(energy ) ;
-             break ;
-           case 4:
-             hPbarPPSD->Fill(energy ) ;
-             break ;
-           default:
-             break ;
-           }
+             if(primaryType == -2122) //Anti proton
+               primaryCode = 4 ;
+             else {
+               TParticle tempo(*primary) ; 
+               if(tempo.GetPDG()->Charge())
+                 primaryCode = 3 ;
+             }
+
+       //==========Now look at the type of RecParticle
+       Float_t energy = CorrectedEnergy(recParticle->Energy()) ;
+       if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){
+         hPhot->Fill(energy ) ;        
+         switch(primaryCode){
+         case 0:
+           hPhotReg->Fill(energy ) ; 
+           break ;
+         case 1:
+           hNReg->Fill(energy ) ; 
+           break ;
+         case 2:
+           hNBarReg->Fill(energy ) ; 
+           break ;
+         case 3:
+           hChargedReg->Fill(energy ) ;
+           break ;
+         case 4:
+           hPbarReg->Fill(energy ) ;
+           break ;
+         default:
+           break ;
          }
-         if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
-            (recParticle->GetType() == AliPHOSFastRecParticle::kELECTRON)||
-            (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)||
-            (recParticle->GetType() == AliPHOSFastRecParticle::kABSURDEM) ){ //with EM shower
-           hShape->Fill(energy ) ;
-           switch(primaryCode){
-           case 0:
-             hPhotEM->Fill(energy ) ; 
-             break ;
-           case 1:
-             hNEM->Fill(energy ) ; 
-             break ;
-           case 2:
-             hNBarEM->Fill(energy ) ; 
-             break ;
-           case 3:
-             hChargedEM->Fill(energy ) ; 
-             break ;
-           case 4:
-             hPbarEM->Fill(energy ) ; 
-             break ;
-           default:
-             break ;
-           }
+       }
+       if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)||
+          (recParticle->GetType() == AliPHOSFastRecParticle::kCHARGEDEMFAST)||
+          (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)||
+          (recParticle->GetType() == AliPHOSFastRecParticle::kCHARGEDEMSLOW) ){ //with EM shower
+         hShape->Fill(energy ) ;
+         switch(primaryCode){
+         case 0:
+           hPhotEM->Fill(energy ) ; 
+           break ;
+         case 1:
+           hNEM->Fill(energy ) ; 
+           break ;
+         case 2:
+           hNBarEM->Fill(energy ) ; 
+           break ;
+         case 3:
+           hChargedEM->Fill(energy ) ; 
+           break ;
+         case 4:
+           hPbarEM->Fill(energy ) ; 
+           break ;
+         default:
+           break ;
          }
-         
-         if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
-            (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHA) ||
-            (recParticle->GetType() == AliPHOSFastRecParticle::kGAMMAHA) ||
-            (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM) ) //nuetral
-           hVeto->Fill(energy ) ;
-         
-         //fill number of primaries identified as ...
-         if(primaryCode >= 0) // Primary code defined
-           counter[recParticle->GetType()][primaryCode]++ ; 
-         
        }
        
-      } // no closest primary found
+       if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)||
+          (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHAFAST) ||
+          (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW) ||
+          (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHASLOW) ) //nuetral
+         hVeto->Fill(energy ) ;
+       
+       //fill number of primaries identified as ...
+       if(primaryCode >= 0) // Primary code defined
+         counter[recParticle->GetType()][primaryCode]++ ; 
+       
+      }
+      
+    } // no closest primary found
   }     
   
   
@@ -1074,24 +1094,18 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){
   hPrimary->Write(0,kOverwrite); 
   hAllRP->Write(0,kOverwrite);  
   hPhot->Write(0,kOverwrite);  
-  hPPSD->Write(0,kOverwrite);  
   hShape->Write(0,kOverwrite); 
   hVeto->Write(0,kOverwrite);  
   hPhotReg->Write(0,kOverwrite); 
   hPhotEM->Write(0,kOverwrite);   
-  hPhotPPSD->Write(0,kOverwrite); 
   hNReg ->Write(0,kOverwrite);  
   hNEM  ->Write(0,kOverwrite); 
-  hNPPSD->Write(0,kOverwrite);  
   hNBarReg ->Write(0,kOverwrite); 
   hNBarEM  ->Write(0,kOverwrite); 
-  hNBarPPSD->Write(0,kOverwrite); 
   hChargedReg ->Write(0,kOverwrite); 
   hChargedEM  ->Write(0,kOverwrite); 
-  hChargedPPSD->Write(0,kOverwrite); 
   hPbarReg ->Write(0,kOverwrite); 
   hPbarEM  ->Write(0,kOverwrite); 
-  hPbarPPSD->Write(0,kOverwrite); 
   
   cfile->Write(0,kOverwrite); 
   cfile->Close();