]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOSAnalyze.cxx
Variable format fixed (e.g., fx -> fX)
[u/mrichter/AliRoot.git] / PHOS / AliPHOSAnalyze.cxx
index 232eada7c91293115085bdf617d48d46de9169f4..5226cb4c2588d3a9f996268c0e8f241c2b599b69 100644 (file)
@@ -74,12 +74,20 @@ AliPHOSAnalyze::AliPHOSAnalyze(Text_t * name)
     cout << " AliPHOSAnalyze > Error opening " << name << endl ; 
   }
   else { 
-    gAlice = (AliRun*) fRootFile->Get("gAlice");
-    fPHOS  = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
-    fGeom  = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() ) ;
-   
-    fEvt = -999 ; 
+      //========== Get AliRun object from file 
+      gAlice = (AliRun*) fRootFile->Get("gAlice") ;
+
+      //=========== Get the PHOS object and associated geometry from the file      
+      fPHOS  = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
+      fGeom  = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
+      //========== Initializes the Index to Object converter
+      fObjGetter = AliPHOSIndexToObject::GetInstance(fPHOS) ; 
+      //========== Current event number 
+      fEvt = -999 ; 
+
   }
+  ResetHistograms() ;
 }
 
 //____________________________________________________________________________
@@ -104,120 +112,138 @@ 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) ; 
-  
-  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
-
+void AliPHOSAnalyze::ActivePPSD(Int_t Nevents=1){
   
-    // =========== Finish
-
-    cout << "AnalyzeOneEvent > event # " << fEvt << " processed" << endl ;   
-  } // ok
-  else
-    cout << "AnalyzeOneEvent > filed to process event # " << evt << endl ;   
+  fhEnergyCorrelations  = new TH2F("hEnergyCorrelations","hEnergyCorrelations",40,  0., 0.15, 30, 0., 3.e-5);
+  //========== Create the Clusterizer
+  fClu = new AliPHOSClusterizerv1() ; 
+  fClu->SetEmcEnergyThreshold(0.05) ; 
+  fClu->SetEmcClusteringThreshold(0.20) ; 
+  fClu->SetPpsdEnergyThreshold    (0.0000002) ; 
+  fClu->SetPpsdClusteringThreshold(0.0000001) ; 
+  fClu->SetLocalMaxCut(0.03) ;
+  fClu->SetCalibrationParameters(0., 0.00000001) ; 
+
+  Int_t ievent;
   
-  ts.Stop() ;   cout << "CPU  time = " << ts.CpuTime()   << endl ; 
-  cout << "Real time = " << ts.RealTime()  << endl ; 
+  for ( ievent=0; ievent<Nevents; ievent++)
+    {  
+      
+      //========== Event Number>         
+      if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. ) 
+       cout <<  "AnalyzeResolutions > " << "Event is " << ievent << endl ;  
+      
+      //=========== Connects the various Tree's for evt
+      gAlice->GetEvent(ievent);
+
+      //=========== Gets the Kine TTree
+      gAlice->TreeK()->GetEvent(0) ;
+            
+      //=========== Get the Digit Tree
+      gAlice->TreeD()->GetEvent(0) ;
+      
+      //========== Creating branches ===================================       
+      AliPHOSRecPoint::RecPointsList ** EmcRecPoints =  fPHOS->EmcRecPoints() ;
+      gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP", EmcRecPoints ) ;
+      
+      AliPHOSRecPoint::RecPointsList ** PpsdRecPoints = fPHOS->PpsdRecPoints() ;
+      gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", PpsdRecPoints ) ;
+      
+      AliPHOSTrackSegment::TrackSegmentsList **  TrackSegmentsList = fPHOS->TrackSegments() ;
+      if( (*TrackSegmentsList) )
+       (*TrackSegmentsList)->Clear() ;
+      gAlice->TreeR()->SetBranchAddress( "PHOSTS", TrackSegmentsList ) ;
+      
+      AliPHOSRecParticle::RecParticlesList ** RecParticleList  = fPHOS->RecParticles() ;
+      if( (*RecParticleList) )
+       (*RecParticleList)->Clear() ;
+      gAlice->TreeR()->SetBranchAddress( "PHOSRP", RecParticleList ) ;
+      
+      
+      //=========== Gets the Reconstraction TTree
+      gAlice->TreeR()->GetEvent(0) ;
+            
+      AliPHOSPpsdRecPoint * RecPoint ;
+      Int_t relid[4] ; 
+      TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
+      while( ( RecPoint = (AliPHOSPpsdRecPoint *)nextRP() ) )
+       {
+         if(!(RecPoint->GetUp()) ) {
+           AliPHOSDigit *digit ;
+           Int_t iDigit ;
+           for(iDigit = 0; iDigit < fPHOS->Digits()->GetEntries(); iDigit++) 
+             {
+               digit = (AliPHOSDigit *) fPHOS->Digits()->At(iDigit) ;
+               fGeom->AbsToRelNumbering(digit->GetId(), relid) ;    
+                if((relid[2]==1)&&(relid[3]==1)&&(relid[0]==RecPoint->GetPHOSMod())){
+                 Float_t ConvertorEnergy = fClu->Calibrate(digit->GetAmp()) ;
+                 fhEnergyCorrelations->Fill(ConvertorEnergy,RecPoint->GetTotalEnergy() );  
+                 break ; 
+               } 
+             }
+           break ;
+         }
+       }
+    }
+  SaveHistograms() ;
+  fhEnergyCorrelations->Draw("BOX") ;
 }
 
 //____________________________________________________________________________
  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 ;  
   else
     {
-      //========== Get AliRun object from file 
-      gAlice = (AliRun*) fRootFile->Get("gAlice") ;
-      //=========== Get the PHOS object and associated geometry from the file      
-      fPHOS  = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
-      fGeom  = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
       //========== 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 ;  
+
          //=========== 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() ) )
@@ -225,43 +251,44 @@ void AliPHOSAnalyze::AnalyzeOneEvent(Int_t evt)
              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()));
