]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOSAnalyze.cxx
in dtor check the existence of fFile before to close it (the case occurs when merger...
[u/mrichter/AliRoot.git] / PHOS / AliPHOSAnalyze.cxx
index 5ee244b20084209a4abf4b619918bf228b4fd7e2..853ce1b29c0b0bf74f1111caa3f8142abcc056fe 100644 (file)
@@ -88,9 +88,7 @@
 #include "AliPHOSSDigitizer.h"
 #include "AliPHOSTrackSegment.h"
 #include "AliPHOSRecParticle.h"
-#include "AliPHOSIndexToObject.h"
 #include "AliPHOSCpvRecPoint.h"
-#include "AliPHOSPpsdRecPoint.h"
 #include "AliPHOSGetter.h"
 
 
@@ -129,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 ;
@@ -155,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 ;
@@ -169,161 +202,155 @@ 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 ;
-  const AliPHOSDigit * sdigit ;
-  
-  for(iSDigit = 0; iSDigit < gime->NSDigits(); iSDigit++)
+  AliPHOSDigit * sdigit ;
+  const TClonesArray * sdigits = gime->SDigits() ;
+  Int_t nsdig[5] = {0,0,0,0,0} ;
+  if(sdigits){
+    for(iSDigit = 0; iSDigit < sdigits->GetEntriesFast() ; iSDigit++)
       {
-       sdigit = gime->SDigit(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()->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 ;
-  const AliPHOSDigit * digit ;
-  for(iDigit = 0; iDigit < gime->NDigits(); iDigit++)
-    {
-      digit = gime->Digit(iDigit) ;
-      Int_t relid[4];
-      phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
-      Float_t x,z ;
-      phosgeom->RelPosInModule(relid,x,z) ;
-      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) ;
+  AliPHOSDigit * digit ;
+  const TClonesArray * digits = gime->Digits(); 
+  if(digits) {
+    for(iDigit = 0; iDigit < digits->GetEntriesFast(); 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()->Calibrate(digit->GetAmp()) ;
+       if(relid[0]==Nmod){
+         if(relid[1]==0)  //EMC
+           emcDigits->Fill(x,z,e) ;
+         if( relid[1]!=0 )
+           cpvDigits->Fill(x,z,e) ;
+       }
       }
-    }
-
-
+  }
+  
+  
   //Plot RecPoints
   Int_t irecp ;
   TVector3 pos ;
-  
-  for(irecp = 0; irecp < gime->NEmcRecPoints() ; irecp ++){
-    const AliPHOSEmcRecPoint * emc= gime->EmcRecPoint(irecp) ;
-    if(emc->GetPHOSMod()==Nmod){
-      emc->GetLocalPosition(pos) ;
-      emcOccupancy->Fill(pos.X(),pos.Z(),emc->GetEnergy());
-    }
-  }
-
-
-  for(irecp = 0; irecp < gime->NCpvRecPoints() ; irecp ++){
-    const AliPHOSRecPoint * cpv = gime->CpvRecPoint(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());
+  TObjArray * emcrp = gime->EmcRecPoints() ;
+  if(emcrp) {
+    for(irecp = 0; irecp < emcrp->GetEntriesFast() ; irecp ++){
+      AliPHOSEmcRecPoint * emc = (AliPHOSEmcRecPoint *) emcrp->At(irecp) ;
+      if(emc->GetPHOSMod()==Nmod){
+       emc->GetLocalPosition(pos) ;
+       emcRecPoints->Fill(pos.X(),pos.Z(),emc->GetEnergy());
       }
     }
-    else{
-      AliPHOSCpvRecPoint * cpv1 = (AliPHOSCpvRecPoint*) cpv ;
-      if(cpv1->GetPHOSMod()==Nmod){
-       cpv1->GetLocalPosition(pos) ;
-       ppsdUpCl->Fill(pos.X(),pos.Z(),cpv1->GetEnergy());
+  }
+  
+  TObjArray * cpvrp = gime->CpvRecPoints() ;
+  if(cpvrp) {
+    for(irecp = 0; irecp < cpvrp->GetEntriesFast() ; irecp ++){
+      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 ;
-  for(iRecParticle = 0; iRecParticle < gime->NRecParticles() ;iRecParticle++ )
-    {
-      recParticle =  gime->RecParticle(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 ;
-       
-
-       //extract list of primaries: it is stored at EMC RecPoints
-       Int_t emcIndex = gime->TrackSegment(recParticle->GetPHOSTSIndex())->GetEmcIndex() ;
-       Int_t numberofprimaries ;
-       Int_t * listofprimaries  = gime->EmcRecPoint(emcIndex)->GetPrimaries(numberofprimaries)  ;
-       Int_t index ;
-       const TParticle * primary ;
-       Double_t distance = minDistance ;
-       
-       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)
-           distance = TMath::Sqrt((recX-primX)*(recX-primX)+(recZ-primZ)*(recZ-primZ) ) ;
-         if(minDistance > distance)
-           {
-             minDistance = distance ;
-             closestPrimary = listofprimaries[index] ;
-           }
-       }
-       
-       if(closestPrimary >=0 ){
+  TClonesArray * rp = gime->RecParticles() ;
+  TClonesArray * ts = gime->TrackSegments() ;
+  if(rp && ts && emcrp) {
+    for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ; 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){
          
-         Int_t primaryType = gime->Primary(closestPrimary)->GetPdgCode() ;
+         Double_t minDistance = 5. ;
+         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  = ((AliPHOSRecPoint*) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries)  ;
+         Int_t index ;
+         const TParticle * primary ;
+         Double_t distance = minDistance ;
          
-         if(primaryType==22)
-           recPhot->Fill(recZ,recX,recParticle->Energy()) ;
-         else
-           if(primaryType==-2112)
-             recNbar->Fill(recZ,recX,recParticle->Energy()) ; 
+         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)
+             distance = TMath::Sqrt((recX-primX)*(recX-primX)+(recZ-primZ)*(recZ-primZ) ) ;
+           if(minDistance > distance)
+             {
+               minDistance = distance ;
+               closestPrimary = listofprimaries[index] ;
+             }
+         }
+         
+         if(closestPrimary >=0 ){
+           
+           Int_t primaryType = gime->Primary(closestPrimary)->GetPdgCode() ;
+           
+           if(primaryType==22)
+             recPhot->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") ;
   
 }
 //____________________________________________________________________________
@@ -414,27 +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 ;
-    for(iRecParticle = 0; iRecParticle < gime->NRecParticles()  ;iRecParticle++ )
+    TClonesArray * rp = gime->RecParticles() ;
+    if(!rp){
+      cout << "AliPHOSAnalyze::InvariantMass --> Can't find RecParticles " << endl ;
+      return ;
+    }
+
+    for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ )
       {
-       recParticle = gime->RecParticle(iRecParticle) ;
-       if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
-          (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM))
+       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 ;
@@ -464,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
          
@@ -533,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
@@ -541,8 +573,27 @@ void AliPHOSAnalyze::Ls(){
 
     const AliPHOSRecParticle * recParticle ;
     Int_t iRecParticle ;
-    for(iRecParticle = 0; iRecParticle < gime->NRecParticles() ;iRecParticle++ ){
-      recParticle = gime->RecParticle(iRecParticle) ;
+    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 ;
@@ -553,10 +604,10 @@ void AliPHOSAnalyze::Ls(){
       Int_t closestPrimary = -1 ;
       
       //extract list of primaries: it is stored at EMC RecPoints
-      Int_t emcIndex = gime->TrackSegment(recParticle->GetPHOSTSIndex())->GetEmcIndex() ;
+      Int_t emcIndex = ((AliPHOSTrackSegment*) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
       Int_t numberofprimaries ;
-      Int_t * listofprimaries  = gime->EmcRecPoint(emcIndex)->GetPrimaries(numberofprimaries)  ;
-
+      Int_t * listofprimaries  = ((AliPHOSEmcRecPoint*) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries)  ;
+      
       Int_t index ;
       const TParticle * primary ;
       Double_t distance = minDistance ;
@@ -586,12 +637,12 @@ void AliPHOSAnalyze::Ls(){
        const TParticle * primary = gime->Primary(closestPrimary) ;
        if(primary->GetPdgCode() == 22){
          hAllEnergy->Fill(primary->Energy(), recParticle->Energy()) ;
-         if(recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA){
+         if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){
            hPhotEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; 
            hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; 
          }
          else
-           if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)
+           if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)
              hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; 
        }
       }
@@ -652,17 +703,35 @@ 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 ;
-    for(iRecParticle = 0; iRecParticle < gime->NRecParticles() ;iRecParticle++ ){
-      recParticle = gime->RecParticle(iRecParticle) ;
+    for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ ){
+      recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
       
       //find the closest primary
       Int_t moduleNumberRec ;
@@ -673,9 +742,9 @@ void AliPHOSAnalyze::PositionResolution(const char * branchTitle)
       Int_t closestPrimary = -1 ;
       
       //extract list of primaries: it is stored at EMC RecPoints
-      Int_t emcIndex = gime->TrackSegment(recParticle->GetPHOSTSIndex())->GetEmcIndex() ;
+      Int_t emcIndex = ((AliPHOSTrackSegment*) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
       Int_t numberofprimaries ;
-      Int_t * listofprimaries  = gime->EmcRecPoint(emcIndex)->GetPrimaries(numberofprimaries)  ;
+      Int_t * listofprimaries  = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries)  ;
 
       Int_t index ;
       const TParticle * primary ;
@@ -709,12 +778,12 @@ void AliPHOSAnalyze::PositionResolution(const char * branchTitle)
          hAllPosition->Fill(primary->Energy(), minDistance) ;
          hAllPositionX->Fill(primary->Energy(), dX) ;
          hAllPositionZ->Fill(primary->Energy(), dZ) ;
-         if(recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA){
+         if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){
            hPhotPosition->Fill(primary->Energy(), minDistance ) ; 
            hEMPosition->Fill(primary->Energy(), minDistance ) ; 
          }
          else
-           if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)
+           if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)
              hEMPosition->Fill(primary->Energy(), minDistance ) ; 
        }
       }
@@ -744,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
 
@@ -752,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
@@ -788,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.);
@@ -806,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");
@@ -817,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");
@@ -828,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");
@@ -839,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)
@@ -850,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
@@ -870,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 ;
@@ -892,11 +956,12 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){
       }
       
     }
+    
     //========== Now scan over RecParticles            
     const AliPHOSRecParticle * recParticle ;
     Int_t iRecParticle ;
-    for(iRecParticle = 0; iRecParticle < gime->NRecParticles(); iRecParticle++ ){
-      recParticle = gime->RecParticle(iRecParticle) ;
+    for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ ){
+      recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
       //fill histo spectrum of all RecParticles
       hAllRP->Fill(CorrectedEnergy(recParticle->Energy())) ;
       
@@ -909,9 +974,9 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){
       Int_t closestPrimary = -1 ;
       
       //extract list of primaries: it is stored at EMC RecPoints
-      Int_t emcIndex = gime->TrackSegment(recParticle->GetPHOSTSIndex())->GetEmcIndex() ;
+      Int_t emcIndex = ((AliPHOSTrackSegment *) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
       Int_t numberofprimaries ;
-      Int_t * listofprimaries  = gime->EmcRecPoint(emcIndex)->GetPrimaries(numberofprimaries)  ;
+      Int_t * listofprimaries  = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries)  ;
       Int_t index ;
       const TParticle * primary ;
       Double_t distance = minDistance ;
@@ -960,7 +1025,7 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){
 
        //==========Now look at the type of RecParticle
        Float_t energy = CorrectedEnergy(recParticle->Energy()) ;
-       if(recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA){
+       if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){
          hPhot->Fill(energy ) ;        
          switch(primaryCode){
          case 0:
@@ -982,33 +1047,10 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){
            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((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
-          (recParticle->GetType() == AliPHOSFastRecParticle::kELECTRON)||
-          (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)||
-          (recParticle->GetType() == AliPHOSFastRecParticle::kABSURDEM) ){ //with EM shower
+       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:
@@ -1031,10 +1073,10 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){
          }
        }
        
-       if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
-          (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHA) ||
-          (recParticle->GetType() == AliPHOSFastRecParticle::kGAMMAHA) ||
-          (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM) ) //nuetral
+       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 ...
@@ -1052,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();