]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
Dimitri just makes it work
authorschutz <schutz@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 30 Oct 2000 15:36:56 +0000 (15:36 +0000)
committerschutz <schutz@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 30 Oct 2000 15:36:56 +0000 (15:36 +0000)
PHOS/AliPHOSAnalyze.cxx
PHOS/AliPHOSAnalyze.h
PHOS/AliPHOSEmcRecPoint.cxx
PHOS/AliPHOSPID.cxx
PHOS/AliPHOSPID.h
PHOS/AliPHOSPIDv1.cxx
PHOS/AliPHOSPIDv1.h
PHOS/AliPHOSRecPoint.cxx
PHOS/AliPHOSTrackSegmentMakerv1.cxx
PHOS/AliPHOSTrackSegmentMakerv1.h
PHOS/AliPPSDGeometry.cxx

index d42b41c162f0775b37620b7a9c20b6ba5b121b0f..b8c9a84a07ce9358985f659d8afa857798ed4635 100644 (file)
@@ -34,7 +34,6 @@
 #include "TTree.h"
 #include "TMath.h"
 #include "TCanvas.h" 
-#include "TStyle.h" 
 
 // --- Standard library ---
 
@@ -53,7 +52,6 @@
 #include "AliPHOSTrackSegment.h"
 #include "AliPHOSRecParticle.h"
 #include "AliPHOSIndexToObject.h"
-#include "AliPHOSCPV.h"
 
 ClassImp(AliPHOSAnalyze)
 
@@ -89,7 +87,10 @@ AliPHOSAnalyze::AliPHOSAnalyze(Text_t * name)
       fEvt = -999 ; 
 
   }
-  fDebugLevel = 0;
+  fClu = 0 ;
+  fPID = 0 ;
+  fTrs = 0 ;
+  fRec = 0 ;
   ResetHistograms() ;
 }
 
@@ -115,12 +116,23 @@ AliPHOSAnalyze::~AliPHOSAnalyze()
 
   if (fRootFile->IsOpen() ) 
     fRootFile->Close() ; 
-  if(fRootFile) {delete fRootFile ; fRootFile=0 ;}
-  if(fPHOS)     {delete fPHOS     ; fPHOS    =0 ;}
-  if(fClu)      {delete fClu      ; fClu     =0 ;}
-  if(fPID)      {delete fPID      ; fPID     =0 ;}
-  if(fRec)      {delete fRec      ; fRec     =0 ;}
-  if(fTrs)      {delete fTrs      ; fTrs     =0 ;}
+  if(fRootFile)
+    delete fRootFile ;  
+
+  if(fPHOS)
+    delete fPHOS ; 
+
+  if(fClu)
+    delete fClu ; 
+
+  if(fPID)
+    delete fPID ; 
+
+  if(fRec)
+    delete fRec ; 
+
+  if(fTrs)
+    delete fTrs ; 
 
 }
 
