]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOSAnalyze.cxx
Dimitri just makes it work
[u/mrichter/AliRoot.git] / PHOS / AliPHOSAnalyze.cxx
index 6a4f97ad92a24895bb5373eb6451d7c6f098333f..b8c9a84a07ce9358985f659d8afa857798ed4635 100644 (file)
@@ -87,6 +87,11 @@ AliPHOSAnalyze::AliPHOSAnalyze(Text_t * name)
       fEvt = -999 ; 
 
   }
+  fClu = 0 ;
+  fPID = 0 ;
+  fTrs = 0 ;
+  fRec = 0 ;
+  ResetHistograms() ;
 }
 
 //____________________________________________________________________________
@@ -111,73 +116,110 @@ AliPHOSAnalyze::~AliPHOSAnalyze()
 
   if (fRootFile->IsOpen() ) 
     fRootFile->Close() ; 
-  delete fRootFile ; 
-  fRootFile = 0 ; 
+  if(fRootFile)
+    delete fRootFile ;  
 
-  delete fPHOS ; 
-  fPHOS = 0 ; 
+  if(fPHOS)
+    delete fPHOS ; 
 
-  delete fClu ; 
-  fClu = 0 ;
+  if(fClu)
+    delete fClu ; 
 
-  delete fPID ; 
-  fPID = 0 ;
+  if(fPID)
+    delete fPID ; 
 
-  delete fRec ; 
-  fRec = 0 ;
+  if(fRec)
+    delete fRec ; 
 
-  delete fTrs ; 
-  fTrs = 0 ;
+  if(fTrs)
+    delete fTrs ; 
 
 }
 
 //____________________________________________________________________________
-void AliPHOSAnalyze::AnalyzeOneEvent(Int_t evt)
-{
-  // analyze one single event with id=evt
-
-  TStopwatch ts ; 
-  ts.Start() ; 
-  Bool_t ok = Init(evt) ; 
+void AliPHOSAnalyze::ActivePPSD(Int_t Nevents=1){
   
-  if ( ok ) {
-    //=========== Get the number of entries in the Digits array
-    
-    Int_t nId = fPHOS->Digits()->GetEntries();    
-    printf("AnalyzeOneEvent > Number of entries in the Digit array is %d \n",nId);
-    
-    //=========== Do the reconstruction
-    
-    cout << "AnalyzeOneEvent > Found  " << nId << "  digits in PHOS"   << endl ;  
-    
-    fPHOS->Reconstruction(fRec);  
-
-    // =========== End of reconstruction
-
-    // Deleting fClu, fTrs, fPID et fRec
-    fClu->Delete();
-    fTrs->Delete();
-    fPID->Delete();
-    fRec->Delete();
-
-    // =========== Write the root file
+  fhEnergyCorrelations  = new TH2F("hEnergyCorrelations","hEnergyCorrelations",40,  0., 0.15, 30, 0., 3.e-5);
+  //========== 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;
   
-    // =========== Finish
+  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);
 
-    cout << "AnalyzeOneEvent > event # " << fEvt << " processed" << endl ;   
-  } // ok
-  else
-    cout << "AnalyzeOneEvent > filed to process event # " << evt << endl ;   
-  
-  ts.Stop() ;   cout << "CPU  time = " << ts.CpuTime()   << endl ; 
-  cout << "Real time = " << ts.RealTime()  << endl ; 
+      //=========== 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") ;
 }
 
+
 //____________________________________________________________________________
- void AliPHOSAnalyze::AnalyzeManyEvents(Int_t Nevents, Int_t module)    
+void AliPHOSAnalyze::AnalyzeManyEvents(Int_t Nevents, Int_t module)    
 {
   // analyzes Nevents events in a single PHOS module  
+  // Events should be reconstructed by Reconstruct()
 
   if ( fRootFile == 0 ) 
     cout << "AnalyzeManyEvents > " << "Root File not openned" << endl ;  
@@ -186,84 +228,70 @@ void AliPHOSAnalyze::AnalyzeOneEvent(Int_t evt)
       //========== Booking Histograms
       cout << "AnalyzeManyEvents > " << "Booking Histograms" << endl ; 
       BookingHistograms();
+
       Int_t ievent;
       Int_t relid[4] ; 
       AliPHOSDigit * digit ;
-//       AliPHOSEmcRecPoint * emc ;
-//       AliPHOSPpsdRecPoint * ppsd ;
+      AliPHOSEmcRecPoint * emc ;
+      AliPHOSPpsdRecPoint * ppsd ;
       //      AliPHOSTrackSegment * tracksegment ;
       AliPHOSRecParticle * recparticle;
+
       for ( ievent=0; ievent<Nevents; ievent++)
        {  
-          if (ievent==0) {
-           cout << "AnalyzeManyEvents > " << "Starting Analyzing " << endl ; 
-           //========== Create the Clusterizer
-           fClu = new AliPHOSClusterizerv1() ; 
-           fClu->SetEmcEnergyThreshold(0.05) ; 
-           fClu->SetEmcClusteringThreshold(0.50) ; 
-           fClu->SetPpsdEnergyThreshold    (0.0000002) ; 
-           fClu->SetPpsdClusteringThreshold(0.0000001) ; 
-           fClu->SetLocalMaxCut(0.03) ;
-           fClu->SetCalibrationParameters(0., 0.00000001) ;  
-           //========== Creates the track segment maker
-           fTrs = new AliPHOSTrackSegmentMakerv1()  ;
-           //    fTrs->UnsetUnfoldFlag() ; 
-           //========== Creates the particle identifier
-           fPID = new AliPHOSPIDv1() ;
-           fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ; 
-           fPID->Print() ;         
-           //========== Creates the Reconstructioner  
-           fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
-         }
          //========== Event Number>         
-         //  if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. ) 
-         cout <<  "AnalyzeManyEvents > " << "Event is " << ievent << endl ;  
+         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()));
-             }
+             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()));
+               }
            }
-         //=========== Do the reconstruction
-         fPHOS->Reconstruction(fRec);
          
-         //      //=========== Cluster in module
-         //      TIter nextEmc(fPHOS->EmcRecPoints()  ) ;
-         //      while((emc = (AliPHOSEmcRecPoint *)nextEmc())) 
-         //        {
-         //          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()) ;
-//             }
-//         }
+
+         //=========== Cluster in module
+         TIter nextEmc(*fPHOS->EmcRecPoints()  ) ;
+         while((emc = (AliPHOSEmcRecPoint *)nextEmc())) 
+           {
+             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()) ;
+               }
+           }
+
          //========== TRackSegments in the event
          TIter nextRecParticle(*fPHOS->RecParticles() ) ; 
          while((recparticle = (AliPHOSRecParticle *)nextRecParticle())) 
@@ -277,69 +305,29 @@ void AliPHOSAnalyze::AnalyzeOneEvent(Int_t evt)
                  Int_t index ;
                  for ( index = 0 ; index < numberofprimaries ; index++)
                    cout << "    primary # " << index << " =  " << listofprimaries[index] << endl ;  