+             }
+           }
+         
+
+         //=========== Cluster in module
+         TIter nextEmc(*fPHOS->EmcRecPoints()  ) ;
+         while((emc = (AliPHOSEmcRecPoint *)nextEmc())) 
+           {
+             if ( emc->GetPHOSMod() == module )
                {  
-                 if (relid[1]<17) fhVetoDigit->Fill(fClu->Calibrate(digit->GetAmp())); 
-                 if (relid[1]>16) fhConvertorDigit->Fill(fClu->Calibrate(digit->GetAmp()));
+                 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()) ;
                }
            }
-         //=========== 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()) ;
-//             }
-//         }
          //========== TRackSegments in the event
-         TIter nextRecParticle(fPHOS->RecParticles() ) ; 
+         TIter nextRecParticle(*fPHOS->RecParticles() ) ; 
          while((recparticle = (AliPHOSRecParticle *)nextRecParticle())) 
            {
              if ( recparticle->GetPHOSTrackSegment()->GetPHOSMod() == module )
@@ -273,210 +300,415 @@ 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 ;                 
-//                 }
                }
            }
-         // Deleting fClu, fTrs, fPID et fRec
-         fClu->Delete();
-         fTrs->Delete();
-         fPID->Delete();
-         fRec->Delete();
-
        }   // endfor
-      SavingHistograms();
+      SaveHistograms();
     }       // endif
 }           // endfunction
 
 //____________________________________________________________________________
