]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliEMCALPIDv1.cxx
Station2 detailed geometry builder (Sanjoy, Sukalyan, Shakeel)
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALPIDv1.cxx
index 0eef8756ae94bc1235d5967914702674b516a4de..391696736bbc4137c4248740c0ea5c57ca86127c 100644 (file)
@@ -29,8 +29,8 @@
   // --- Standard library ---
   
   // --- AliRoot header files ---
+#include "AliGenerator.h"
 #include "AliEMCALPIDv1.h"
-#include "AliEMCALTrackSegment.h"
 #include "AliEMCALRecParticle.h"
 #include "AliEMCALGetter.h"
   
@@ -79,6 +79,47 @@ const TString AliEMCALPIDv1::BranchName() const
   return GetName() ;
 }
  
+//____________________________________________________________________________
+Float_t  AliEMCALPIDv1::GetCalibratedEnergy(Float_t e) const
+{
+//      It calibrates Energy depending on the recpoint energy.
+//      The energy of the reconstructed cluster is corrected with 
+//      the formula A + B* E  + C* E^2, whose parameters where obtained 
+//      through the study of the reconstructed energy distribution of 
+//      monoenergetic photons.
+  //Float_t p[]={0.,0.,0.};
+  //for (Int_t i=0; i<3; i++) p[i] = GetParameterCalibration(i);
+  Float_t enerec = e ; // p[0] +  p[1]*e + p[2]*e*e;
+  return enerec ;
+
+}
+//____________________________________________________________________________
+TVector3 AliEMCALPIDv1::GetMomentumDirection(AliEMCALRecPoint * emc)const 
+{ 
+  // Calculates the momentum direction:
+  // direction is given by IP and this RecPoint
+  
+
+  TVector3 dir(0,0,0) ; 
+  TVector3 emcglobalpos ;
+  // TMatrix  dummy ;
+  
+  emc->GetGlobalPosition(emcglobalpos) ;
+  
+
+  dir = emcglobalpos ;  
+  // dir.SetMag(1.) ; Removed to avoid warings !!!!!!!!!!!!!! TO BE REVISED
+
+  //account correction to the position of IP
+  Float_t xo,yo,zo ; //Coordinates of the origin
+  gAlice->Generator()->GetOrigin(xo,yo,zo) ;
+  TVector3 origin(xo,yo,zo);
+  dir = dir - origin ;
+
+  return dir ;  
+}
+
 //____________________________________________________________________________
 void AliEMCALPIDv1::Init()
 {
@@ -115,28 +156,29 @@ void  AliEMCALPIDv1::Exec(Option_t * option)
     return ; 
   }
   AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
+  
+  if (fLastEvent == -1) 
+    fLastEvent = gime->MaxEvent() - 1 ;
+  else 
+    fLastEvent = TMath::Min(fLastEvent,gime->MaxEvent());
+  Int_t nEvents   = fLastEvent - fFirstEvent + 1;
 
-  Int_t nevents = gime->MaxEvent() ;      
   Int_t ievent ;
 
-
-  for(ievent = 0; ievent < nevents; ievent++){
-    gime->Event(ievent,"TR") ;
-    if(gime->TrackSegments() && //Skip events, where no track segments made
-       gime->TrackSegments()->GetEntriesFast()) {   
-      MakeRecParticles() ;
-      WriteRecParticles();
-      if(strstr(option,"deb"))
-       PrintRecParticles(option) ;
-      //increment the total number of rec particles per run 
-      fRecParticlesInRun += gime->RecParticles()->GetEntriesFast() ; 
-    }
+  for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
+    gime->Event(ievent,"R") ;
+    MakeRecParticles() ;
+    WriteRecParticles();
+    if(strstr(option,"deb"))
+      PrintRecParticles(option) ;
+    //increment the total number of rec particles per run 
+    fRecParticlesInRun += gime->RecParticles()->GetEntriesFast() ;  
   }
   if(strstr(option,"tim")){
     gBenchmark->Stop("EMCALPID");
-    Info("Exec", "took %f seconds for PID %f seconds per event", 
+    printf("Exec: took %f seconds for PID %f seconds per event", 
         gBenchmark->GetCpuTime("EMCALPID"),  
-        gBenchmark->GetCpuTime("EMCALPID")/nevents) ;
+        gBenchmark->GetCpuTime("EMCALPID")/nEvents) ;
   } 
 
   Unload();