-//               switch(recparticle->GetType())
-//                 {
-//                 case AliPHOSFastRecParticle::kGAMMA:
-//                   fhPhotonEnergy->Fill(recparticle->Energy() ) ; 
-//                   //fhPhotonPositionX->Fill(recpart. ) ;
-//                   //fhPhotonPositionY->Fill(recpart. ) ;                 
-//                   cout << "PHOTON" << endl;
-//                   break;
-//                 case  AliPHOSFastRecParticle::kELECTRON:
-//                   fhElectronEnergy->Fill(recparticle->Energy() ) ; 
-//                   //fhElectronPositionX->Fill(recpart. ) ;
-//                   //fhElectronPositionY->Fill(recpart. ) ; 
-//                   cout << "ELECTRON" << endl;
-//                   break;
-//                 case  AliPHOSFastRecParticle::kNEUTRALHA:
-//                   fhNeutralHadronEnergy->Fill(recparticle->Energy() ) ; 
-//                   //fhNeutralHadronPositionX->Fill(recpart. ) ;
-//                   //fhNeutralHadronPositionY->Fill(recpart. ) ; 
-//                   cout << "NEUTRAl HADRON" << endl;
-//                   break ;
-//                 case  AliPHOSFastRecParticle::kNEUTRALEM:
-//                   fhNeutralEMEnergy->Fill(recparticle->Energy() ) ; 
-//                   //fhNeutralEMPositionX->Fill(recpart. ) ;
-//                   //fhNeutralEMPositionY->Fill(recpart. ) ; 
-//                   //cout << "NEUTRAL EM" << endl;
-//                   break ;
-//                 case  AliPHOSFastRecParticle::kCHARGEDHA:
-//                   fhChargedHadronEnergy->Fill(recparticle->Energy() ) ; 
-//                   //fhChargedHadronPositionX->Fill(recpart. ) ;
-//                   //fhChargedHadronPositionY->Fill(recpart. ) ; 
-//                   cout << "CHARGED HADRON" << endl;
-//                   break ;
-//                 case  AliPHOSFastRecParticle::kGAMMAHA:
-//                   fhPhotonHadronEnergy->Fill(recparticle->Energy() ) ; 
-//                   //fhPhotonHadronPositionX->Fill(recpart. ) ;
-//                   //fhPhotonHadronPositionY->Fill(recpart. ) ; 
-//                   cout << "PHOTON HADRON" << endl;
-//                   break ;                 
-//                 }
                }
            }
        }   // endfor
-      SavingHistograms();
+      SaveHistograms();
     }       // endif
 }           // endfunction
 
 //____________________________________________________________________________
- void AliPHOSAnalyze::Reconstruct(Int_t Nevents )    
-{  
-   
+ void AliPHOSAnalyze::Reconstruct(Int_t Nevents,Int_t FirstEvent )    
+{     
   Int_t ievent ;   
-  for ( ievent=0; ievent<Nevents; ievent++)
+  for ( ievent=FirstEvent; ievent<Nevents; ievent++)
     {  
-      if (ievent==0
+      if (ievent==FirstEvent
        {
-         cout << "Analyze > " << "Starting Analyzing " << endl ; 
+         cout << "Analyze > Starting Reconstructing " << endl ; 
          //========== Create the Clusterizer
          fClu = new AliPHOSClusterizerv1() ; 
-         fClu->SetEmcEnergyThreshold(0.05) ; 
+         fClu->SetEmcEnergyThreshold(0.03) ; 
          fClu->SetEmcClusteringThreshold(0.20) ; 
-         fClu->SetPpsdEnergyThreshold    (0.0000002) ; 
-         fClu->SetPpsdClusteringThreshold(0.0000001) ; 
-         fClu->SetLocalMaxCut(0.03) ;
+         fClu->SetPpsdEnergyThreshold    (0.0000001) ; 
+         fClu->SetPpsdClusteringThreshold(0.0000001) ;
+         fClu->SetLocalMaxCut(0.02) ;
          fClu->SetCalibrationParameters(0., 0.00000001) ; 
          
          //========== Creates the track segment maker
@@ -349,16 +337,18 @@ void AliPHOSAnalyze::AnalyzeOneEvent(Int_t evt)
          //========== Creates the particle identifier
          fPID = new AliPHOSPIDv1() ;
          fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ; 
+         fPID->SetDispersionCutOff(2.0) ;
+         fPID->SetRelativeDistanceCut(3.) ;
          
          //========== Creates the Reconstructioner  
          fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ; 
-         //  fRec -> SetDebugReconstruction(kTRUE);     
+         //              fRec -> SetDebugReconstruction(kTRUE);     
 
        }
       
       //========== Event Number>         
-      if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. ) 
-       cout <<  "Analyze > " << "Event is " << ievent << endl ;  
+      //      if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. ) 
+       cout <<  "Reconstruct > Event is " << ievent << endl ;  
       
       //=========== Connects the various Tree's for evt
       gAlice->GetEvent(ievent);
@@ -371,91 +361,132 @@ void AliPHOSAnalyze::AnalyzeOneEvent(Int_t evt)
     }
 
     fClu->Delete();
+    fClu=0 ;
     fTrs->Delete();
+    fTrs = 0 ;
     fPID->Delete();
+    fPID = 0 ;
     fRec->Delete();
+    fRec = 0 ;
 
 }
-//-------------------------------------------------------------------------------------
-
-//   TClonesArray AllDigitArray = TClonesArray("AliPHOSDigit",1000) ;
-//   TClonesArray * PhotonsList ;
-//   TClonesArray * FalsDigitsList ;
-//   TClonesArray AllPrimary = TClonesArray("TParticle",5000) ;
-//   TFile * file2 = new TFile("ph100.root") ; // file with added photons
-//   gAlice = (AliRun*) file2->Get("gAlice") ;
-//   Int_t ievent;
-//   Int_t NDigits[Nevents+1] ;
-//   NDigits[0]=0 ;
-//   Int_t NAllDigits = 0;
-//   Int_t NprimPerEvent = 20 ;
-//   for (ievent=0; ievent <Nevents; ievent++)
-//     {
-//       PhotonsList  = gAlice->Particles();  //Primary
-//       FalsDigitsList  = ((AliPHOSv1 *)gAlice->GetDetector("PHOS"))->Digits();  //Digits
-//       gAlice->GetEvent(ievent) ;
-//       gAlice->TreeD()->GetEvent(0) ;
-//       gAlice->TreeK()->GetEvent(0) ;
-//       //Copy Primary
-//       Int_t Nprim ;
-//       for(Nprim = 0 ;Nprim < NprimPerEvent ; Nprim++)
-//     new (AllPrimary[Nprim+ievent*NprimPerEvent])  TParticle(*((TParticle *) PhotonsList->At(Nprim))) ;
-
-//       //Copy Digits
-//       TIter nextDigit(FalsDigitsList) ;
-//       AliPHOSDigit * FalseDigit ;
-//       NDigits[ievent+1] = NDigits[ievent]+ FalsDigitsList->GetEntriesFast() ; 
-//       while( (FalseDigit = (AliPHOSDigit *) nextDigit()))
-//     {        
-//       new (AllDigitArray[NAllDigits])  AliPHOSDigit(FalseDigit->GetPrimary(1),FalseDigit->GetId(),FalseDigit->GetAmp()) ;
-//       NAllDigits++ ;
-//     }         
-//     }
-//   file2->Close() ;
-
-
-
-//       //Add primary particles
-//       cout << "# of Primaries before add " << PrimaryList->GetEntriesFast() << endl;
-//      Int_t NTruePrimary = 0 ;  //PrimaryList->GetEntriesFast() ;
-//       Int_t Nprim ;
-//       for(Nprim = 0; Nprim < NprimPerEvent; Nprim++)
-//     new ((*PrimaryList)[NTruePrimary+Nprim])  TParticle(*((TParticle *) AllPrimary.At(Nprim+ievent*NprimPerEvent))) ;
+//____________________________________________________________________________
+ void AliPHOSAnalyze::InvariantMass(Int_t Nevents )    
+{
+  // Calculates Real and Mixed invariant mass distributions
+  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.) ;
+  TH2D * hRealPhot = new TH2D("hRealPhot", "Real for kPhoton particles", 250,0.,1.,40,0.,4.) ;
+  TH2D * hMixedEM  = new TH2D("hMixedEM",  "Mixed for EM particles",     250,0.,1.,40,0.,4.) ;
+  TH2D * hMixedPhot= new TH2D("hMixedPhot","Mixed for kPhoton particles",250,0.,1.,40,0.,4.) ;
+  
+  Int_t ievent;
+  Int_t EventInMixedLoop ;
+  
+  Int_t NRecParticles[NMixedEvents] ;
+  
+  AliPHOSRecParticle::RecParticlesList * AllRecParticleList  = new TClonesArray("AliPHOSRecParticle", NMixedEvents*1000) ;
+  
+  for(EventInMixedLoop = 0; EventInMixedLoop < MixedLoops; EventInMixedLoop++  ){
+    Int_t iRecPhot = 0 ;
+    
+    for ( ievent=0; ievent < NMixedEvents; ievent++){        
       
-//       cout << "# of Primaries after add " << PrimaryList->GetEntriesFast() <<endl;
-
-//       cout << "Digits before add " << DigitsList->GetEntries() << endl ;
-//       cout << "Digits to add " <<  NDigits[ievent+1]-  NDigits[ievent]<< endl ;
+      Int_t AbsEventNumber = EventInMixedLoop*NMixedEvents + ievent ;
       
-      //=========== Add fals digits ==============================
-//       TIter nextDigit(DigitsList) ;
-//       AliPHOSDigit * FalseDigit ;
-//       AliPHOSDigit * RealDigit ;
-//       Int_t NTrueDigits = DigitsList->GetEntriesFast() ; 
-//       Int_t Ndigit ;
-//       for(Ndigit=NDigits[ievent];Ndigit<NDigits[ievent+1];Ndigit++)
-//     {        
-//       FalseDigit = (AliPHOSDigit*) AllDigitArray.At(Ndigit) ;
-//       Bool_t Add = kTRUE ; 
-//       AliPHOSDigit tmpDigit=AliPHOSDigit(FalseDigit->GetPrimary(1)+NTruePrimary,FalseDigit->GetId(),FalseDigit->GetAmp()) ;
-         
-//       while( (RealDigit = (AliPHOSDigit *) nextDigit()) && Add)
-//         {
-//           if((*RealDigit) == (tmpDigit)) 
-//             {
-//               *RealDigit=*RealDigit+tmpDigit ;
-//               Add = kFALSE ;
-//             }
-//         }
-//       if(Add)
-//         {
-//           new ((*DigitsList)[NTrueDigits])  AliPHOSDigit(FalseDigit->GetPrimary(1)+NTruePrimary,FalseDigit->GetId(),FalseDigit->GetAmp()) ;
-//           ((AliPHOSDigit *)DigitsList->At(NTrueDigits))->SetIndexInList(NTrueDigits) ;
-//           NTrueDigits++ ;
-//         }     
-//     }
-//       cout << "Digits after add " << DigitsList->GetEntries() << endl ;
+      //=========== Connects the various Tree's for evt
+      gAlice->GetEvent(AbsEventNumber);
+      
+      //=========== Get the Digit Tree
+      gAlice->TreeD()->GetEvent(0) ;
+      
+      //========== Creating branches ===================================       
+      
+      AliPHOSRecParticle::RecParticlesList ** RecParticleList  = fPHOS->RecParticles() ;
+      if( (*RecParticleList) )
+       (*RecParticleList)->Clear() ;
+      gAlice->TreeR()->SetBranchAddress( "PHOSRP", RecParticleList ) ;
+      
+      //=========== Gets the Reconstraction TTree
+      gAlice->TreeR()->GetEvent(0) ;
+      
+      AliPHOSRecParticle * RecParticle ;
+      Int_t 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) ;
+           iRecPhot++;
+         }
+       }
+      
+       NRecParticles[ievent] = iRecPhot-1 ;  
+    }
+    
+
+    //Now calculate invariant mass:
+    Int_t irp1,irp2 ;
+    Int_t NCurEvent = 0 ;
+
+    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) ;
+           
+       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);
+       
+       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(irp2 <= NRecParticles[NCurEvent]){ //'Real' event
+         hRealEM->Fill(InvMass,Pt);
+         if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
+           hRealPhot->Fill(InvMass,Pt);
+       }
+       else{
+         hMixedEM->Fill(InvMass,Pt);
+         if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
+           hMixedPhot->Fill(InvMass,Pt);
+       } //real-mixed
+           
+      } //loop over second rp
+    }//loop over first rp
+
+
+    AllRecParticleList->Delete() ;
+  } //Loop over events
+  
+  delete AllRecParticleList ;
+  
+  //writing output
+  TFile output("invmass.root","RECREATE");
+  output.cd();
+  
+  hRealEM->Write() ;
+  hRealPhot->Write() ;
+  hMixedEM->Write() ;
+  hMixedPhot->Write() ;
+  
+  output.Write();
+  output.Close();
 