- void AliPHOSAnalyze::AnalyzeResolutions(Int_t Nevents )    
-{
-  // analyzes Nevents events and calculate Energy and Position resolution as well as
-  // probaility of correct indentifiing of the incident particle
-
-  if ( fRootFile == 0 ) 
-    cout << "AnalyzeResolutions > " << "Root File not openned" << endl ;  
-  else
-    {
-      //========== Get AliRun object from file 
-      gAlice = (AliRun*) fRootFile->Get("gAlice") ;
-      //=========== Get the PHOS object and associated geometry from the file      
-      fPHOS  = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
-      fGeom  = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
-
-      //========== Initializes the Index to Object converter
-      fObjGetter = AliPHOSIndexToObject::GetInstance(fPHOS) ; 
-
-      //========== Booking Histograms
-      cout << "AnalyzeResolutions > " << "Booking Histograms" << endl ; 
-      BookResolutionHistograms();
-
-      Int_t ievent;
-      for ( ievent=0; ievent<Nevents; ievent++)
-       {  
-          if (ievent==0)  cout << "AnalyzeResolutions > " << "Starting Analyzing " << endl ; 
-
+ void AliPHOSAnalyze::Reconstruct(Int_t Nevents )    
+{     
+  Int_t ievent ;   
+  for ( ievent=0; ievent<Nevents; ievent++)
+    {  
+      if (ievent==0) 
+       {
+         cout << "Analyze > Starting Reconstructing " << endl ; 
          //========== Create the Clusterizer
          fClu = new AliPHOSClusterizerv1() ; 
          fClu->SetEmcEnergyThreshold(0.05) ; 
-         fClu->SetEmcClusteringThreshold(0.50) ; 
+         fClu->SetEmcClusteringThreshold(0.20) ; 
          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 ) ; 
-           
+         
          //========== Creates the Reconstructioner  
          fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ; 
+         //  fRec -> SetDebugReconstruction(kTRUE);     
 
-         //========== Event Number>         
-         if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. ) 
-           cout <<  "AnalyzeResolutions > " << "Event is " << ievent << endl ;  
+       }
+      
+      //========== Event Number>         
+      if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. ) 
+       cout <<  "Analyze > Event is " << ievent << endl ;  
+      
+      //=========== Connects the various Tree's for evt
+      gAlice->GetEvent(ievent);
 
-         //=========== Connects the various Tree's for evt
-         gAlice->GetEvent(ievent);
+      //=========== Gets the Digit TTree
+      gAlice->TreeD()->GetEvent(0) ;
 
-          //=========== Get the Digit Tree
-         gAlice->TreeD()->GetEvent(0) ;
+      //=========== Do the reconstruction
+      fPHOS->Reconstruction(fRec);
+    }
 
-         //=========== Do the reconstruction
-         fPHOS->Reconstruction(fRec);
+    fClu->Delete();
+    fTrs->Delete();
+    fPID->Delete();
+    fRec->Delete();
 
-          //========== 
-         TClonesArray * RecParticleList     = new  TClonesArray("AliPHOSRecParticle", 2000) ;
-         TBranch * RecBranch = gAlice->TreeR()->GetBranch("PHOSRP");
-         RecBranch-> SetAddress(&RecParticleList);
+}
+//-------------------------------------------------------------------------------------
+
+//   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))) ;
+      
+//       cout << "# of Primaries after add " << PrimaryList->GetEntriesFast() <<endl;
 
-         //=========== Gets the Reconstraction TTree
-         gAlice->TreeR()->GetEvent(0) ;
+//       cout << "Digits before add " << DigitsList->GetEntries() << endl ;
+//       cout << "Digits to add " <<  NDigits[ievent+1]-  NDigits[ievent]<< endl ;
+      
+      //=========== 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 ;
 
-         //=========== Gets the Kine TTree
-         gAlice->TreeK()->GetEvent(0) ;
 