@@ -149,48 +191,27 @@ void  AliEMCALPIDv1::MakeRecParticles(){
   
   AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
   TObjArray * aECARecPoints = gime->ECARecPoints() ; 
-  TObjArray * aPRERecPoints = gime->PRERecPoints() ; 
-  TObjArray * aHCARecPoints = gime->HCARecPoints() ; 
-  TClonesArray * trackSegments = gime->TrackSegments() ; 
-  if ( !aECARecPoints || !aPRERecPoints || !aHCARecPoints || !trackSegments ) {
+  if ( !aECARecPoints ) {
     Fatal("MakeRecParticles", "RecPoints or TrackSegments not found !") ;  
   }
   TClonesArray * recParticles  = gime->RecParticles() ; 
   recParticles->Clear();
 
-  TIter next(trackSegments) ; 
-  AliEMCALTrackSegment * ts ; 
+  TIter next(aECARecPoints) ; 
+  AliEMCALRecPoint * eca ; 
   Int_t index = 0 ; 
   AliEMCALRecParticle * rp ; 
-  while ( (ts = (AliEMCALTrackSegment *)next()) ) {
+  while ( (eca = (AliEMCALRecPoint *)next()) ) {
     
     new( (*recParticles)[index] ) AliEMCALRecParticle() ;
     rp = (AliEMCALRecParticle *)recParticles->At(index) ; 
-    rp->SetTrackSegment(index) ;
+    rp->SetRecPoint(index) ;
     rp->SetIndexInList(index) ;
        
-    AliEMCALTowerRecPoint * eca = 0 ;
-    if(ts->GetECAIndex()>=0)
-      eca = dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(ts->GetECAIndex())) ;
-    
-    AliEMCALTowerRecPoint * pre = 0 ;
-    if(ts->GetPREIndex()>=0)
-      pre = dynamic_cast<AliEMCALTowerRecPoint *>(aPRERecPoints->At(ts->GetPREIndex())) ;
-    
-    AliEMCALTowerRecPoint * hca = 0 ;
-    if(ts->GetHCAIndex()>=0)
-      hca = dynamic_cast<AliEMCALTowerRecPoint *>(aHCARecPoints->At(ts->GetHCAIndex())) ;
-
     // Now set type (reconstructed) of the particle
 
     // Choose the cluster energy range
     
-    if (!eca) {
-      Fatal("MakeRecParticles", "-> emcal(%d) = %d", ts->GetECAIndex(), eca ) ;
-    }
-
-    Float_t    e = eca->GetEnergy() ;   
-    
     Float_t  lambda[2] ;
     eca->GetElipsAxis(lambda) ;
     
@@ -207,13 +228,13 @@ void  AliEMCALPIDv1::MakeRecParticles(){
       emaxdtotal=eca->GetMaximalEnergy()/eca->GetEnergy(); 
     }
     
-    //    Float_t time = ecal->GetTime() ;
+    //    Float_t time = eca->GetTime() ;
       
     //Set momentum, energy and other parameters 
-    Float_t  enca = e;
-    TVector3 dir(0., 0., 0.) ; 
-    dir.SetMag(enca) ;
-    rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),enca) ;
+    Float_t  encal = GetCalibratedEnergy(eca->GetEnergy());
+    TVector3 dir   = GetMomentumDirection(eca) ; 
+    // dir.SetMag(encal) ;Removed to avoid warings !!!!!!!!!!!!!! TO BE REVISED
+    rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),encal) ;
     rp->SetCalcMass(0);
     rp->Name(); //If photon sets the particle pdg name to gamma
     rp->SetProductionVertex(0,0,0,0);
@@ -222,6 +243,14 @@ void  AliEMCALPIDv1::MakeRecParticles(){
     rp->SetFirstDaughter(-1);
     rp->SetLastDaughter(-1);
     rp->SetPolarisation(0,0,0);
+    //Set the position in global coordinate system from the RecPoint
+    //AliEMCALGeometry * geom = gime->EMCALGeometry() ; 
+    //AliEMCALTowerRecPoint  * erp = gime->ECARecPoint(rp->GetEMCALRPIndex()) ; 
+    TVector3 pos ; 
+    //geom->GetGlobal(erp, pos) ; !!!!!!!!!! to check 
+    rp->SetPos(pos);
+
+
     index++ ; 
   }
   
@@ -232,12 +261,20 @@ void  AliEMCALPIDv1:: Print(Option_t * /*option*/) const
 {
   // Print the parameters used for the particle type identification
 
-    Info("Print", "=============== AliEMCALPID1 ================") ;
+    printf("Print: =============== AliEMCALPID1 ================") ;
     printf("Making PID\n");
     printf("    Pricipal analysis file from 0.5 to 100 %s\n", fFileName.Data() ) ; 
     printf("    Name of parameters file     %s\n", fFileNamePar.Data() )  ;
 }
 
+//____________________________________________________________________________
+void  AliEMCALPIDv1::Print() const
+{
+  // Print the parameters used for the particle type identification
+
+    Info("Print", "=============== AliEMCALPIDv1 ================") ;
+}
+
 //____________________________________________________________________________
 void AliEMCALPIDv1::PrintRecParticles(Option_t * option)
 {
@@ -247,28 +284,22 @@ void AliEMCALPIDv1::PrintRecParticles(Option_t * option)
 
   TClonesArray * recParticles = gime->RecParticles() ; 
 
-  TString message ; 
-  message  = "\nevent " ;
-  message += gAlice->GetEvNumber() ; 
-  message += "       found " ; 
-  message += recParticles->GetEntriesFast(); 
-  message += " RecParticles\n" ; 
+  printf("\nevent %i", gAlice->GetEvNumber()); 
+  printf("       found %i", recParticles->GetEntriesFast()); 
+  printf(" RecParticles\n"); 
 
   if(strstr(option,"all")) {  // printing found TS
-    message += "\n  PARTICLE         Index    \n" 
+    printf("\n  PARTICLE         Index    \n")
     
     Int_t index ;
     for (index = 0 ; index < recParticles->GetEntries() ; index++) {
       AliEMCALRecParticle * rp = (AliEMCALRecParticle * ) recParticles->At(index) ;       
-      message += "\n" ;
-      message += rp->Name().Data() ;  
-      message += " " ;
-      message += rp->GetIndexInList() ;  
-      message += " " ;
-      message += rp->GetType()  ;
+      printf("\n");
+      printf(rp->Name().Data());  
+      printf(" %i", rp->GetIndexInList());  
+      printf(" %i", rp->GetType());
     }
   }
-  Info("Print", message.Data() ) ; 
 }
 
 //____________________________________________________________________________