+}
 
 //____________________________________________________________________________
  void AliPHOSAnalyze::AnalyzeResolutions(Int_t Nevents )    
@@ -483,15 +514,13 @@ void AliPHOSAnalyze::AnalyzeOneEvent(Int_t evt)
     {  
       
       //========== Event Number>         
-      if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. ) 
+      //      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
+      //=========== Gets the Kine TTree
       gAlice->TreeK()->GetEvent(0) ;
       
       //=========== Gets the list of Primari Particles
@@ -500,23 +529,31 @@ void AliPHOSAnalyze::AnalyzeOneEvent(Int_t evt)
       TParticle * Primary ;
       Int_t iPrimary ;
       for ( iPrimary = 0 ; iPrimary < PrimaryList->GetEntries() ; iPrimary++)
-           {
-             Primary = (TParticle*)PrimaryList->UncheckedAt(iPrimary) ;
-             Int_t PrimaryType = Primary->GetPdgCode() ;
-             if( PrimaryType == 22 ) 
-               fhPrimary->Fill(Primary->Energy()) ;
-           } 
+       {
+         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() ;
@@ -526,22 +563,22 @@ void AliPHOSAnalyze::AnalyzeOneEvent(Int_t evt)
       if( (*RecParticleList) )
        (*RecParticleList)->Clear() ;
       gAlice->TreeR()->SetBranchAddress( "PHOSRP", RecParticleList ) ;
-           
-
+      
       //=========== Gets the Reconstraction TTree
       gAlice->TreeR()->GetEvent(0) ;
+      
       AliPHOSRecParticle * RecParticle ;
       Int_t iRecParticle ;
       for(iRecParticle = 0; iRecParticle < (*RecParticleList)->GetEntries() ;iRecParticle++ )
        {
          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) ;
          
-         Double_t MinDistance = 10000 ;
+         Double_t MinDistance = 5. ;
          Int_t ClosestPrimary = -1 ;
          
          Int_t numberofprimaries ;
@@ -549,160 +586,209 @@ void AliPHOSAnalyze::AnalyzeOneEvent(Int_t evt)
          Int_t index ;
          TParticle * Primary ;
          Double_t Distance = MinDistance ;
-         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] ;
-               }
-           }
+         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] ;
+             }
+         }
          TotalRecPart++ ;
