updating

diff --git a/PHOS/AnaESD.C b/PHOS/AnaESD.C
index e5444710d4747fb85c4e7124f686c2ed4cbee5a1..165561e2551a33468075b58fc1cb4b53e7302d18 100644 (file)
--- a/PHOS/AnaESD.C
+++ b/PHOS/AnaESD.C
@@ -25,6 +25,7 @@
 #if !defined(__CINT__) || defined(__MAKECINT__)
 #include "TFile.h"
 #include "TMath.h"
+#include "TH1D.h"
 #include "AliPHOSGetter.h"
 #include "AliPHOSGeometry.h"
 #include "Riostream.h"
@@ -36,79 +37,98 @@
 #include "AliKalmanTrack.h"
 #endif
 
+Double_t Match(TParticle * pp, AliESDtrack * cp) ; 
+TH1D * heta = new TH1D("heta", "Eta correlation", 100, -2., 2.) ; 
+TH1D * hphi = new TH1D("hphi", "Phi correlation", 360, 0., 360.) ; 
+ 
 void Ana() 
-{
+{ 
   AliPHOSGetter * gime = AliPHOSGetter::Instance("galice.root") ; 
   Int_t nEvent = gime->MaxEvent() ;  
   Int_t event ; 
   AliESD * esd = 0 ;
   for (event = 0 ; event < nEvent; event++) {
     esd = gime->ESD(event) ; 
-    esd->Print();  
-    Int_t index ;
+    //esd->Print();  
+    Int_t caloindex ;
     // Calorimeter tracks 
     AliESDCaloTrack * ct ; 
-    AliPHOSRecParticle * pp ; 
-    AliEMCALRecParticle * ep ; 
-    for (index = 0 ; index < esd->GetNumberOfCaloTracks() ; index++) {
-      ct = esd->GetCaloTrack(index) ;
-      pp = dynamic_cast<AliPHOSRecParticle*>(ct->GetRecParticle()) ; 
-      ep = dynamic_cast<AliEMCALRecParticle*>(ct->GetRecParticle()) ; 
-      if (pp) { 
-       TVector3 pos = pp->GetPos() ;
-       cout << "PHOS particle # " << index << " pos (" 
-            << pos.X() << ", " << pos.Y() << ", " <<pos.Z() << ") : (" << pos.Eta() 
-            << ", " << pos.Phi()*TMath::RadToDeg() << ") : " 
-            << TMath::Sqrt(pos.X()*pos.X() + pos.Y()*pos.Y() ) << endl ; 
-       Int_t n ; 
-       Double_t z,x ; 
-       gime->PHOSGeometry()->ImpactOnEmc(pos.Theta(), pos.Phi(), n, z, x) ; 
-       if (n) 
-         cout << "Matching: " << n << " " << z << " " << x << endl ; 
-      }
-      if(ep) { 
-       TVector3 pos = ep->GetPos() ;
-       //cout << "EMCAL particle # " << index << " pos " << 
-       //  TMath::Sqrt(pos.X()*pos.X() + pos.Y()*pos.Y() ) << endl ; 
+    for (caloindex = 0 ; caloindex < esd->GetNumberOfCaloTracks() ; caloindex++) {
+      // get the calorimeter type of particles (PHOS or EMCAL)
+      ct = esd->GetCaloTrack(caloindex) ;
+      TParticle * part = ct->GetRecParticle() ; 
+
+      AliESDtrack * cp ; 
+      Int_t cpindex ; 
+      for (cpindex = 0 ; cpindex < esd->GetNumberOfTracks() ; cpindex++) {
+       // get the charged tracks from central tracking
+       cp = esd->GetTrack(cpindex) ;
+       Double_t dist = Match(part, cp) ; 
+       
+       if (dist < 99999.) 
+         cout << "================ Distance = " << dist << endl ; 
       }
     }
-    //Charged tracks from central tracking
-    AliESDtrack * cp ; 
-    for (index = 0 ; index < esd->GetNumberOfTracks() ; index++) {
-      cp = esd->GetTrack(index) ;
-      Double_t xyzAtPHOS[3] ; 
-      cp->GetOuterXYZ(xyzAtPHOS) ; 
-      // check if the track reaches PHOS
-      if ( (xyzAtPHOS[0] +  xyzAtPHOS[1] + xyzAtPHOS[2]) == 0.)
-       continue;
-      // the next check are only if we want high quality tracks 
+  }
+  //  heta->Draw() ; 
+  hphi->Draw() ; 
+}
+Double_t Match(TParticle * part, AliESDtrack * cp) 
+{
+  // Calculates the distance (x,z) between  the particle detected by PHOS and 
+  // the charged particle reconstructed by the global tracking 
+
+  Double_t dist = 99999. ;
+   
+  AliPHOSRecParticle * pp  = dynamic_cast<AliPHOSRecParticle*>(part) ;  
+  AliEMCALRecParticle * ep = dynamic_cast<AliEMCALRecParticle*>(part) ;  
+   
+  Int_t phN ; 
+  Double_t phZ, phX ; 
+  
+  if (pp) { // it is a PHOS particle 
+    
+ //    cout << "PHOS particle # " << " pos (" 
+//      << pp->GetPos().X() << ", " << pp->GetPos().Y() << ", " << pp->GetPos().Z() << ")" << endl ;
+    
+    AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
+    gime->PHOSGeometry()->ImpactOnEmc(*pp, phN, phZ, phX) ; 
+    Double_t xyzAtPHOS[3] ; 
+    cp->GetOuterXYZ(xyzAtPHOS) ; 
+    if ( (xyzAtPHOS[0] +  xyzAtPHOS[1] + xyzAtPHOS[2]) != 0.) { //it has reached PHOS
+      //the next check are only if we want high quality tracks 
       //       ULong_t status = cp->GetStatus() ;  
       //       if ((status & AliESDtrack::kTRDput)==0) 
-      //       continue;
+      //       do not continue;
       //       if ((status & AliESDtrack::kTRDStop)!=0) 
-      //       continue;
-
-      // Gets the Global coordinate of the track at the entrance of PHOS 
-     
-      TVector3 poscp(xyzAtPHOS[0], xyzAtPHOS[1], xyzAtPHOS[2]) ; 
-      cout << "Charged particle # " << index << " pos (" 
-          << poscp.X() << ", " << poscp.Y() << ", " <<poscp.Z() << ") : (" << poscp.Eta() 
-          << ", " << poscp.Phi()*TMath::RadToDeg() << ") : " 
-          << TMath::Sqrt(poscp.X()*poscp.X() + poscp.Y()*poscp.Y() ) << endl ;
-      Int_t n ;
-      Double_t z,x ; 
-      gime->PHOSGeometry()->ImpactOnEmc(poscp.Theta(), poscp.Phi(), n, z, x) ; 
-      if (n) 
-       cout << "Matching: " << n << " " << z << " " << x << endl ; 
-      // Does the matching with PHOS/EMCAL
-//       for (index = 0 ; index < esd->GetNumberOfCaloTracks() ; index++) {
-//     ct = esd->GetCaloTrack(index) ;
-//     pp = dynamic_cast<AliPHOSRecParticle*>(ct->GetRecParticle()) ; 
-//     ep = dynamic_cast<AliEMCALRecParticle*>(ct->GetRecParticle()) ; 
-//     if (pp) {
-//       TVector3 pos = pp->GetPos() ; 
-//     }
-//       }
+      //       do not continue;  
+//       cout << "Charged particle # " << " pos (" 
+//        << xyzAtPHOS[0] << ", " << xyzAtPHOS[1] << ", " << xyzAtPHOS[2] << ")" <<  endl ;     
+      TVector3 poscp(xyzAtPHOS[0], xyzAtPHOS[1], xyzAtPHOS[2]) ;
+      Int_t cpN ;
+      Double_t cpZ,cpX ; 
+      gime->PHOSGeometry()->ImpactOnEmc(poscp, cpN, cpZ, cpX) ; 
+      if (cpN) {// we are inside the PHOS acceptance 
+//     cout << "Charged Matching 1: " << cpN << " " << cpZ << " " << cpX << endl ; 
+//     cout << "Charged Matching 2: " << phN << " " << phZ << " " << phX << endl ; 
+       dist = TMath::Sqrt( (cpZ-phZ)*(cpZ-phZ) + (cpX-phX)*(cpX-phX)) ;  
+      } 
+      Double_t phTheta = pp->Theta() ; 
+      Double_t phPhi   = pp->Phi() ;
+      TParticle tempo ; 
+      tempo.SetMomentum(xyzAtPHOS[0], xyzAtPHOS[1], xyzAtPHOS[2], 0.) ;  
+      Double_t cpTheta = tempo.Theta() ; 
+      Double_t cpPhi   = tempo.Phi() ;
+      //cout << phTheta << " " << phPhi << " " << endl 
+      //cout <<         cpTheta << " " << cpPhi-phPhi << " " << endl ; 
+      heta->Fill( (phTheta - cpTheta)*TMath::RadToDeg() ) ; 
+      hphi->Fill( (phPhi - cpPhi)*TMath::RadToDeg() ) ; 
     }
   }
+  
+  if (ep) {
+    //cout << "EMCAL particle # " << endl ; 
+  }
+  return dist ; 
 }
+