-         //=========== Gets the list of Primari Particles
-         TClonesArray * PrimaryList  = gAlice->Particles();    
+//____________________________________________________________________________
+ void AliPHOSAnalyze::AnalyzeResolutions(Int_t Nevents )    
+{
+  // analyzes Nevents events and calculate Energy and Position resolution as well as
+  // probaility of correct indentifiing of the incident particle
 
-         //========== TRackSegments in the event
-         TIter nextRecParticle(RecParticleList) ; 
-          AliPHOSRecParticle * RecParticle ;
+  //========== Booking Histograms
+  cout << "AnalyzeResolutions > " << "Booking Histograms" << endl ; 
+  BookResolutionHistograms();
+
+  Int_t Counter[9][5] ;     
+  Int_t i1,i2,TotalInd = 0 ;
+  for(i1 = 0; i1<9; i1++)
+    for(i2 = 0; i2<5; i2++)
+      Counter[i1][i2] = 0 ;
+  
+  Int_t TotalPrimary = 0 ;
+  Int_t TotalRecPart = 0 ;
+  Int_t TotalRPwithPrim = 0 ;
+  Int_t ievent;
+
+  cout << "Start Analysing"<< endl ;
+  for ( ievent=0; ievent<Nevents; ievent++)
+    {  
+      
+      //========== Event Number>         
+      if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. ) 
+       cout <<  "AnalyzeResolutions > " << "Event is " << ievent << endl ;  
+      
+      //=========== Connects the various Tree's for evt
+      gAlice->GetEvent(ievent);
+
+      
+
+            //=========== Gets the Kine TTree
+      gAlice->TreeK()->GetEvent(0) ;
+      
+      //=========== Gets the list of Primari Particles
+      TClonesArray * PrimaryList  = gAlice->Particles();     
 
-         while( (RecParticle = (AliPHOSRecParticle *) nextRecParticle()))
+      TParticle * Primary ;
+      Int_t iPrimary ;
+      for ( iPrimary = 0 ; iPrimary < PrimaryList->GetEntries() ; iPrimary++)
            {
-              Int_t ModuleNumberRec ;
-              Double_t RecX, RecZ ;
-             fGeom->ImpactOnEmc(RecParticle->Theta(), RecParticle->Phi(), ModuleNumberRec, RecX, RecZ) ;
+             Primary = (TParticle*)PrimaryList->UncheckedAt(iPrimary) ;
+             Int_t PrimaryType = Primary->GetPdgCode() ;
+             if( PrimaryType == 22 ) 
+               fhPrimary->Fill(Primary->Energy()) ;
+           } 
+      
+      //=========== Get the Digit Tree
+      gAlice->TreeD()->GetEvent(0) ;
+
+      //========== Creating branches ===================================       
+      AliPHOSRecPoint::RecPointsList ** EmcRecPoints =  fPHOS->EmcRecPoints() ;
+      gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP", EmcRecPoints ) ;
 
-              Double_t MinDistance = 10000 ;
-             Int_t ClosestPrimary = -1 ;
+      AliPHOSRecPoint::RecPointsList ** PpsdRecPoints = fPHOS->PpsdRecPoints() ;
+      gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", PpsdRecPoints ) ;
 
-             Int_t numberofprimaries ;
-             Int_t * listofprimaries  = RecParticle->GetPrimaries(numberofprimaries)  ;
+      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) ;
+      cout << ievent << " " << (*EmcRecPoints) << " " <<(*PpsdRecPoints) <<fPHOS->Digits()<< endl ; 
+      cout << "    " << " " << (*EmcRecPoints)->GetEntries() << " " <<(*PpsdRecPoints)->GetEntries() <<fPHOS->Digits()->GetEntries()<< endl ; 
 
-             Int_t index ;
-              TParticle * Primary ;
-             Double_t Distance = MinDistance ;
-             for ( index = 0 ; index < numberofprimaries ; index++)
+      AliPHOSRecParticle * RecParticle ;
+      Int_t iRecParticle ;
+      for(iRecParticle = 0; iRecParticle < (*RecParticleList)->GetEntries() ;iRecParticle++ )
+       {
+         RecParticle = (AliPHOSRecParticle *) (*RecParticleList)->At(iRecParticle) ;
+         
+         Int_t ModuleNumberRec ;
+         Double_t RecX, RecZ ;
+         fGeom->ImpactOnEmc(RecParticle->Theta(), RecParticle->Phi(), ModuleNumberRec, RecX, RecZ) ;
+         
+         Double_t MinDistance = 10000 ;
+         Int_t ClosestPrimary = -1 ;
+         
+         Int_t numberofprimaries ;
+         Int_t * listofprimaries  = RecParticle->GetPrimaries(numberofprimaries)  ;
+         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)
                {
-                 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] ;
-                   }
+                 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;
                }
 