-         if(ClosestPrimary >=0 )
-           {
-             fhPhotonAllEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ; 
-             fhPhotonAllPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),Distance) ;
-             TotalRPwithPrim++;
-             Int_t PrimaryType = ((TParticle *)PrimaryList->At(ClosestPrimary))->GetPdgCode() ;
-             TParticlePDG* PDGparticle = ((TParticle *)PrimaryList->At(ClosestPrimary))->GetPDG();
-             Double_t charge =  PDGparticle->Charge() ;
-             Int_t PrimaryCode ;
-             switch(PrimaryType)
-               {
-               case 22:
-                 PrimaryCode = 0;  //Photon
-                 break;
-               case 11 :
-                 PrimaryCode = 1;  //Electron
-                 break;
-               case -11 :
-                 PrimaryCode = 1;  //positron
-                 break;
-               case 321 :
-                 PrimaryCode = 4;  //K+
-                 break;
-               case -321 :
-                 PrimaryCode = 4;  //K-
-                 break;
-               case 310 :
-                 PrimaryCode = 4;  //K0s
-                 break;
-               case 130 :
-                 PrimaryCode = 4;  //K0l
-                 break;
-               default:
-                 if(charge)
-                   PrimaryCode = 2; //Charged hadron
-                 else
-                   PrimaryCode = 3; //Neutral hadron
-                 break;
-               }
 
-             switch(RecParticle->GetType())
-               {
-               case AliPHOSFastRecParticle::kGAMMA:
-                 if(PrimaryType == 22){
-                   fhPhotonEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ; 
-                   fhPhotonPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),Distance) ;
-                   fhPhotonReg->Fill(RecParticle->Energy() ) ;
-                   fhPhotonEM->Fill(RecParticle->Energy() ) ;
-                   fhPhotPhot->Fill(RecParticle->Energy() ) ;
+         if(ClosestPrimary >=0 ){
+           TotalRPwithPrim++;
+           
+           Int_t PrimaryType = ((TParticle *)PrimaryList->At(ClosestPrimary))->GetPdgCode() ;
+           TParticlePDG* PDGparticle = ((TParticle *)PrimaryList->At(ClosestPrimary))->GetPDG();
+           Double_t charge =  PDGparticle->Charge() ;
+           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) ;
+               break;
+             case 11 :
+               PrimaryCode = 1;  //Electron
+               break;
+             case -11 :
+               PrimaryCode = 1;  //positron
+               break;
+             case 321 :
+               PrimaryCode = 4;  //K+
+               break;
+             case -321 :
+               PrimaryCode = 4;  //K-
+               break;
+             case 310 :
+               PrimaryCode = 4;  //K0s
+               break;
+             case 130 :
+               PrimaryCode = 4;  //K0l
+               break;
+             default:
+               if(charge)
+                 PrimaryCode = 2; //Charged hadron
+               else
+                 PrimaryCode = 3; //Neutral hadron
+               break;
+             }
+           
+           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() ) ; 
+
+                 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()) ) ;
+
+                 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 ~
+                 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()) ) ;
+               }
+               
+               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 == 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(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]++;
+               break;
+             case  AliPHOSFastRecParticle::kNEUTRALHA:
+               if(PrimaryType == 22) 
+                 fhPhotNeuH->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+
+               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 ) ;
+               
+                 fhPhotNuEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+                 fhPhotEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+               }
+               if(PrimaryType == 2112) //neutron
+                 fhNEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+               
+               if(PrimaryType == -2112) //neutron ~
+                 fhNBarEM->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()) ) ;
+
+               Counter[3][PrimaryCode]++;
+               break ;
+             case  AliPHOSFastRecParticle::kCHARGEDHA:
+               if(PrimaryType == 22) //photon
+                 fhPhotChHa->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+               
+               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 == 2112){ //neutron
-                   fhNReg->Fill(RecParticle->Energy() ) ;
-                   fhNEM->Fill(RecParticle->Energy() ) ;
+                   fhNPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
                  }
-                 
+               
                  if(PrimaryType == -2112){ //neutron ~
-                   fhNBarReg->Fill(RecParticle->Energy() ) ;
-                   fhNBarEM->Fill(RecParticle->Energy() ) ;
-
+                   fhNBarPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ; 
                  }
                  if(PrimaryCode == 2){
-                   fhChargedReg->Fill(RecParticle->Energy() ) ;
-                   fhChargedEM->Fill(RecParticle->Energy() ) ;
+                   fhChargedPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
                  }
+               
+                 fhAllPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+                 fhVeto->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+                 fhPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 
-                 fhAllReg->Fill(RecParticle->Energy() ) ;
-                 fhAllEM->Fill(RecParticle->Energy() ) ;
-                 Counter[0][PrimaryCode]++;
-                 break;
-               case  AliPHOSFastRecParticle::kELECTRON:
-                 if(PrimaryType == 11 || PrimaryType == -11){
-                   fhElectronEnergy->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ; 
-                   fhElectronPosition->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),Distance ) ;
-                 }
-                 if(PrimaryType == 22) 
-                   fhPhotElec->Fill(RecParticle->Energy() ) ;
-
-                 Counter[1][PrimaryCode]++;
-                 break;
-               case  AliPHOSFastRecParticle::kNEUTRALHA:
-                 if(PrimaryType == 22) 
-                   fhPhotNeuH->Fill(RecParticle->Energy() ) ;
-
-                 fhNeutralHadronEnergy->Fill( ((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy()) ; 
-                 fhNeutralHadronPosition->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy() ,Distance  ) ;
-                 Counter[2][PrimaryCode]++;
-                 break ;
-               case  AliPHOSFastRecParticle::kNEUTRALEM:
-                 if(PrimaryType == 22 || PrimaryType == 11 || PrimaryType == -11){
-                   fhNeutralEMEnergy->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),RecParticle->Energy() ) ; 
-                   fhNeutralEMPosition->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),Distance ) ;
-                 }
-
-                 if(PrimaryType == 22){ //photon
-                   fhPhotNuEM->Fill(RecParticle->Energy() ) ;
-                   fhPhotonEM->Fill(RecParticle->Energy() ) ;
-                 }
-                 if(PrimaryType == 2112) //neutron
-                   fhNEM->Fill(RecParticle->Energy() ) ;
-                 
-                 if(PrimaryType == -2112) //neutron ~
-                   fhNBarEM->Fill(RecParticle->Energy() ) ;
-
-                 if(PrimaryCode == 2)
-                   fhChargedEM->Fill(RecParticle->Energy() ) ;
-
-                 fhAllEM->Fill(RecParticle->Energy() ) ;
-
-                 Counter[3][PrimaryCode]++;
-                 break ;
-               case  AliPHOSFastRecParticle::kCHARGEDHA:
-                 if(PrimaryType == 22) //photon
-                   fhPhotChHa->Fill(RecParticle->Energy() ) ;
-                 
-                 fhChargedHadronEnergy->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),RecParticle->Energy() ) ; 
-                 fhChargedHadronPosition->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),Distance ) ;
-                 Counter[4][PrimaryCode]++ ;
-                 break ;
-               case  AliPHOSFastRecParticle::kGAMMAHA:
-                 if(PrimaryType == 22) //photon
-                   fhPhotGaHa->Fill(RecParticle->Energy() ) ;
-                 fhPhotonHadronEnergy->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy()) ; 
-                 fhPhotonHadronPosition->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),Distance ) ;
                  Counter[5][PrimaryCode]++ ;
                  break ;       