@@ -130,11 +142,11 @@ void AliPHOSAnalyze::ActivePPSD(Int_t Nevents=1){
   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->SetEmcEnergyThreshold(0.01) ; 
   fClu->SetEmcClusteringThreshold(0.20) ; 
   fClu->SetPpsdEnergyThreshold    (0.0000002) ; 
   fClu->SetPpsdClusteringThreshold(0.0000001) ; 
-  fClu->SetLocalMaxCut(0.03) ;
+  fClu->SetLocalMaxCut(0.02) ;
   fClu->SetCalibrationParameters(0., 0.00000001) ; 
 
   Int_t ievent;
@@ -151,7 +163,7 @@ void AliPHOSAnalyze::ActivePPSD(Int_t Nevents=1){
 
       //=========== Gets the Kine TTree
       gAlice->TreeK()->GetEvent(0) ;
-      
+            
       //=========== Get the Digit Tree
       gAlice->TreeD()->GetEvent(0) ;
       
@@ -202,6 +214,7 @@ void AliPHOSAnalyze::ActivePPSD(Int_t Nevents=1){
   fhEnergyCorrelations->Draw("BOX") ;
 }
 
+
 //____________________________________________________________________________
 void AliPHOSAnalyze::AnalyzeManyEvents(Int_t Nevents, Int_t module)    
 {
@@ -300,83 +313,21 @@ void AliPHOSAnalyze::AnalyzeManyEvents(Int_t Nevents, Int_t module)
 }           // endfunction
 
 //____________________________________________________________________________
-void AliPHOSAnalyze::ReconstructCPV(Int_t Nevents )    
+ void AliPHOSAnalyze::Reconstruct(Int_t Nevents,Int_t FirstEvent )    
 {     
-
-  // Perform reconstruction of EMC and CPV (GPS2 or IHEP) for <Nevents> events
-  // Yuri Kharlov. 19 October 2000
-
   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.20) ; 
-      fClu->SetLocalMaxCut(0.03) ;
-      if      (strcmp(fGeom->GetName(),"GPS2") == 0) {
-       fClu->SetPpsdEnergyThreshold    (0.0000002) ; 
-       fClu->SetPpsdClusteringThreshold(0.0000001) ; 
-      }
-      else if (strcmp(fGeom->GetName(),"IHEP") == 0) {
-       fClu->SetLocalMaxCutCPV(0.03) ;
-       fClu->SetLogWeightCutCPV(4.0) ;
-       fClu->SetPpsdEnergyThreshold    (0.09) ;
-      }
-      fClu->SetCalibrationParameters(0., 0.00000001) ; 
-      
-      //========== Creates the track segment maker
-      fTrs = new AliPHOSTrackSegmentMakerv1()  ;
-      
-      //========== Creates the particle identifier for GPS2 only
-      if      (strcmp(fGeom->GetName(),"GPS2") == 0) {
-       fPID = new AliPHOSPIDv1() ;
-       fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ; 
-      }          
-      
-      //========== Creates the Reconstructioner
-      fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ; 
-      if (fDebugLevel != 0) fRec -> SetDebugReconstruction(kTRUE);     
-    }
-      
-    if (fDebugLevel != 0 ||
-       (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
-      cout <<  "======= Analyze ======> Event " << ievent+1 << endl ;
-    
-    //=========== Connects the various Tree's for evt
-    gAlice->GetEvent(ievent);
-    
-    //=========== Gets the Digit TTree
-    gAlice->TreeD()->GetEvent(0) ;
-    
-    //=========== Do the reconstruction
-    fPHOS->Reconstruction(fRec);
-  }
-
-  if(fClu)      {delete fClu      ; fClu     =0 ;}
-  if(fPID)      {delete fPID      ; fPID     =0 ;}
-  if(fRec)      {delete fRec      ; fRec     =0 ;}
-  if(fTrs)      {delete fTrs      ; fTrs     =0 ;}
-  
-}
-//-------------------------------------------------------------------------------------
-
-void AliPHOSAnalyze::Reconstruct(Int_t Nevents )    
-{     
-  Int_t ievent ;   
-  for ( ievent=0; ievent<Nevents; ievent++)
+  for ( ievent=FirstEvent; ievent<Nevents; ievent++)
     {  
-      if (ievent==0
+      if (ievent==FirstEvent
        {
          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
@@ -386,16 +337,18 @@ void AliPHOSAnalyze::Reconstruct(Int_t Nevents )
          //========== 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 ((ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
-      cout <<  "======= Analyze ======> Event " << ievent+1 << 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);
@@ -407,268 +360,131 @@ void AliPHOSAnalyze::Reconstruct(Int_t Nevents )
       fPHOS->Reconstruction(fRec);
     }
 
-  if(fClu)      {delete fClu      ; fClu     =0 ;}
-  if(fPID)      {delete fPID      ; fPID     =0 ;}
-  if(fRec)      {delete fRec      ; fRec     =0 ;}
-  if(fTrs)      {delete fTrs      ; fTrs     =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))) ;
-      
-//       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 ;
-      
-      //=========== 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 ;
-
-
-//____________________________________________________________________________
-void AliPHOSAnalyze::ReadAndPrintCPV(Int_t Nevents)
-{
-  //
-  // Read and print generated and reconstructed hits in CPV
-  // Author: Yuri Kharlov
-  // 12 October 2000
-  //
-
-  cout << "Start CPV Analysis"<< endl ;
-  for ( Int_t ievent=0; ievent<Nevents; ievent++) {  
-      
-    //========== Event Number>         
-    cout << endl <<  "==== ReadAndPrintCPV ====> Event is " << ievent+1 << endl ;
-    
-    //=========== Connects the various Tree's for evt
-    gAlice->GetEvent(ievent);
-
-    //=========== Get the Hits Tree
-    gAlice->ResetHits();
-    gAlice->TreeH()->GetEvent(0);
-    
-    //========== Creating branches ===================================
-    AliPHOSRecPoint::RecPointsList ** EmcRecPoints = fPHOS->EmcRecPoints() ;
-    gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , EmcRecPoints  ) ;
-    
-    AliPHOSRecPoint::RecPointsList ** CpvRecPoints = fPHOS->PpsdRecPoints() ;
-    gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", CpvRecPoints ) ;
-
-    //=========== Gets the Reconstruction TTree
-    gAlice->TreeR()->GetEvent(0) ;
-
-    // Read and print CPV hits
-
-    TClonesArray *CPVhits;
-    for (Int_t iModule=0; iModule < fGeom->GetNModules(); iModule++) {
-      CPVModule   cpvModule = fPHOS->GetCPVModule(iModule);
-      CPVhits   = cpvModule.Hits();
-      Int_t nCPVhits  = CPVhits->GetEntriesFast();
-      for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
-       CPVHit        *cpvHit = (CPVHit*)CPVhits->UncheckedAt(ihit);
-       TLorentzVector p      = cpvHit->GetMomentum();
-       Float_t        xgen   = cpvHit->GetX();
-       Float_t        zgen   = cpvHit->GetY();
-       Int_t          ipart  = cpvHit->GetIpart();
-       printf("CPV hit in module %d: ",iModule+1);
-       printf(" p = (%f, %f, %f, %f) GeV,\n",
-              p.Px(),p.Py(),p.Pz(),p.Energy());
-       printf("               xy = (%8.4f, %8.4f) cm, ipart = %d\n",
-              xgen,zgen,ipart);
-      }
-    }
-
-    // Read and print CPV reconstructed points
+    fClu->Delete();
+    fClu=0 ;
+    fTrs->Delete();
+    fTrs = 0 ;
+    fPID->Delete();
+    fPID = 0 ;
+    fRec->Delete();
+    fRec = 0 ;
 
-    TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
-    AliPHOSPpsdRecPoint *cpvRecPoint ;
-    while( ( cpvRecPoint = (AliPHOSPpsdRecPoint *)nextRP() ) ) {
-      TVector3  locpos;
-      cpvRecPoint->GetLocalPosition(locpos);
-      Int_t PHOSModule = cpvRecPoint->GetPHOSMod();
-      printf("CPV recpoint in module %d: (X,Y,Z) = (%f,%f,%f) cm\n",
-            PHOSModule,locpos.X(),locpos.Y(),locpos.Z());
-    }
-  }
 }
-
 //____________________________________________________________________________