-             if(ClosestPrimary >=0 )
+             switch(RecParticle->GetType())
                {
-                 switch(RecParticle->GetType())
-                   {
-                   case AliPHOSFastRecParticle::kGAMMA:
-                     fhPhotonEnergy->Fill(RecParticle->Energy(),((TParticle *) PrimaryList->At(ClosestPrimary))->Energy() ) ; 
-                     fhPhotonPosition->Fill(RecParticle->Energy(),Distance) ;
-                     break;
-                   case  AliPHOSFastRecParticle::kELECTRON:
-                     fhElectronEnergy->Fill(RecParticle->Energy(),((TParticle *)PrimaryList->At(ClosestPrimary))->Energy() ) ; 
-                     fhElectronPosition->Fill(RecParticle->Energy(),Distance ) ;
-                     break;
-                   case  AliPHOSFastRecParticle::kNEUTRALHA:
-                     fhNeutralHadronEnergy->Fill( RecParticle->Energy(),((TParticle *)PrimaryList->At(ClosestPrimary))->Energy()) ; 
-                     fhNeutralHadronPosition->Fill(RecParticle->Energy(),Distance  ) ;
-                     break ;
-                   case  AliPHOSFastRecParticle::kNEUTRALEM:
-                     fhNeutralEMEnergy->Fill( RecParticle->Energy(),((TParticle *)PrimaryList->At(ClosestPrimary))->Energy() ) ; 
-                     fhNeutralEMPosition->Fill(RecParticle->Energy(),Distance ) ;
-                     break ;
-                   case  AliPHOSFastRecParticle::kCHARGEDHA:
-                     fhChargedHadronEnergy->Fill(RecParticle->Energy(),((TParticle *)PrimaryList->At(ClosestPrimary))->Energy() ) ; 
-                     fhChargedHadronPosition->Fill(RecParticle->Energy(),Distance ) ;
-                     break ;
-                   case  AliPHOSFastRecParticle::kGAMMAHA:
-                     fhPhotonHadronEnergy->Fill( RecParticle->Energy(),((TParticle *)PrimaryList->At(ClosestPrimary))->Energy()) ; 
-                     fhPhotonHadronPosition->Fill(RecParticle->Energy(),Distance ) ;
-                     break ;                 
-                   }
+               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(PrimaryType == 2112){ //neutron
+                   fhNReg->Fill(RecParticle->Energy() ) ;
+                   fhNEM->Fill(RecParticle->Energy() ) ;
+                 }
+                 
+                 if(PrimaryType == -2112){ //neutron ~
+                   fhNBarReg->Fill(RecParticle->Energy() ) ;
+                   fhNBarEM->Fill(RecParticle->Energy() ) ;
+
+                 }
+                 if(PrimaryCode == 2){
+                   fhChargedReg->Fill(RecParticle->Energy() ) ;
+                   fhChargedEM->Fill(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;
                }
            }
-
-         // Deleting fClu, fTrs, fPID et fRec
-         fClu->Delete();
-         fTrs->Delete();
-         fPID->Delete();
-         fRec->Delete();
-         PrimaryList->Delete() ;
-
-       }   // endfor
-      SaveResolutionHistograms();
-    }       // endif
+       }  
+    }   // endfor
+  SaveHistograms();
+  cout << "Resolutions: Analyzed " << Nevents << " event(s)" << endl ;
+  cout << "Resolutions: Total primary       " << TotalPrimary << endl ;
+  cout << "Resoluitons: Total reconstracted " << TotalRecPart << endl ;
+  cout << "TotalReconstructed with Primarie " << TotalRPwithPrim << endl ;
+  cout << "                        Primary:   Photon   Electron   Ch. Hadr.  Neutr. Hadr  Kaons" << endl ; 
+  cout << "             Detected as photon       " << Counter[0][0] << "          " << Counter[0][1] << "          " << Counter[0][2] << "          " <<Counter[0][3] << "          " << Counter[0][4] << endl ;
+  cout << "           Detected as electron       " << Counter[1][0] << "          " << Counter[1][1] << "          " << Counter[1][2] << "          " <<Counter[1][3] << "          " << Counter[1][4] << endl ; 
+  cout << "     Detected as neutral hadron       " << Counter[2][0] << "          " << Counter[2][1] << "          " << Counter[2][2] << "          " <<Counter[2][3] << "          " << Counter[2][4] << endl ;
+  cout << "         Detected as neutral EM       " << Counter[3][0] << "          " << Counter[3][1] << "          " << Counter[3][2] << "          " <<Counter[3][3] << "          " << Counter[3][4] << endl ;
+  cout << "     Detected as charged hadron       " << Counter[4][0] << "          " << Counter[4][1] << "          " << Counter[4][2] << "          " <<Counter[4][3] << "          " << Counter[4][4] << endl ;
+  cout << "       Detected as gamma-hadron       " << Counter[5][0] << "          " << Counter[5][1] << "          " << Counter[5][2] << "          " <<Counter[5][3] << "          " << Counter[5][4] << endl ;
+  cout << "          Detected as Absurd EM       " << Counter[6][0] << "          " << Counter[6][1] << "          " << Counter[6][2] << "          " <<Counter[6][3] << "          " << Counter[6][4] << endl ;
+  cout << "      Detected as absurd hadron       " << Counter[7][0] << "          " << Counter[7][1] << "          " << Counter[7][2] << "          " <<Counter[7][3] << "          " << Counter[7][4] << endl ;
+  cout << "          Detected as undefined       " << Counter[8][0] << "          " << Counter[8][1] << "          " << Counter[8][2] << "          " <<Counter[8][3] << "          " << Counter[8][4] << endl ;
+      
+      for(i1 = 0; i1<9; i1++)
+       for(i2 = 0; i2<5; i2++)
+         TotalInd+=Counter[i1][i2] ;
+      cout << "Indentified particles            " << TotalInd << endl ;
+      
 }           // endfunction
 
 