-               case  AliPHOSFastRecParticle::kABSURDEM:              
-                 Counter[6][PrimaryCode]++ ;
-                 break;
-               case  AliPHOSFastRecParticle::kABSURDHA:
-                 Counter[7][PrimaryCode]++ ;
-                 break;
-               default:
-                 Counter[8][PrimaryCode]++ ;
-                 break;
-               }
-           }
+             case  AliPHOSFastRecParticle::kABSURDEM:        
+               Counter[6][PrimaryCode]++ ;
+               fhShape->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+               break;
+             case  AliPHOSFastRecParticle::kABSURDHA:
+               Counter[7][PrimaryCode]++ ;
+               break;
+             default:
+               Counter[8][PrimaryCode]++ ;
+               break;
+             }
+         }
        }  
     }   // endfor
-  SaveResolutionHistograms();
+  SaveHistograms();
   cout << "Resolutions: Analyzed " << Nevents << " event(s)" << endl ;
   cout << "Resolutions: Total primary       " << TotalPrimary << endl ;
   cout << "Resoluitons: Total reconstracted " << TotalRecPart << endl ;
@@ -739,32 +825,13 @@ void  AliPHOSAnalyze::BookingHistograms()
   delete fhConvertorCluster  ;
   delete fhConvertorEmc  ;
   
-//   fhEmcDigit                = new TH1F("hEmcDigit",      "hEmcDigit",         1000,  0. ,  25.);
-//   fhVetoDigit               = new TH1F("hVetoDigit",     "hVetoDigit",         500,  0. ,  3.e-5);
-//   fhConvertorDigit          = new TH1F("hConvertorDigit","hConvertorDigit",    500,  0. ,  3.e-5);
-//   fhEmcCluster              = new TH1F("hEmcCluster",    "hEmcCluster",       1000,  0. ,  30.);
-//   fhVetoCluster             = new TH1F("hVetoCluster",   "hVetoCluster",       500,  0. ,  3.e-5);
-//   fhConvertorCluster        = new TH1F("hConvertorCluster","hConvertorCluster",500,  0. ,  3.e-5);
-//   fhConvertorEmc            = new TH2F("hConvertorEmc",  "hConvertorEmc",      200,  1. ,  3., 200, 0., 3.e-5);
-//   fhPhotonEnergy            = new TH1F("hPhotonEnergy",  "hPhotonEnergy",     1000,  0. ,  30.);
-//   fhElectronEnergy          = new TH1F("hElectronEnergy","hElectronEnergy",   1000,  0. ,  30.);
-//   fhNeutralHadronEnergy     = new TH1F("hNeutralHadronEnergy", "hNeutralHadronEnergy",    1000,  0. ,  30.);
-//   fhNeutralEMEnergy         = new TH1F("hNeutralEMEnergy", "hNeutralEMEnergy",    1000,  0. ,  30.);
-//   fhChargedHadronEnergy     = new TH1F("hChargedHadronEnergy", "hChargedHadronEnergy",    1000,  0. ,  30.);
-//   fhPhotonHadronEnergy      = new TH1F("hPhotonHadronEnergy","hPhotonHadronEnergy",500,-80. , 80.);
-//   fhPhotonPositionX         = new TH1F("hPhotonPositionX","hPhotonPositionX",   500,-80. , 80.);
-//   fhElectronPositionX       = new TH1F("hElectronPositionX","hElectronPositionX",500,-80. , 80.);
-//   fhNeutralHadronPositionX  = new TH1F("hNeutralHadronPositionX","hNeutralHadronPositionX",500,-80. , 80.);
-//   fhNeutralEMPositionX      = new TH1F("hNeutralEMPositionX","hNeutralEMPositionX",500,-80. , 80.);
-//   fhChargedHadronPositionX  = new TH1F("hChargedHadronPositionX","hChargedHadronPositionX",500,-80. , 80.);
-//   fhPhotonHadronPositionX   = new TH1F("hPhotonHadronPositionX","hPhotonHadronPositionX",500,-80. , 80.);
-//   fhPhotonPositionY         = new TH1F("hPhotonPositionY","hPhotonPositionY",   500,-80. , 80.);
-//   fhElectronPositionY       = new TH1F("hElectronPositionY","hElectronPositionY",500,-80. , 80.);
-//   fhNeutralHadronPositionY  = new TH1F("hNeutralHadronPositionY","hNeutralHadronPositionY",500,-80. , 80.);
-//   fhNeutralEMPositionY      = new TH1F("hNeutralEMPositionY","hNeutralEMPositionY",500,-80. , 80.);
-//   fhChargedHadronPositionY  = new TH1F("hChargedHadronPositionY","hChargedHadronPositionY",500,-80. , 80.);
-//   fhPhotonHadronPositionY   = new TH1F("hPhotonHadronPositionY","hPhotonHadronPositionY",500,-80. , 80.);
-
+  fhEmcDigit                = new TH1F("hEmcDigit",      "hEmcDigit",         1000,  0. ,  25.);
+  fhVetoDigit               = new TH1F("hVetoDigit",     "hVetoDigit",         500,  0. ,  3.e-5);
+  fhConvertorDigit          = new TH1F("hConvertorDigit","hConvertorDigit",    500,  0. ,  3.e-5);
+  fhEmcCluster              = new TH1F("hEmcCluster",    "hEmcCluster",       1000,  0. ,  30.);
+  fhVetoCluster             = new TH1F("hVetoCluster",   "hVetoCluster",       500,  0. ,  3.e-5);
+  fhConvertorCluster        = new TH1F("hConvertorCluster","hConvertorCluster",500,  0. ,  3.e-5);
+  fhConvertorEmc            = new TH2F("hConvertorEmc",  "hConvertorEmc",      200,  1. ,  3., 200, 0., 3.e-5);
 
 }
 //____________________________________________________________________________
