// --- Standard library ---
// --- AliRoot header files ---
+#include "AliGenerator.h"
#include "AliEMCALPIDv1.h"
-#include "AliEMCALTrackSegment.h"
#include "AliEMCALRecParticle.h"
#include "AliEMCALGetter.h"
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()
{
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();
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) ;
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);
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++ ;
}
{
// 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)
{
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() ) ;
}
//____________________________________________________________________________