@@ -485,47 +717,21 @@ void  AliPHOSAnalyze::BookingHistograms()
 {
   // Books the histograms where the results of the analysis are stored (to be changed)
 
-  if (fhEmcDigit )         
-    delete fhEmcDigit  ;
-  if (fhVetoDigit )     
-    delete fhVetoDigit  ;
-  if (fhConvertorDigit ) 
-    delete fhConvertorDigit   ;
-  if (fhEmcCluster   )  
-    delete  fhEmcCluster   ;
-  if (fhVetoCluster )     
-    delete fhVetoCluster   ;
-  if (fhConvertorCluster )
-    delete fhConvertorCluster  ;
-  if (fhConvertorEmc )    
-    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.);
-
+  delete fhEmcDigit  ;
+  delete fhVetoDigit  ;
+  delete fhConvertorDigit   ;
+  delete  fhEmcCluster   ;
+  delete fhVetoCluster   ;
+  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);
 
 }
 //____________________________________________________________________________
@@ -535,41 +741,118 @@ void  AliPHOSAnalyze::BookResolutionHistograms()
 
   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",                100, 0., 5., 100, 0., 5.);
-  fhElectronPosition        = new TH2F("hElectronPosition","hElectronPosition",            100, 0., 5., 100, 0., 5.);
-  fhNeutralHadronPosition   = new TH2F("hNeutralHadronPosition","hNeutralHadronPosition",  100, 0., 5., 100, 0., 5.);
-  fhNeutralEMPosition       = new TH2F("hNeutralEMPosition","hNeutralEMPosition",          100, 0., 5., 100, 0., 5.);
-  fhChargedHadronPosition   = new TH2F("hChargedHadronPosition","hChargedHadronPosition",  100, 0., 5., 100, 0., 5.);
-  fhPhotonHadronPosition    = new TH2F("hPhotonHadronPosition","hPhotonHadronPosition",    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.);
+
+  if(fhPhotonReg)
+    delete fhPhotonReg ;
+  if(fhAllReg)
+    delete fhAllReg ;
+  if(fhNReg)
+    delete fhNReg ;
+  if(fhNReg)
+    delete fhNReg ;
+  if(fhNReg)
+    delete fhNReg ;
+  
+  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.);
+  
+  if(fhPhotonEM)
+    delete fhPhotonEM ;
+  if(fhAllEM)
+    delete fhAllEM ;
+  if(fhNEM)
+    delete fhNEM ;
+  if(fhNBarEM)
+    delete fhNBarEM ;
+  if(fhChargedEM)
+    delete fhChargedEM ;
+  
+  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.);
+  
+  if(fhPrimary)
+    delete fhPrimary ;
+  fhPrimary= new TH1F("hPrimary", "hPrimary",  20, 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", 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
+
 
 }
 //____________________________________________________________________________