@@ -772,81 +839,125 @@ void  AliPHOSAnalyze::BookResolutionHistograms()
 {
   // Books the histograms where the results of the Resolution analysis are stored
 
-  if(fhPhotonEnergy)
-    delete fhPhotonEnergy ;
-  if(fhPhotonAllEnergy)
-    delete fhPhotonAllEnergy ;
-  if(fhElectronEnergy)
-    delete fhElectronEnergy ;
-  if(fhElectronAllEnergy)
-    delete fhElectronAllEnergy ;
-  if(fhNeutralHadronEnergy)
-    delete fhNeutralHadronEnergy ;
-  if(fhNeutralEMEnergy)
-    delete fhNeutralEMEnergy ;
-  if(fhNeutralEMAllEnergy)
-    delete fhNeutralEMAllEnergy ;
-  if(fhChargedHadronEnergy)
-    delete fhChargedHadronEnergy ;
-  if(fhPhotonHadronEnergy)
-    delete fhPhotonHadronEnergy ;
-  if(fhPhotonPosition)
-    delete fhPhotonPosition ;
-  if(fhPhotonAllPosition)
-    delete fhPhotonAllPosition ;
-  if(fhElectronPosition)
-    delete fhElectronPosition ;
-  if(fhElectronAllPosition)
-    delete fhElectronAllPosition ;
-  if(fhNeutralHadronPosition)
-    delete fhNeutralHadronPosition ;
-  if(fhNeutralEMPosition)
-    delete fhNeutralEMPosition ;
-  if(fhNeutralEMAllPosition)
-    delete fhNeutralEMAllPosition ;
-  if(fhChargedHadronPosition)
-    delete fhChargedHadronPosition ;
-  if(fhPhotonHadronPosition)
-    delete fhPhotonHadronPosition ;
-
-  fhPhotonEnergy            = new TH2F("hPhotonEnergy",  "hPhotonEnergy",              100, 0., 5., 100, 0., 5.);
-  fhPhotonAllEnergy         = new TH2F("hPhotonAllEnergy",  "hPhotonAllEnergy",        100, 0., 5., 100, 0., 5.);
-  fhElectronEnergy          = new TH2F("hElectronEnergy","hElectronEnergy",            100, 0., 5., 100, 0., 5.);
-  fhElectronAllEnergy       = new TH2F("hElectronAllEnergy","hElectronAllEnergy",      100, 0., 5., 100, 0., 5.);
-  fhNeutralHadronEnergy     = new TH2F("hNeutralHadronEnergy", "hNeutralHadronEnergy", 100, 0., 5., 100, 0., 5.);
-  fhNeutralEMEnergy         = new TH2F("hNeutralEMEnergy", "hNeutralEMEnergy",         100, 0., 5., 100, 0., 5.);
-  fhNeutralEMAllEnergy      = new TH2F("hNeutralEMAllEnergy", "hNeutralEMAllEnergy",   100, 0., 5., 100, 0., 5.);
-  fhChargedHadronEnergy     = new TH2F("hChargedHadronEnergy", "hChargedHadronEnergy", 100, 0., 5., 100, 0., 5.);
-  fhPhotonHadronEnergy      = new TH2F("hPhotonHadronEnergy","hPhotonHadronEnergy",    100, 0., 5., 100, 0., 5.);
-  fhPhotonPosition          = new TH2F("hPhotonPosition","hPhotonPosition",                20, 0., 5., 100, 0., 5.);
-  fhPhotonAllPosition       = new TH2F("hPhotonAllPosition","hPhotonAllPosition",          20, 0., 5., 100, 0., 5.);
-  fhElectronPosition        = new TH2F("hElectronPosition","hElectronPosition",            20, 0., 5., 100, 0., 5.);
-  fhElectronAllPosition     = new TH2F("hElectronAllPosition","hElectronAllPosition",      20, 0., 5., 100, 0., 5.);
-  fhNeutralHadronPosition   = new TH2F("hNeutralHadronPosition","hNeutralHadronPosition",  20, 0., 5., 100, 0., 5.);
-  fhNeutralEMPosition       = new TH2F("hNeutralEMPosition","hNeutralEMPosition",          20, 0., 5., 100, 0., 5.);
-  fhNeutralEMAllPosition    = new TH2F("hNeutralEMAllPosition","hNeutralEMAllPosition",    20, 0., 5., 100, 0., 5.);
-  fhChargedHadronPosition   = new TH2F("hChargedHadronPosition","hChargedHadronPosition",  20, 0., 5., 100, 0., 5.);
-  fhPhotonHadronPosition    = new TH2F("hPhotonHadronPosition","hPhotonHadronPosition",    20, 0., 5., 100, 0., 5.);
-  fhPhotonReg = new TH1F("hPhotonReg","hPhotonReg", 20, 0., 5.);
-  fhAllReg    = new TH1F("hAllReg", "hAllReg",  20, 0., 5.);
-  fhNReg      = new TH1F("hNReg", "hNReg",  20, 0., 5.);
-  fhNBarReg   = new TH1F("hNBarReg", "hNBarReg",  20, 0., 5.);
-  fhChargedReg= new TH1F("hChargedReg", "hChargedReg",  20, 0., 5.);
-
-  fhPhotonEM = new TH1F("hPhotonEM","hPhotonEM", 20, 0., 5.);
-  fhAllEM    = new TH1F("hAllEM", "hAllEM",  20, 0., 5.);
-  fhNEM      = new TH1F("hNEM", "hNEM",  20, 0., 5.);
-  fhNBarEM   = new TH1F("hNBarEM", "hNBarEM",  20, 0., 5.);
-  fhChargedEM= new TH1F("hChargedEM", "hChargedEM",  20, 0., 5.);
-
-  fhPrimary= new TH1F("hPrimary", "hPrimary",  20, 0., 5.);
-
-  fhPhotPhot = new TH1F("hPhotPhot","hPhotPhot", 20, 0., 5.);   //Photon registered as photon
-  fhPhotElec = new TH1F("hPhotElec","hPhotElec", 20, 0., 5.);   //Photon registered as Electron
-  fhPhotNeuH = new TH1F("hPhotNeuH","hPhotNeuH", 20, 0., 5.);   //Photon registered as Neutral Hadron
-  fhPhotNuEM = new TH1F("hPhotNuEM","hPhotNuEM", 20, 0., 5.);   //Photon registered as Neutral EM
-  fhPhotChHa = new TH1F("hPhotChHa","hPhotChHa", 20, 0., 5.);   //Photon registered as Charged Hadron
-  fhPhotGaHa = new TH1F("hPhotGaHa","hPhotGaHa", 20, 0., 5.);   //Photon registered as Gamma-Hadron
+//   if(fhAllEnergy)
+//     delete fhAllEnergy ;
+//   if(fhPhotEnergy)
+//     delete fhPhotEnergy ;
+//   if(fhEMEnergy)
+//     delete fhEMEnergy ;
+//   if(fhPPSDEnergy)
+//     delete fhPPSDEnergy ;
+
+
+  fhAllEnergy  = new TH2F("hAllEnergy",  "Energy of any RP with primary photon",100, 0., 5., 100, 0., 5.);
+  fhPhotEnergy = new TH2F("hPhotEnergy", "Energy of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
+  fhEMEnergy   = new TH2F("hEMEnergy",   "Energy of EM with primary photon",    100, 0., 5., 100, 0., 5.);
+  fhPPSDEnergy = new TH2F("hPPSDEnergy", "Energy of PPSD with primary photon",  100, 0., 5., 100, 0., 5.);
+
+//   if(fhAllPosition)
+//     delete fhAllPosition ;
+//   if(fhPhotPosition)
+//     delete fhPhotPosition ;
+//   if(fhEMPosition)
+//     delete fhEMPosition ;
+//   if(fhPPSDPosition)
+//     delete fhPPSDPosition ;
+
+
+  fhAllPosition  = new TH2F("hAllPosition",  "Position of any RP with primary photon",100, 0., 5., 100, 0., 5.);
+  fhPhotPosition = new TH2F("hPhotPosition", "Position of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
+  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.);
+
+//   if(fhAllReg)
+//     delete fhAllReg ;
+//   if(fhPhotReg)
+//     delete fhPhotReg ;
+//   if(fhNReg)
+//     delete fhNReg ;
+//   if(fhNBarReg)
+//     delete fhNBarReg ;
+//   if(fhChargedReg)
+//     delete fhChargedReg ;
+  
+  fhAllReg    = new TH1F("hAllReg",    "All primaries registered as photon",  100, 0., 5.);
+  fhPhotReg   = new TH1F("hPhotReg",   "Photon registered as photon",         100, 0., 5.);
+  fhNReg      = new TH1F("hNReg",      "N registered as photon",              100, 0., 5.);
+  fhNBarReg   = new TH1F("hNBarReg",   "NBar registered as photon",           100, 0., 5.);
+  fhChargedReg= new TH1F("hChargedReg", "Charged hadron registered as photon",100, 0., 5.);
+  
+//   if(fhAllEM)
+//     delete fhAllEM ;
+//   if(fhPhotEM)
+//     delete fhPhotEM ;
+//   if(fhNEM)
+//     delete fhNEM ;
+//   if(fhNBarEM)
+//     delete fhNBarEM ;
+//   if(fhChargedEM)
+//     delete fhChargedEM ;
+  
+  fhAllEM    = new TH1F("hAllEM",    "All primary registered as EM",100, 0., 5.);
+  fhPhotEM   = new TH1F("hPhotEM",   "Photon registered as EM", 100, 0., 5.);
+  fhNEM      = new TH1F("hNEM",      "N registered as EM",      100, 0., 5.);
+  fhNBarEM   = new TH1F("hNBarEM",   "NBar registered as EM",   100, 0., 5.);
+  fhChargedEM= new TH1F("hChargedEM","Charged registered as EM",100, 0., 5.);
+
+//   if(fhAllPPSD)
+//     delete fhAllPPSD ;
+//   if(fhPhotPPSD)
+//     delete fhPhotPPSD ;
+//   if(fhNPPSD)
+//     delete fhNPPSD ;
+//   if(fhNBarPPSD)
+//     delete fhNBarPPSD ;
+//   if(fhChargedPPSD)
+//     delete fhChargedPPSD ;
+  
+  fhAllPPSD    = new TH1F("hAllPPSD",    "All primary registered as PPSD",100, 0., 5.);
+  fhPhotPPSD   = new TH1F("hPhotPPSD",   "Photon registered as PPSD", 100, 0., 5.);
+  fhNPPSD      = new TH1F("hNPPSD",      "N registered as PPSD",      100, 0., 5.);
+  fhNBarPPSD   = new TH1F("hNBarPPSD",   "NBar registered as PPSD",   100, 0., 5.);
+  fhChargedPPSD= new TH1F("hChargedPPSD","Charged registered as PPSD",100, 0., 5.);
+  
+//   if(fhPrimary)
+//     delete fhPrimary ;
+  fhPrimary= new TH1F("hPrimary", "hPrimary",  100, 0., 5.);
+
+//   if(fhAllRP)
+//     delete fhAllRP ;
+//   if(fhVeto)
+//     delete fhVeto ;
+//   if(fhShape)
+//     delete fhShape ;
+//   if(fhPPSD)
+//     delete fhPPSD ;
+
+  fhAllRP = new TH1F("hAllRP","All Reconstructed particles",  100, 0., 5.);
+  fhVeto  = new TH1F("hVeto", "All uncharged particles",      100, 0., 5.);
+  fhShape = new TH1F("hShape","All particles with EM shaower",100, 0., 5.);
+  fhPPSD  = new TH1F("hPPSD", "All PPSD photon particles",    100, 0., 5.);
+
+
+//   if(fhPhotPhot)
+//     delete fhPhotPhot ;
+//   if(fhPhotElec)
+//     delete fhPhotElec ;
+//   if(fhPhotNeuH)
+//     delete fhPhotNeuH ;
+//   if(fhPhotNuEM)
+//     delete fhPhotNuEM ;
+//   if(fhPhotChHa)
+//     delete fhPhotChHa ;
+//   if(fhPhotGaHa)
+//     delete fhPhotGaHa ;
+
+  fhPhotPhot = new TH1F("hPhotPhot","hPhotPhot", 100, 0., 5.);   //Photon registered as photon
+  fhPhotElec = new TH1F("hPhotElec","hPhotElec", 100, 0., 5.);   //Photon registered as Electron
+  fhPhotNeuH = new TH1F("hPhotNeuH","hPhotNeuH", 100, 0., 5.);   //Photon registered as Neutral Hadron
+  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
 
 
 }
@@ -1039,8 +1150,8 @@ void AliPHOSAnalyze::DisplayRecParticles()
     cout << "DisplayRecParticles > Analyze an event first ... (y/n) " ; 
     Text_t answer[1] ; 
     cin >> answer ; cout << answer ; 
-    if ( answer == "y" ) 
-      AnalyzeOneEvent() ;
+//     if ( answer == "y" ) 
+//       AnalyzeOneEvent() ;
   } 
     if (fEvt != -999) {
       
@@ -1118,8 +1229,8 @@ void AliPHOSAnalyze::DisplayRecPoints()
     cout << "DisplayRecPoints > Analyze an event first ... (y/n) " ; 
     Text_t answer[1] ; 
     cin >> answer ; cout << answer ; 
-    if ( answer == "y" ) 
-      AnalyzeOneEvent() ;
+//     if ( answer == "y" ) 
+//       AnalyzeOneEvent() ;
   } 
     if (fEvt != -999) {
       
@@ -1278,8 +1389,8 @@ void AliPHOSAnalyze::DisplayTrackSegments()
     cout << "DisplayTrackSegments > Analyze an event first ... (y/n) " ; 
     Text_t answer[1] ; 
     cin >> answer ; cout << answer ; 
-    if ( answer == "y" ) 
-      AnalyzeOneEvent() ;
+//     if ( answer == "y" ) 
+//       AnalyzeOneEvent() ;
   } 
     if (fEvt != -999) {
 
@@ -1349,7 +1460,7 @@ Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
   return  fRootFile->IsOpen() ; 
 }
 //____________________________________________________________________________