-void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
+ void AliPHOSAnalyze::InvariantMass(Int_t Nevents )    
 {
-  //
-  // Analyzes CPV characteristics
-  // Author: Yuri Kharlov
-  // 9 October 2000
-  //
-
-  // Book histograms
-
-  TH1F *hDx   = new TH1F("hDx"  ,"CPV x-resolution@reconstruction",100,-5. , 5.);
-  TH1F *hDz   = new TH1F("hDz"  ,"CPV z-resolution@reconstruction",100,-5. , 5.);
-  TH1S *hNrp  = new TH1S("hNrp" ,"CPV rec.point multiplicity",      21,-0.5,20.5);
-
-  cout << "Start CPV Analysis"<< endl ;
-  for ( Int_t ievent=0; ievent<Nevents; ievent++) {  
-      
-    //========== Event Number>         
-    if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
-      cout << endl <<  "==== AnalyzeCPV ====> Event is " << ievent+1 << endl ;
-    
-    //=========== Connects the various Tree's for evt
-    gAlice->GetEvent(ievent);
-
-    //=========== Get the Hits Tree
-    gAlice->ResetHits();
-    gAlice->TreeH()->GetEvent(0);
-    
-    //========== Creating branches ===================================
-    AliPHOSRecPoint::RecPointsList ** EmcRecPoints = fPHOS->EmcRecPoints() ;
-    gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , EmcRecPoints  ) ;
+  // 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 ;
     
-    AliPHOSRecPoint::RecPointsList ** CpvRecPoints = fPHOS->PpsdRecPoints() ;
-    gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", CpvRecPoints ) ;
-
-    //=========== Gets the Reconstruction TTree
-    gAlice->TreeR()->GetEvent(0) ;
-
-    TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
-    AliPHOSEmcRecPoint *cpvRecPoint ;
-    CPVModule           cpvModule;
-    TClonesArray        *CPVhits;
-    while( ( cpvRecPoint = (AliPHOSEmcRecPoint *)nextRP() ) ) {
-      TVector3  locpos;
-      cpvRecPoint->GetLocalPosition(locpos);
-      Int_t PHOSModule = cpvRecPoint->GetPHOSMod();
-      Int_t rpMult     = cpvRecPoint->GetDigitsMultiplicity();
-      Float_t xrec  = locpos.X();
-      Float_t zrec  = locpos.Z();
-      Float_t dxmin = 1.e+10;
-      Float_t dzmin = 1.e+10;
-
-      cpvModule = fPHOS->GetCPVModule(PHOSModule-1);
-      CPVhits   = cpvModule.Hits();
-      Int_t nCPVhits  = CPVhits->GetEntriesFast();
-      for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
-       CPVHit        *cpvHit = (CPVHit*)CPVhits->UncheckedAt(ihit);
-       Float_t        xgen   = cpvHit->GetX();
-       Float_t        zgen   = cpvHit->GetY();
-       if ( TMath::Abs(xgen-xrec) < TMath::Abs(dxmin) ) dxmin = xgen-xrec;
-       if ( TMath::Abs(zgen-zrec) < TMath::Abs(dzmin) ) dzmin = zgen-zrec;
-      }
-      cpvModule.Clear();
-      hDx  ->Fill(dxmin);
-      hDz  ->Fill(dzmin);
-      hNrp ->Fill(rpMult);
+    for ( ievent=0; ievent < NMixedEvents; ievent++){        
+      
+      Int_t AbsEventNumber = EventInMixedLoop*NMixedEvents + ievent ;
+      
+      //=========== Connects the various Tree's for evt
+      gAlice->GetEvent(AbsEventNumber);
+      
+      //=========== Get the Digit Tree
+      gAlice->TreeD()->GetEvent(0) ;
+      
+      //========== 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 ;  
     }
+    
 
-  }
-
-  // Save histograms
-
-  Text_t outputname[80] ;
-  sprintf(outputname,"%s.analyzed",fRootFile->GetName());
-  TFile output(outputname,"RECREATE");
-  output.cd();
-
-  hDx  ->Write() ;
-  hDz  ->Write() ;
-  hNrp ->Write() ;
-
-  // Plot histograms
-
-  TCanvas *CPVcanvas = new TCanvas("CPV","CPV analysis",20,20,300,900);
-  gStyle->SetOptStat(111111);
-  gStyle->SetOptFit(1);
-  gStyle->SetOptDate(1);
-  CPVcanvas->Divide(3,1);
+    //Now calculate invariant mass:
+    Int_t irp1,irp2 ;
+    Int_t NCurEvent = 0 ;
 
-  CPVcanvas->cd(1);
-  gPad->SetFillColor(10);
-  hNrp->SetFillColor(16);
-  hNrp->Draw();
+    for(irp1 = 0; irp1 < AllRecParticleList->GetEntries()-1; irp1++){
+      AliPHOSRecParticle * rp1 = (AliPHOSRecParticle *)AllRecParticleList->At(irp1) ;
 
-  CPVcanvas->cd(2);
-  gPad->SetFillColor(10);
-  hDx->SetFillColor(16);
-  hDx->Fit("gaus");
-  hDx->Draw();
+      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
 
-  CPVcanvas->cd(3);
-  gPad->SetFillColor(10);
-  hDz->SetFillColor(16);
-  hDz->Fit("gaus");
-  hDz->Draw();
 
-  CPVcanvas->Print("CPV.ps");
+    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();
 
 }
 
@@ -698,15 +514,13 @@ void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
     {  
       
       //========== 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
@@ -715,23 +529,31 @@ void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
       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() ;
@@ -741,25 +563,22 @@ void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
       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 ; 
-
+      
       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 ;
@@ -767,157 +586,206 @@ void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
          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
   SaveHistograms();
@@ -971,119 +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.);
-
-  if(fhPhotonReg)
-    delete fhPhotonReg ;
-  if(fhAllReg)
-    delete fhAllReg ;
-  if(fhNReg)
-    delete fhNReg ;
-  if(fhNReg)
-    delete fhNReg ;
-  if(fhNReg)
-    delete fhNReg ;
+//   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 ;
   
-  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.);
+  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(fhPhotonEM)
-    delete fhPhotonEM ;
-  if(fhAllEM)
-    delete fhAllEM ;
-  if(fhNEM)
-    delete fhNEM ;
-  if(fhNBarEM)
-    delete fhNBarEM ;
-  if(fhChargedEM)
-    delete fhChargedEM ;
+//   if(fhAllEM)
+//     delete fhAllEM ;
+//   if(fhPhotEM)
+//     delete fhPhotEM ;
+//   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.);
+  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 ;
   
-  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
+  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
 
 
 }
@@ -1150,8 +1024,7 @@ Bool_t AliPHOSAnalyze::Init(Int_t evt)
     //========== Creates the Reconstructioner  
     
     fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
-//      fRec -> SetDebugReconstruction(kFALSE);     
-    fRec -> SetDebugReconstruction(kTRUE);     
+    fRec -> SetDebugReconstruction(kFALSE);     
     
     //=========== Connect the various Tree's for evt
     
@@ -1587,69 +1460,6 @@ 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();
-//   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::SaveHistograms()
 {
   // Saves the histograms in a root file named "name.analyzed" 
@@ -1659,64 +1469,62 @@ void AliPHOSAnalyze::SaveHistograms()
   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)
@@ -1737,6 +1545,12 @@ void AliPHOSAnalyze::SaveHistograms()
   output.Write();
   output.Close();
 }