@@ -607,11 +890,6 @@ Bool_t AliPHOSAnalyze::Init(Int_t evt)
   
   if ( ok ) {
     
-
-    //========== Initializes the Index to Object converter
-
-    fObjGetter = AliPHOSIndexToObject::GetInstance(fPHOS) ; 
-
     //========== Create the Clusterizer
 
     fClu =  new AliPHOSClusterizerv1() ; 
@@ -766,8 +1044,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) {
       
@@ -791,7 +1069,7 @@ void AliPHOSAnalyze::DisplayRecParticles()
       Text_t canvasname[80] ; 
       sprintf(canvasname, "Reconstructed particles in PHOSmodule # %d", module) ;
       TCanvas * rparticlecanvas = new TCanvas("RparticleCanvas", canvasname, 650, 500) ; 
-      AliPHOSRecParticle::RecParticlesList * rpl = fPHOS->RecParticles() ; 
+      AliPHOSRecParticle::RecParticlesList * rpl = *fPHOS->RecParticles() ; 
       Int_t nRecParticles = rpl->GetEntries() ; 
       Int_t nRecParticlesInModule = 0 ; 
       TIter nextRecPart(rpl) ; 
@@ -845,8 +1123,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) {
       
@@ -915,7 +1193,7 @@ void AliPHOSAnalyze::DisplayRecPoints()
       //=========== Cluster in module
 
       //      TClonesArray * emcRP = fPHOS->EmcClusters() ; 
-      TObjArray * emcRP = fPHOS->EmcRecPoints() ; 
+      TObjArray * emcRP = *(fPHOS->EmcRecPoints()) ; 
       
       etot = 0.; 
       Int_t totalnClusters = 0 ; 
@@ -950,7 +1228,7 @@ void AliPHOSAnalyze::DisplayRecPoints()
       //=========== Cluster in module PPSD Down
 
       //      TClonesArray * ppsdRP = fPHOS->PpsdClusters() ;
-      TObjArray * ppsdRP = fPHOS->PpsdRecPoints() ;
+      TObjArray * ppsdRP = *(fPHOS->PpsdRecPoints() );
  
       etot = 0.; 
       TIter nextPpsd(ppsdRP) ;
@@ -972,7 +1250,7 @@ void AliPHOSAnalyze::DisplayRecPoints()
 
       //=========== Cluster in module PPSD Up
   
-      ppsdRP = fPHOS->PpsdRecPoints() ;
+      ppsdRP = *(fPHOS->PpsdRecPoints()) ;
      
       etot = 0.; 
       TIter nextPpsdUp(ppsdRP) ;
@@ -1005,8 +1283,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) {
 
@@ -1042,7 +1320,7 @@ void AliPHOSAnalyze::DisplayTrackSegments()
       TCanvas * trackcanvas = new TCanvas("TrackSegmentCanvas", canvasname, 650, 500) ; 
       histotrack->Draw() ; 
 
-      AliPHOSTrackSegment::TrackSegmentsList * trsegl = fPHOS->TrackSegments() ;
+      AliPHOSTrackSegment::TrackSegmentsList * trsegl = *(fPHOS->TrackSegments()) ;
       AliPHOSTrackSegment * trseg ;
  
       Int_t nTrackSegments = trsegl->GetEntries() ;
@@ -1076,14 +1354,14 @@ Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
   return  fRootFile->IsOpen() ; 
 }
 //____________________________________________________________________________
-void AliPHOSAnalyze::SavingHistograms()
-{
-  // 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();
+// void AliPHOSAnalyze::SavingHistograms()
+// {
+//   // 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 (fhEmcDigit )         
 //     fhEmcDigit->Write()  ;
 //   if (fhVetoDigit )  
@@ -1135,11 +1413,11 @@ void AliPHOSAnalyze::SavingHistograms()
 //   if (fhPhotonHadronPositionY)
 //     fhPhotonHadronPositionX->Write() ;
 
-  output.Write();
-  output.Close();
-}
+//   output.Write();
+//   output.Close();
+// }
 //____________________________________________________________________________
-void AliPHOSAnalyze::SaveResolutionHistograms()
+void AliPHOSAnalyze::SaveHistograms()
 {
   // Saves the histograms in a root file named "name.analyzed" 
 
@@ -1147,22 +1425,35 @@ void AliPHOSAnalyze::SaveResolutionHistograms()
   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) 
@@ -1171,10 +1462,105 @@ void AliPHOSAnalyze::SaveResolutionHistograms()
     fhPhotonHadronEnergy->Write() ;
   if (fhPhotonHadronPosition) 
     fhPhotonHadronPosition->Write() ;