-void AliPHOSAnalyze::SavingHistograms()
+void AliPHOSAnalyze::SaveHistograms()
 {
   // Saves the histograms in a root file named "name.analyzed" 
 
@@ -1357,127 +1468,63 @@ void AliPHOSAnalyze::SavingHistograms()
   sprintf(outputname,"%s.analyzed",fRootFile->GetName());
   TFile output(outputname,"RECREATE");
   output.cd();
-//   if (fhEmcDigit )         
-//     fhEmcDigit->Write()  ;
-//   if (fhVetoDigit )  
-//     fhVetoDigit->Write()  ;
-//   if (fhConvertorDigit ) 
-//     fhConvertorDigit->Write()   ;
-//   if (fhEmcCluster   )
-//     fhEmcCluster->Write()   ;
-//   if (fhVetoCluster ) 
-//     fhVetoCluster->Write()   ;
-//   if (fhConvertorCluster )
-//     fhConvertorCluster->Write()  ;
-//   if (fhConvertorEmc ) 
-//     fhConvertorEmc->Write()  ;
-//   if (fhPhotonEnergy)    
-//     fhPhotonEnergy->Write() ;
-//   if (fhPhotonPositionX)  
-//     fhPhotonPositionX->Write() ;
-//   if (fhPhotonPositionY)  
-//     fhPhotonPositionX->Write() ;
-//   if (fhElectronEnergy)  
-//     fhElectronEnergy->Write() ;
-//   if (fhElectronPositionX)
-//     fhElectronPositionX->Write() ;
-//   if (fhElectronPositionY) 
-//     fhElectronPositionX->Write() ;
-//   if (fhNeutralHadronEnergy) 
-//     fhNeutralHadronEnergy->Write() ;
-//   if (fhNeutralHadronPositionX)
-//     fhNeutralHadronPositionX->Write() ;
-//   if (fhNeutralHadronPositionY) 
-//     fhNeutralHadronPositionX->Write() ;
-//   if (fhNeutralEMEnergy)   
-//     fhNeutralEMEnergy->Write() ;
-//   if (fhNeutralEMPositionX)
-//     fhNeutralEMPositionX->Write() ;
-//   if (fhNeutralEMPositionY) 
-//     fhNeutralEMPositionX->Write() ;
-//   if (fhChargedHadronEnergy) 
-//     fhChargedHadronEnergy->Write() ;
-//   if (fhChargedHadronPositionX) 
-//     fhChargedHadronPositionX->Write() ;
-//   if (fhChargedHadronPositionY)
-//     fhChargedHadronPositionX->Write() ;
-//   if (fhPhotonHadronEnergy) 
-//     fhPhotonHadronEnergy->Write() ;
-//   if (fhPhotonHadronPositionX) 
-//     fhPhotonHadronPositionX->Write() ;
-//   if (fhPhotonHadronPositionY)
-//     fhPhotonHadronPositionX->Write() ;
-
-  output.Write();
-  output.Close();
-}
-//____________________________________________________________________________
-void AliPHOSAnalyze::SaveResolutionHistograms()
-{
-  // Saves the histograms in a root file named "name.analyzed" 
 
-  Text_t outputname[80] ;
-  sprintf(outputname,"%s.analyzed",fRootFile->GetName());
-  TFile output(outputname,"RECREATE");
-  output.cd();
-  if (fhPhotonEnergy)    
-    fhPhotonEnergy->Write() ;
-  if (fhPhotonAllEnergy)    
-    fhPhotonAllEnergy->Write() ;
-  if (fhPhotonPosition)  
-    fhPhotonPosition->Write() ;
-  if (fhPhotonAllPosition)  
-    fhPhotonAllPosition->Write() ;
-  if (fhElectronEnergy)  
-    fhElectronEnergy->Write() ;
-  if (fhElectronAllEnergy)  
-    fhElectronAllEnergy->Write() ;
-  if (fhElectronPosition)
-    fhElectronPosition->Write() ;
-  if (fhElectronAllPosition)
-    fhElectronAllPosition->Write() ;
-  if (fhNeutralHadronEnergy) 
-    fhNeutralHadronEnergy->Write() ;
-  if (fhNeutralHadronPosition)
-    fhNeutralHadronPosition->Write() ;
-  if (fhNeutralEMEnergy)   
-    fhNeutralEMEnergy->Write() ;
-  if (fhNeutralEMAllEnergy)   
-    fhNeutralEMAllEnergy->Write() ;
-  if (fhNeutralEMPosition)
-    fhNeutralEMPosition->Write() ;
-  if (fhNeutralEMAllPosition)
-    fhNeutralEMAllPosition->Write() ;
-  if (fhChargedHadronEnergy) 
-    fhChargedHadronEnergy->Write() ;
-  if (fhChargedHadronPosition) 
-    fhChargedHadronPosition->Write() ;
-  if (fhPhotonHadronEnergy) 
-    fhPhotonHadronEnergy->Write() ;
-  if (fhPhotonHadronPosition) 
-    fhPhotonHadronPosition->Write() ;
-  if (fhPhotonReg) 
-    fhPhotonReg->Write() ;
+  if (fhAllEnergy)    
+    fhAllEnergy->Write() ;
+  if (fhPhotEnergy)    
+    fhPhotEnergy->Write() ;
+  if(fhEMEnergy)
+    fhEMEnergy->Write()  ;
+  if(fhPPSDEnergy)
+    fhPPSDEnergy->Write() ;
+  if(fhAllPosition)
+    fhAllPosition->Write() ;
+  if(fhPhotPosition)
+    fhPhotPosition->Write() ;
+  if(fhEMPosition)
+    fhEMPosition->Write() ;
+  if(fhPPSDPosition)
+    fhPPSDPosition->Write() ;
   if (fhAllReg) 
     fhAllReg->Write() ;
+  if (fhPhotReg) 
+    fhPhotReg->Write() ;
   if(fhNReg)
     fhNReg->Write() ;
   if(fhNBarReg)
     fhNBarReg->Write() ;
   if(fhChargedReg)
     fhChargedReg->Write() ;
-  if (fhPhotonEM) 
-    fhPhotonEM->Write() ;
   if (fhAllEM) 
     fhAllEM->Write() ;
+  if (fhPhotEM) 
+    fhPhotEM->Write() ;
   if(fhNEM)
     fhNEM->Write() ;
   if(fhNBarEM)
     fhNBarEM->Write() ;
   if(fhChargedEM)
     fhChargedEM->Write() ;
+  if (fhAllPPSD) 
+    fhAllPPSD->Write() ;
+  if (fhPhotPPSD) 
+    fhPhotPPSD->Write() ;
+  if(fhNPPSD)
+    fhNPPSD->Write() ;
+  if(fhNBarPPSD)
+    fhNBarPPSD->Write() ;
+  if(fhChargedPPSD)
+    fhChargedPPSD->Write() ;
   if(fhPrimary)
     fhPrimary->Write() ;
+  if(fhAllRP)
+    fhAllRP->Write()  ;
+  if(fhVeto)
+    fhVeto->Write()  ;
+  if(fhShape)
+    fhShape->Write()  ;
+  if(fhPPSD)
+    fhPPSD->Write()  ;
   if(fhPhotPhot)
     fhPhotPhot->Write() ;
   if(fhPhotElec)
@@ -1492,10 +1539,66 @@ void AliPHOSAnalyze::SaveResolutionHistograms()
     fhPhotChHa->Write() ;
   if(fhPhotGaHa)
     fhPhotGaHa->Write() ;
-
+  if(fhEnergyCorrelations)
+    fhEnergyCorrelations->Write() ;
+  
   output.Write();
   output.Close();
 }