+//____________________________________________________________________________
+Float_t AliPHOSAnalyze::CorrectEnergy(Float_t ERecPart)
+{
+  return ERecPart/0.8783 ;
+}
+
 //____________________________________________________________________________
 void AliPHOSAnalyze::ResetHistograms()
 {
@@ -1750,40 +1564,31 @@ void AliPHOSAnalyze::ResetHistograms()
    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 ;          
+   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 ;          
-   fhPhotonEM = 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 ;
index 043c91522862697010645b8461a315582ab88c0c..796d22d3fe5252830a7196cf857a0783cd957623 100644 (file)
@@ -35,14 +35,13 @@ public:
 
   void ActivePPSD(Int_t Nevents) ;
   void AnalyzeManyEvents(Int_t Nevtents = 100, Int_t Module=0) ;  // analyzes many events   ;
-  void Reconstruct(Int_t Nevtents = 100) ;
-  void ReconstructCPV(Int_t Nevents = 1) ;  // reconstruct points in EMC and CPV
-  void AnalyzeResolutions(Int_t Nevtents) ; // analyzes Energy and Position resolutions   ;
-  void ReadAndPrintCPV(Int_t Nevents);      // Read & print generated and reconstructed hits in CPV
-  void AnalyzeCPV(Int_t Nevents);           // analyzes various CPV characteristics
+  void InvariantMass(Int_t Nevents = 100) ; 
+  void Reconstruct(Int_t Nevtents = 100,Int_t FirstEvent = 0) ;
+  void AnalyzeResolutions(Int_t Nevtents) ;  // analyzes Energy and Position resolutions   ;
   void BookingHistograms() ;                // booking histograms for the ManyEvent analysis ;
   void BookResolutionHistograms() ;         // booking histograms for the Resoluion analysis ;
   void Copy(TObject & obj) ;                // copies an analysis into an other one   
+  Float_t CorrectEnergy(Float_t ERecPart) ;   //Corrects reconstracted energy
   Bool_t Init(Int_t evt) ;                  // does various initialisations
   void DisplayKineEvent(Int_t evt = -999) ; // displays the Kine events in ALICE coordinate 
   void DisplayRecParticles() ;              // displays RecParticles in ALICE coordinate  
@@ -57,7 +56,6 @@ public:
     assert(0==1) ;
     return *this ; 
   }
-  void SetDebugLevel(Int_t flag) { fDebugLevel = flag; }
  
  private:
   
@@ -71,8 +69,6 @@ public:
   TFile * fRootFile ;                 // the root file that contains the data
   AliPHOSTrackSegmentMaker * fTrs ;   // a tracksegmentmaker ;
 
-  Int_t fDebugLevel;                  // debug level for analysis
-
   TH2F * fhEnergyCorrelations ;     //Energy correlations between Eloss in Convertor and PPSD(2)
 
 
@@ -84,41 +80,44 @@ public:
   TH1F * fhConvertorCluster ;       // Histo of Cluster energies in Convertor
   TH2F * fhConvertorEmc ;           // 2d Convertor versus Emc energies
 
-  TH2F * fhPhotonEnergy ;           // Spectrum of detected photons with photon primary
-  TH2F * fhPhotonAllEnergy ;        // Total spectrum of detected photons
-  TH2F * fhElectronEnergy ;         // Spectrum of detected electrons with electron primary
-  TH2F * fhElectronAllEnergy ;      // Total spectrum of detected electrons
-  TH2F * fhNeutralHadronEnergy ;    // Spectrum of detected neutral hadron
-  TH2F * fhNeutralEMEnergy ;        // Spectrum of detected neutral EM with EM primary
-  TH2F * fhNeutralEMAllEnergy ;     // Spectrum of detected neutral EM
-  TH2F * fhChargedHadronEnergy ;    // Spectrum of detected charged
-  TH2F * fhPhotonHadronEnergy ;     // Spectrum of detected Photon-Hadron
-  TH2F * fhPhotonPosition ;        // Position Resolution of  photons with photon primary
-  TH2F * fhPhotonAllPosition ;     // Position Resolution of  photons
-  TH2F * fhElectronPosition ;      // Position Resolution of electrons with electron primary
-  TH2F * fhElectronAllPosition ;   // Position Resolution of electrons
-  TH2F * fhNeutralHadronPosition ; // Position Resolution of neutral hadron
-  TH2F * fhNeutralEMPosition ;     // Position Resolution of neutral EM with EM primary
-  TH2F * fhNeutralEMAllPosition ;  // Position Resolution of neutral EM
-  TH2F * fhChargedHadronPosition ; // Position Resolution of charged
-  TH2F * fhPhotonHadronPosition ;  // Position Resolution of Photon-Hadron
+  TH2F * fhAllEnergy ;           // Spectrum of detected photons with photon primary
+  TH2F * fhPhotEnergy ;        // Total spectrum of detected photons
+  TH2F * fhEMEnergy ;        // Spectrum of detected neutral EM with EM primary
+  TH2F * fhPPSDEnergy ;
+
+  TH2F * fhAllPosition ;        // Position Resolution of  photons with photon primary
+  TH2F * fhPhotPosition ;     // Position Resolution of  photons
+  TH2F * fhEMPosition ;     // Position Resolution of neutral EM with EM primary
+  TH2F * fhPPSDPosition ;  // Position Resolution of neutral EM
+
   TH1F * fhPhotonPositionY ;        // Y distribution of detected photons
   TH1F * fhElectronPositionY ;      // Y distribution of detected electrons
   TH1F * fhNeutralHadronPositionY ; // Y distribution of detected neutral hadron
   TH1F * fhNeutralEMPositionY ;     // Y distribution of detected neutral EM
   TH1F * fhChargedHadronPositionY ; // Y distribution of detected charged
   TH1F * fhPhotonHadronPositionY ;  // Y distribution of detected Photon-Hadron
-  TH1F * fhPhotonReg ;          
+
+  TH1F * fhPhotReg ;          
   TH1F * fhAllReg ;          
   TH1F * fhNReg ;          
   TH1F * fhNBarReg ;          
   TH1F * fhChargedReg ;          
-  TH1F * fhPhotonEM ;          
+  TH1F * fhPhotEM ;          
   TH1F * fhAllEM ;          
   TH1F * fhNEM ;          
   TH1F * fhNBarEM ;          
   TH1F * fhChargedEM ;          