-
+  if (fhPhotonReg) 
+    fhPhotonReg->Write() ;
+  if (fhAllReg) 
+    fhAllReg->Write() ;
+  if(fhNReg)
+    fhNReg->Write() ;
+  if(fhNBarReg)
+    fhNBarReg->Write() ;
+  if(fhChargedReg)
+    fhChargedReg->Write() ;
+  if (fhPhotonEM) 
+    fhPhotonEM->Write() ;
+  if (fhAllEM) 
+    fhAllEM->Write() ;
+  if(fhNEM)
+    fhNEM->Write() ;
+  if(fhNBarEM)
+    fhNBarEM->Write() ;
+  if(fhChargedEM)
+    fhChargedEM->Write() ;
+  if(fhPrimary)
+    fhPrimary->Write() ;
+  if(fhPhotPhot)
+    fhPhotPhot->Write() ;
+  if(fhPhotElec)
+    fhPhotElec->Write() ;
+  if(fhPhotNeuH)
+    fhPhotNeuH->Write() ;
+  if(fhPhotNuEM)
+    fhPhotNuEM->Write() ;
+  if(fhPhotNuEM)
+    fhPhotNuEM->Write() ;
+  if(fhPhotChHa)
+    fhPhotChHa->Write() ;
+  if(fhPhotGaHa)
+    fhPhotGaHa->Write() ;
+  if(fhEnergyCorrelations)
+    fhEnergyCorrelations->Write() ;
+  
   output.Write();
   output.Close();
 }
+//____________________________________________________________________________
+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
+
+   fhPhotonEnergy = 0 ;           // Spectrum of detected photons with photon primary
+   fhPhotonAllEnergy = 0 ;        // Total spectrum of detected photons
+   fhElectronEnergy = 0 ;         // Spectrum of detected electrons with electron primary
+   fhElectronAllEnergy = 0 ;      // Total spectrum of detected electrons
+   fhNeutralHadronEnergy = 0 ;    // Spectrum of detected neutral hadron
+   fhNeutralEMEnergy = 0 ;        // Spectrum of detected neutral EM with EM primary
+   fhNeutralEMAllEnergy = 0 ;     // Spectrum of detected neutral EM
+   fhChargedHadronEnergy = 0 ;    // Spectrum of detected charged
+   fhPhotonHadronEnergy = 0 ;     // Spectrum of detected Photon-Hadron
+   fhPhotonPosition = 0 ;        // Position Resolution of  photons with photon primary
+   fhPhotonAllPosition = 0 ;     // Position Resolution of  photons
+   fhElectronPosition = 0 ;      // Position Resolution of electrons with electron primary
+   fhElectronAllPosition = 0 ;   // Position Resolution of electrons
+   fhNeutralHadronPosition = 0 ; // Position Resolution of neutral hadron
+   fhNeutralEMPosition = 0 ;     // Position Resolution of neutral EM with EM primary
+   fhNeutralEMAllPosition = 0 ;  // Position Resolution of neutral EM
+   fhChargedHadronPosition = 0 ; // Position Resolution of charged
+   fhPhotonHadronPosition = 0 ;  // Position Resolution of Photon-Hadron
+   fhPhotonPositionY = 0 ;        // Y distribution of detected photons
+   fhElectronPositionY = 0 ;      // Y distribution of detected electrons
+   fhNeutralHadronPositionY = 0 ; // Y distribution of detected neutral hadron
+   fhNeutralEMPositionY = 0 ;     // Y distribution of detected neutral EM
+   fhChargedHadronPositionY = 0 ; // Y distribution of detected charged
+   fhPhotonHadronPositionY = 0 ;  // Y distribution of detected Photon-Hadron
+   fhPhotonReg = 0 ;          
+   fhAllReg = 0 ;          
+   fhNReg = 0 ;          
+   fhNBarReg = 0 ;          
+   fhChargedReg = 0 ;          
+   fhPhotonEM = 0 ;          
+   fhAllEM = 0 ;          
+   fhNEM = 0 ;          
+   fhNBarEM = 0 ;          
+   fhChargedEM = 0 ;          
+   fhPrimary = 0 ;          
+
+   fhPhotPhot = 0 ;
+   fhPhotElec = 0 ;
+   fhPhotNeuH = 0 ;
+   fhPhotNuEM = 0 ; 
+   fhPhotChHa = 0 ;
+   fhPhotGaHa = 0 ;
+
 
+}