+//____________________________________________________________________________
+Float_t AliPHOSAnalyze::CorrectEnergy(Float_t ERecPart)
+{
+  return ERecPart/0.8783 ;
+}
+
+//____________________________________________________________________________
+void AliPHOSAnalyze::ResetHistograms()
+{
+   fhEnergyCorrelations = 0 ;     //Energy correlations between Eloss in Convertor and PPSD(2)
+
+   fhEmcDigit = 0 ;               // Histo of digit energies in the Emc 
+   fhVetoDigit = 0 ;              // Histo of digit energies in the Veto 
+   fhConvertorDigit = 0 ;         // Histo of digit energies in the Convertor
+   fhEmcCluster = 0 ;             // Histo of Cluster energies in Emc
+   fhVetoCluster = 0 ;            // Histo of Cluster energies in Veto
+   fhConvertorCluster = 0 ;       // Histo of Cluster energies in Convertor
+   fhConvertorEmc = 0 ;           // 2d Convertor versus Emc energies
+
+   fhAllEnergy = 0 ;       
+   fhPhotEnergy = 0 ;        // Total spectrum of detected photons
+   fhEMEnergy = 0 ;         // Spectrum of detected electrons with electron primary
+   fhPPSDEnergy = 0 ;
+   fhAllPosition = 0 ; 
+   fhPhotPosition = 0 ; 
+   fhEMPosition = 0 ; 
+   fhPPSDPosition = 0 ; 
+
+   fhPhotReg = 0 ;          
+   fhAllReg = 0 ;          
+   fhNReg = 0 ;          
+   fhNBarReg = 0 ;          
+   fhChargedReg = 0 ;          
+   fhPhotEM = 0 ;          
+   fhAllEM = 0 ;          
+   fhNEM = 0 ;          
+   fhNBarEM = 0 ;          
+   fhChargedEM = 0 ;          
+   fhPhotPPSD = 0 ;          
+   fhAllPPSD = 0 ;          
+   fhNPPSD = 0 ;          
+   fhNBarPPSD = 0 ;          
+   fhChargedPPSD = 0 ;          
+
+   fhPrimary = 0 ;          
+
+   fhPhotPhot = 0 ;
+   fhPhotElec = 0 ;
+   fhPhotNeuH = 0 ;
+   fhPhotNuEM = 0 ; 
+   fhPhotChHa = 0 ;
+   fhPhotGaHa = 0 ;
+
 
+}