+  TH1F * fhPhotPPSD ;          
+  TH1F * fhAllPPSD ;          
+  TH1F * fhNPPSD ;          
+  TH1F * fhNBarPPSD ;          
+  TH1F * fhChargedPPSD ;          
+
   TH1F * fhPrimary ;          
+  TH1F * fhAllRP ;
+  TH1F * fhPPSD ;
+  TH1F * fhShape ;
+  TH1F * fhVeto ;
 
   TH1F * fhPhotPhot ;
   TH1F * fhPhotElec ;
index fdc4f892f7eda89e6cc03876bcbd5fa5e00f5a81..24cf5bc29ff7afb7055ca2cf071e25d222626456 100644 (file)
@@ -55,6 +55,7 @@ AliPHOSEmcRecPoint::AliPHOSEmcRecPoint(Float_t W0, Float_t LocMaxCut)
   fW0        = W0 ;          
   fLocMaxCut = LocMaxCut ; 
   fLocPos.SetX(1000000.)  ;      //Local position should be evaluated
+  
 }
 
 //____________________________________________________________________________
@@ -352,6 +353,19 @@ void  AliPHOSEmcRecPoint::GetElipsAxis(Float_t * lambda)
   dxz /= wtot ;
   dxz -= x * z ;
 
+//   //Apply correction due to non-perpendicular incidence
+//   Double_t CosX ;
+//   Double_t CosZ ;
+//   Double_t DistanceToIP= (Double_t ) ((AliPHOSGeometry *) fGeom)->GetIPtoCrystalSurface() ;
+  
+//   CosX = DistanceToIP/TMath::Sqrt(DistanceToIP*DistanceToIP+x*x) ;
+//   CosZ = DistanceToIP/TMath::Sqrt(DistanceToIP*DistanceToIP+z*z) ;
+
+//   dxx = dxx/(CosX*CosX) ;
+//   dzz = dzz/(CosZ*CosZ) ;
+//   dxz = dxz/(CosX*CosZ) ;
+
+
   lambda[0] =  0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
   if(lambda[0] > 0)
     lambda[0] = TMath::Sqrt(lambda[0]) ;
@@ -388,6 +402,7 @@ void AliPHOSEmcRecPoint::GetLocalPosition(TVector3 &LPos)
 
   Int_t iDigit;
 
+
   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
     digit = (AliPHOSDigit *) ( please->GimeDigit(fDigitsList[iDigit]) );
 
@@ -399,17 +414,11 @@ void AliPHOSEmcRecPoint::GetLocalPosition(TVector3 &LPos)
     x    += xi * w ;
     z    += zi * w ;
     wtot += w ;
-  }
 
-  if (wtot != 0) {
-    x /= wtot ;
-    z /= wtot ;
-  } else {
-    x = -1e6 ;
-    z = -1e6 ;
-    if (fMulDigit != 0) cout << "AliPHOSEMCRecPoint: too low log weight factor "
-                            << "to evaluate cluster's center\n";
   }
+
+  x /= wtot ;
+  z /= wtot ;
   fLocPos.SetX(x)  ;
   fLocPos.SetY(0.) ;
   fLocPos.SetZ(z)  ;
index 65456912e09fc5b6c81a11deef7e9ef2c9d7ef8e..c369710eac689f40456d30c5d0ed63033e274c8c 100644 (file)
@@ -40,6 +40,8 @@ ClassImp(AliPHOSPID)
 AliPHOSPID::AliPHOSPID()
 {
   // ctor
+  fGeom = AliPHOSGeometry::GetInstance() ;
+
 }
 
 //____________________________________________________________________________
index ccd709591d5dd1676405f4fc3538db2e07645541..ec875722ef78fd7086db5d8589068504de6a9ef3 100644 (file)
@@ -22,6 +22,7 @@
 
 #include "AliPHOSTrackSegment.h"
 #include "AliPHOSRecParticle.h"
+#include "AliPHOSGeometry.h"
 
 
 
@@ -34,8 +35,16 @@ public:
 
   virtual void MakeParticles(AliPHOSTrackSegment::TrackSegmentsList * trsl, 
                             AliPHOSRecParticle::RecParticlesList * rpl) {} ; 
+  virtual void Print(const char *){} ; 
   virtual void SetShowerProfileCuts(Float_t, Float_t, Float_t, Float_t) {} ; 
-  virtual void SetDispersionCutOff(Float_t ) {}    
+  virtual void SetDispersionCutOff(Float_t ) {} ;   
+  virtual void SetRelativeDistanceCut(Float_t ) {};
+
+ protected:
+  
+  AliPHOSGeometry      * fGeom ;           // pointer to PHOS geometry  
+
+
 
   ClassDef(AliPHOSPID,1)  // Particle Identifier algorithm (base class)
 
index dd16d3ae567b569477db028baa32347fd6db9fbd..8d95e39d3c904dc4e7715aef2d88db476b53a5f6 100644 (file)
 
 ClassImp( AliPHOSPIDv1) 
 
+//____________________________________________________________________________
+Float_t  AliPHOSPIDv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcclu,AliPHOSPpsdRecPoint * PpsdClu, Bool_t &toofar, Option_t *  Axis)
+{
+  // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
+  Float_t r;
+  TVector3 vecEmc ;
+  TVector3 vecPpsd ;
+  
+  emcclu->GetLocalPosition(vecEmc) ;
+  PpsdClu->GetLocalPosition(vecPpsd)  ; 
+  if(emcclu->GetPHOSMod() == PpsdClu->GetPHOSMod())
+    { 
+
+      // Correct to difference in CPV and EMC position due to different distance to center.
+      // we assume, that particle moves from center
+      Float_t dCPV = fGeom->GetIPtoOuterCoverDistance();
+      Float_t dEMC = fGeom->GetIPtoCrystalSurface() ;
+      dEMC         = dEMC / dCPV ;
+      vecPpsd = dEMC * vecPpsd  - vecEmc ; 
+      r = vecPpsd.Mag() ;
+      if (Axis == "X") r = vecPpsd.X();
+      if (Axis == "Y") r = vecPpsd.Y();
+      if (Axis == "Z") r = vecPpsd.Z();
+      if (Axis == "R") r = vecPpsd.Mag();
+      
+    } 
+  else 
+    {
+      toofar = kTRUE ;
+    }
+  return r ;
+}
+
+
+
+
 //____________________________________________________________________________
 void  AliPHOSPIDv1::MakeParticles(AliPHOSTrackSegment::TrackSegmentsList * trsl, 
                                  AliPHOSRecParticle::RecParticlesList * rpl)
@@ -47,6 +84,7 @@ void  AliPHOSPIDv1::MakeParticles(AliPHOSTrackSegment::TrackSegmentsList * trsl,
   AliPHOSTrackSegment * tracksegment ; 
   Int_t index = 0 ; 
   AliPHOSRecParticle * rp ; 
+  Bool_t tDistance;
   Int_t type ; 
   Int_t showerprofile;  // 0 narrow and 1 wide
   Int_t cpvdetector;  // 1 hit and 0 no hit
@@ -55,13 +93,23 @@ void  AliPHOSPIDv1::MakeParticles(AliPHOSTrackSegment::TrackSegmentsList * trsl,
   while ( (tracksegment = (AliPHOSTrackSegment *)next()) ) {
     new( (*rpl)[index] ) AliPHOSRecParticle(tracksegment) ;
     rp = (AliPHOSRecParticle *)rpl->At(index) ; 
-    AliPHOSEmcRecPoint * recp = tracksegment->GetEmcRecPoint() ; 
-    Float_t * lambda = new Float_t[2]; 
-    recp->GetElipsAxis(lambda) ; 
-
- // Looking at the lateral development of the shower
-    if ( ( lambda[0] > fLambda1m && lambda[0] < fLambda1M ) && // shower profile cut
-        ( lambda[1] > fLambda2m && lambda[1] < fLambda2M ) )         
+    AliPHOSEmcRecPoint * recp = tracksegment->GetEmcRecPoint() ;
+    AliPHOSPpsdRecPoint * rpcpv = tracksegment->GetPpsdUpRecPoint() ;
+    AliPHOSPpsdRecPoint * rppc  = tracksegment->GetPpsdLowRecPoint() ;
+//     Float_t * lambda = new Float_t[2]; 
+//     recp->GetElipsAxis(lambda) ; 
+
+//     // Looking at the lateral development of the shower
+//     if ( ( lambda[0] > fLambda1m && lambda[0] < fLambda1M ) && // shower profile cut
+//      ( lambda[1] > fLambda2m && lambda[1] < fLambda2M ) )         
+//       //    Float_t R ;
+//       //R=(lambda[0]-1.386)*(lambda[0]-1.386)+1.707*1.707*(lambda[1]-1.008)*(lambda[1]-1.008) ;
+//       //if(R<0.35*0.35)
+
+    Float_t Dispersion;
+    Dispersion = recp->GetDispersion();
+    if (Dispersion < fCutOnDispersion)
       showerprofile = 0 ;   // NARROW PROFILE   
     else      
       showerprofile = 1 ;// WIDE PROFILE
@@ -70,13 +118,13 @@ void  AliPHOSPIDv1::MakeParticles(AliPHOSTrackSegment::TrackSegmentsList * trsl,
     if( tracksegment->GetPpsdLowRecPoint() == 0 )   
       pcdetector = 0 ;  // No hit
     else      
-      pcdetector = 1 ;  // hit
+      if (GetDistanceInPHOSPlane(recp, rppc, tDistance, "R")< fCutOnRelativeDistance)  pcdetector = 1 ;  // hit
   
     // Looking at the photon conversion detector
     if( tracksegment->GetPpsdUpRecPoint() == 0 )
       cpvdetector = 0 ;  // No hit
     else  
-      cpvdetector = 1 ;  // Hit
+      if (GetDistanceInPHOSPlane(recp, rpcpv, tDistance, "R")< fCutOnRelativeDistance) cpvdetector = 1 ;  // Hit
      
     type = showerprofile + 2 * pcdetector + 4 * cpvdetector ;
     rp->SetType(type) ; 
@@ -106,3 +154,13 @@ void  AliPHOSPIDv1::SetShowerProfileCuts(Float_t l1m, Float_t l1M, Float_t l2m,
   fLambda2m = l2m ; 
   fLambda2M = l2M ; 
 }
+
+//____________________________________________________________________________
+void  AliPHOSPIDv1::SetRelativeDistanceCut(Float_t CutOnRelativeDistance)
+{
+  // Modifies the parameters used for the particle type identification
+
+  fCutOnRelativeDistance = CutOnRelativeDistance ; 
+}
+
index c429f310c5d411c621f9ba70910548737fbbc7eb..ca646436740902e17b485b44601c0291f520d523 100644 (file)
@@ -24,14 +24,23 @@ class  AliPHOSPIDv1 : public AliPHOSPID {
 
 public:
 
-  AliPHOSPIDv1(): fCutOnDispersion(1.5){}                     
+  AliPHOSPIDv1() 
+  { 
+    fCutOnDispersion = 1.5; 
+    fCutOnRelativeDistance = 3.0 ;
+  }
+                     
   virtual ~ AliPHOSPIDv1(){} ; // dtor
 
+
+  Float_t GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcclu, AliPHOSPpsdRecPoint * PpsdClu, Bool_t &toofar, Option_t * Axis) ; // Relative Distance PPSD-EMC
   virtual void MakeParticles(AliPHOSTrackSegment::TrackSegmentsList * trsl, 
                             AliPHOSRecParticle::RecParticlesList * rpl ) ; // does the job
   virtual void Print(const char *) ; 
   virtual void SetDispersionCutOff(Float_t Dcut) {fCutOnDispersion = Dcut ; }    
   virtual void SetShowerProfileCuts(Float_t l1m, Float_t l1M, Float_t l2m, Float_t l2M) ; 
+  virtual void SetRelativeDistanceCut(Float_t CutOnRelativeDistance) ;
 
  private:
 
@@ -41,6 +50,7 @@ public:
   Float_t fLambda2m ;        // minimum value for second elips axis
   Float_t fLambda2M ;        // maximum value for second elips axis
   Float_t fCutOnDispersion ; // cut on the shower dispersion to distinguish hadronic from EM showers
+  Float_t fCutOnRelativeDistance; //Cut on the relative distance between PPSD and EMC
 
   ClassDef( AliPHOSPIDv1,1)  // Particle identifier implementation version 1
 
index c6d090762d88edb78d06f005c79a3f9443cc020f..22838832d442165755c0adfb809dab2bd6828147 100644 (file)
@@ -172,10 +172,6 @@ Int_t AliPHOSRecPoint::GetPHOSMod()
 
   phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
   fPHOSMod = relid[0];
-  if (fPHOSMod<0 || fPHOSMod>phosgeom->GetNModules() ) {
-    cout << "Wrong PHOS module number is found: " << fPHOSMod << endl;
-    return 0;
-  }
   return fPHOSMod ;
 }
 
@@ -186,7 +182,7 @@ Int_t * AliPHOSRecPoint::GetPrimaries(Int_t & number)
   
   AliPHOSDigit * digit ;
   Int_t index ;
-  Int_t maxcounter = 10 ;
+  Int_t maxcounter = 20 ;
   Int_t counter    = 0 ;
   Int_t * tempo    = new Int_t[maxcounter] ;
   AliPHOSIndexToObject * please = AliPHOSIndexToObject::GetInstance() ;
index 1bf4acdb65c94c7fa39540ddfbeb07b227c1ae43..7c0e0a08e12c91c84c23cfed161f4fe7c57dd7f8 100644 (file)
@@ -55,7 +55,7 @@ ClassImp( AliPHOSTrackSegmentMakerv1)
   fR0 = 10. ;   
   //clusters are sorted in "rows" and "columns" of width geom->GetCrystalSize(0),
   fDelta = fR0 + fGeom->GetCrystalSize(0) ;
-  fMinuit = new TMinuit(100) ;
+  if(!gMinuit) gMinuit = new TMinuit(100) ;
   fUnfoldFlag = kTRUE ; 
 }
 
@@ -63,8 +63,8 @@ ClassImp( AliPHOSTrackSegmentMakerv1)
  AliPHOSTrackSegmentMakerv1::~AliPHOSTrackSegmentMakerv1()
 { 
   // dtor
-   delete fMinuit ; 
+
+  //   delete gMinuit ; 
 }
 
 //____________________________________________________________________________
@@ -123,7 +123,7 @@ Bool_t  AliPHOSTrackSegmentMakerv1::FindFit(AliPHOSEmcRecPoint * emcRP, int * ma
   Double_t p1 = 1.0 ;
   Double_t p2 = 0.0 ;
 
-  gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TgMinuit to reduce function calls  
+  gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls  
   gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient  
   gMinuit->SetMaxIterations(5);
   gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
@@ -155,7 +155,7 @@ void  AliPHOSTrackSegmentMakerv1::FillOneModule(AliPHOSRecPoint::RecPointsList *
                                                Int_t & emcStopedAt, 
                                                Int_t & ppsdStopedAt)
 {
-  // Fill xxxOut arrays with clusters from one PHOS module EMC+PPSD
+  // Fill xxxOut arrays with clusters from one PHOS module
  
   AliPHOSEmcRecPoint *  emcRecPoint  ; 
   AliPHOSPpsdRecPoint * ppsdRecPoint ;
@@ -195,48 +195,6 @@ void  AliPHOSTrackSegmentMakerv1::FillOneModule(AliPHOSRecPoint::RecPointsList *
   ppsdOutUp->Set(inPpsdUp);
   ppsdStopedAt = index ;
    
-}
-//____________________________________________________________________________
-void  AliPHOSTrackSegmentMakerv1::FillOneModule(AliPHOSRecPoint::RecPointsList * emcIn, 
-                                               TArrayI * emcOut, 
-                                               AliPHOSRecPoint::RecPointsList * cpvIn, 
-                                               TArrayI * cpvOut,
-                                               Int_t & phosmod, 
-                                               Int_t & emcStopedAt, 
-                                               Int_t & cpvStopedAt)
-{
-  // Fill xxxOut arrays with clusters from one PHOS module EMC+CPV
-  AliPHOSEmcRecPoint * emcRecPoint  ; 
-  AliPHOSEmcRecPoint * cpvRecPoint ;
-  Int_t index ;
-  
-  Int_t nEmcUnfolded = emcIn->GetEntries() ;
-  emcOut->Set(nEmcUnfolded);
-  Int_t inEmcOut = 0 ;
-  for(index = emcStopedAt; index < nEmcUnfolded; index++){
-    emcRecPoint = (AliPHOSEmcRecPoint *) emcIn->At(index) ;
-    if(emcRecPoint->GetPHOSMod() != phosmod )  
-      break ;
-    
-    emcOut->AddAt(emcRecPoint->GetIndexInList(),inEmcOut) ;
-    inEmcOut++ ; 
-  }
-  emcOut->Set(inEmcOut) ;
-  emcStopedAt = index ;
-
-  cpvOut->Set(cpvIn->GetEntries()) ;
-  Int_t inCpvOut = 0;
-  for(index = cpvStopedAt; index < cpvIn->GetEntries(); index++){
-    cpvRecPoint = (AliPHOSEmcRecPoint *) cpvIn->At(index) ;
-    if(cpvRecPoint->GetPHOSMod() != phosmod )   
-      break ;
-    else  
-      cpvOut->AddAt(index,inCpvOut++) ;
-  }
-  cpvOut->Set(inCpvOut);
-  cpvStopedAt = index ;
-   
 }
 //____________________________________________________________________________
 Float_t  AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcclu,AliPHOSPpsdRecPoint * PpsdClu, Bool_t &toofar)
@@ -432,7 +390,7 @@ void  AliPHOSTrackSegmentMakerv1::MakeTrackSegments(DigitsList * dl,
 
 
   if(fUnfoldFlag){
-    UnfoldAll(dl, emcl) ;      // Unfolds all EMC clusters
+    UnfoldAll(dl, emcl) ; // Unfolds all EMC clusters
   }
 
 
@@ -475,8 +433,8 @@ void  AliPHOSTrackSegmentMakerv1::MakeTrackSegments(DigitsList * dl,
 
 //____________________________________________________________________________
 void  AliPHOSTrackSegmentMakerv1::MakeTrackSegmentsCPV(DigitsList * dl, 
-                                                      AliPHOSRecPoint::RecPointsList * emcl, 
-                                                      AliPHOSRecPoint::RecPointsList * cpvl)
+                                                       AliPHOSRecPoint::RecPointsList * emcl, 
+                                                       AliPHOSRecPoint::RecPointsList * cpvl)
 {
   // Unfold clusters in EMC and CPV and refill reconstructed point lists emcl and ppsdl
   // Yuri Kharlov. 19 October 2000
@@ -528,14 +486,14 @@ void  AliPHOSTrackSegmentMakerv1::UnfoldAll(DigitsList * dl, AliPHOSRecPoint::Re
   
   for(index = 0 ; index < nEmcUnfolded; index++){
 
-    emcRecPoint   = (AliPHOSEmcRecPoint *) emcIn->At(index) ;
+    emcRecPoint = (AliPHOSEmcRecPoint *) emcIn->At(index) ;
     
-    Int_t nMultipl        = emcRecPoint->GetMultiplicity() ; 
-    Int_t   * maxAt       = new Int_t[nMultipl] ;
+    Int_t nMultipl = emcRecPoint->GetMultiplicity() ; 
+    Int_t * maxAt = new Int_t[nMultipl] ;
     Float_t * maxAtEnergy = new Float_t[nMultipl] ;
     Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy) ;
     
-    if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0 
+    if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0       
       UnfoldClusters(dl, emcIn, emcRecPoint, nMax, maxAt, maxAtEnergy) ;
       emcIn->Remove(emcRecPoint); 
       emcIn->Compress() ;
@@ -551,7 +509,9 @@ void  AliPHOSTrackSegmentMakerv1::UnfoldAll(DigitsList * dl, AliPHOSRecPoint::Re
 
   // to set index to new and correct index of old RecPoints
   for( index = 0 ; index < emcIn->GetEntries() ; index++){
+    
     ((AliPHOSEmcRecPoint *) emcIn->At(index))->SetIndexInList(index) ;   
+    
   }
 
   emcIn->Sort() ;
@@ -579,7 +539,6 @@ void  AliPHOSTrackSegmentMakerv1::UnfoldClusters(DigitsList * dl,
     return ;
   }
 
-
   Float_t xDigit ;
   Float_t zDigit ;
   Int_t relid[4] ;
index fe37af8aa9a6f81a64ea919449192bef84090b1b..ce203ded756901623a7fc541ac376ad64a407e28 100644 (file)
@@ -47,14 +47,7 @@ public:
                        TArrayI * ppsdOutLow, 
                        Int_t &PHOSModule, 
                        Int_t & emcStopedAt, 
-                       Int_t & ppsdStopedAt) ; // Fills temporary arrais with clusters from one module EMC+PPSD
-  void    FillOneModule(AliPHOSRecPoint::RecPointsList * emcIn, 
-                       TArrayI * emcOut, 
-                       AliPHOSRecPoint::RecPointsList * cpvIn, 
-                       TArrayI * cpvOut, 
-                       Int_t & PHOSModule, 
-                       Int_t & emcStopedAt, 
-                       Int_t & cpvStopedAt) ;  // Fills temporary arrais with clusters from one module EMC+CPV
+                       Int_t & ppsdStopedAt) ; // Fills temporary arrais with clusters from one module  
   Float_t GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * EmcClu , AliPHOSPpsdRecPoint * Ppsd , Bool_t & TooFar ) ; // see R0
 
   void    MakeLinks(TArrayI * EmcRecPoints, TArrayI * PpsdRecPointsUp, TArrayI * PpsdRecPointsLow, 
@@ -69,14 +62,14 @@ public:
                            AliPHOSRecPoint::RecPointsList * emcl, 
                            AliPHOSRecPoint::RecPointsList * ppsdl, 
                            AliPHOSTrackSegment::TrackSegmentsList * trsl ) ; // does the job
-  void    MakeTrackSegmentsCPV(DigitsList * DL, 
-                           AliPHOSRecPoint::RecPointsList * emcl, 
-                           AliPHOSRecPoint::RecPointsList * ppsdl ) ; // just unfold EMC and CPV clusters
+  virtual void MakeTrackSegmentsCPV(DigitsList * DL, 
+                                AliPHOSRecPoint::RecPointsList * emcl, 
+                                AliPHOSRecPoint::RecPointsList * ppsdl ); // just unfold EMC and CPV clusters
   virtual void SetMaxEmcPpsdDistance(Float_t r){ fR0 = r ;}
   virtual void    SetUnfoldFlag() { fUnfoldFlag = kTRUE ; } ; 
   static Double_t ShowerShape(Double_t r) ; // Shape of shower used in unfolding; class member function (not object member function)
   void    UnfoldAll(DigitsList * Dl, AliPHOSRecPoint::RecPointsList * emcIn) ; 
-                                            // Unfolds and sorts all EMC clusters
+                                             // Unfolds and sorts all EMC clusters
   void  UnfoldClusters(DigitsList * DL, 
                       AliPHOSRecPoint::RecPointsList * emcIn, 
                       AliPHOSEmcRecPoint * iniEmc, 
@@ -95,7 +88,6 @@ public:
 private:
 
   Float_t fDelta ;     // parameter used for sorting
-  TMinuit * fMinuit ;  // Minuit object needed by cluster unfolding
   Float_t fR0 ;        // Maximum distance between a EMC RecPoint and a PPSD RecPoint   
   Bool_t fUnfoldFlag ; // Directive to unfold or not the clusters in case of multiple maxima
 
index 2506971bcf7d72390cc935a1fcec05fb7fdaf5be..0c0116540c8e702c6dc4819b4ff5448c8faca1fd 100644 (file)
@@ -15,6 +15,9 @@
 
 /*
   $Log$
+  Revision 1.1  2000/10/25 08:24:33  schutz
+  New CPV(GPS2) geometry class
+
 */
 
 //_________________________________________________________________________
@@ -53,8 +56,8 @@ AliPPSDGeometry::AliPPSDGeometry()
   fMicromegasWallThickness  = 0.6 ; 
   fNumberOfModulesPhi       = 4 ; 
   fNumberOfModulesZ         = 4 ; 
-  fNumberOfPadsPhi          = 24 ; 
-  fNumberOfPadsZ            = 24 ;   
+  fNumberOfPadsPhi          = 32 ; 
+  fNumberOfPadsZ            = 32 ;   
   fPCThickness              = 0.1 ; 
   fPhiDisplacement          = 0.8 ;  
   fZDisplacement            = 0.8 ;