* provided "as is" without express or implied warranty. *
**************************************************************************/
+/* $Id$ */
+
//_________________________________________________________________________
-// Algorythm class to analyze PHOS events
-//*-- Y. Schutz : SUBATECH
+// Algorythm class to analyze PHOSv1 events:
+// Construct histograms and displays them.
+// Use the macro EditorBar.C for best access to the functionnalities
+//
+//*-- Author: Y. Schutz (SUBATECH) & Gines Martinez (SUBATECH)
//////////////////////////////////////////////////////////////////////////////
// --- ROOT system ---
// --- Standard library ---
-#include <iostream>
-#include <cstdio>
+#include <iostream.h>
+#include <stdio.h>
// --- AliRoot header files ---
#include "AliPHOSDigit.h"
#include "AliPHOSTrackSegment.h"
#include "AliPHOSRecParticle.h"
+#include "AliPHOSIndexToObject.h"
ClassImp(AliPHOSAnalyze)
//____________________________________________________________________________
AliPHOSAnalyze::AliPHOSAnalyze()
{
- // ctor
+ // default ctor (useless)
fRootFile = 0 ;
}
//____________________________________________________________________________
AliPHOSAnalyze::AliPHOSAnalyze(Text_t * name)
{
- // ctor
+ // ctor: analyze events from root file "name"
Bool_t ok = OpenRootFile(name) ;
if ( !ok ) {
cout << " AliPHOSAnalyze > Error opening " << name << endl ;
}
else {
- gAlice = (AliRun*) fRootFile->Get("gAlice");
- fPHOS = (AliPHOSv0 *)gAlice->GetDetector("PHOS") ;
- fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() ) ;
- fEvt = -999 ;
+ //========== Get AliRun object from file
+ gAlice = (AliRun*) fRootFile->Get("gAlice") ;
+
+ //=========== Get the PHOS object and associated geometry from the file
+ fPHOS = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
+ fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
+
+ //========== Initializes the Index to Object converter
+ fObjGetter = AliPHOSIndexToObject::GetInstance(fPHOS) ;
+ //========== Current event number
+ fEvt = -999 ;
+
}
+ fClu = 0 ;
+ fPID = 0 ;
+ fTrs = 0 ;
+ fRec = 0 ;
+ ResetHistograms() ;
+}
+
+//____________________________________________________________________________
+AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana)
+{
+ // copy ctor
+ ( (AliPHOSAnalyze &)ana ).Copy(*this) ;
+}
+
+//____________________________________________________________________________
+void AliPHOSAnalyze::Copy(TObject & obj)
+{
+ // copy an analysis into an other one
+ TObject::Copy(obj) ;
+ // I do nothing more because the copy is silly but the Code checkers requires one
}
//____________________________________________________________________________
{
// dtor
- fRootFile->Close() ;
- delete fRootFile ;
- fRootFile = 0 ;
+ if (fRootFile->IsOpen() )
+ fRootFile->Close() ;
+ if(fRootFile)
+ delete fRootFile ;
- delete fPHOS ;
- fPHOS = 0 ;
+ if(fPHOS)
+ delete fPHOS ;
- delete fClu ;
- fClu = 0 ;
+ if(fClu)
+ delete fClu ;
- delete fPID ;
- fPID = 0 ;
+ if(fPID)
+ delete fPID ;
- delete fRec ;
- fRec = 0 ;
+ if(fRec)
+ delete fRec ;
- delete fTrs ;
- fTrs = 0 ;
+ if(fTrs)
+ delete fTrs ;
}
//____________________________________________________________________________
-void AliPHOSAnalyze::AnalyzeOneEvent(Int_t evt)
-{
- Bool_t ok = Init(evt) ;
+void AliPHOSAnalyze::ActivePPSD(Int_t Nevents=1){
- if ( ok ) {
- //=========== Get the number of entries in the Digits array
-
- Int_t nId = fPHOS->Digits()->GetEntries();
- printf("AnalyzeOneEvent > Number of entries in the Digit array is %d \n",nId);
-
- //=========== Do the reconstruction
-
- cout << "AnalyzeOneEvent > Found " << nId << " digits in PHOS" << endl ;
-
- fPHOS->Reconstruction(fRec);
-
- // =========== End of reconstruction
-
- cout << "AnalyzeOneEvent > event # " << fEvt << " processed" << endl ;
- } // ok
- else
- cout << "AnalyzeOneEvent > filed to process event # " << evt << endl ;
+ fhEnergyCorrelations = new TH2F("hEnergyCorrelations","hEnergyCorrelations",40, 0., 0.15, 30, 0., 3.e-5);
+ //========== Create the Clusterizer
+ fClu = new AliPHOSClusterizerv1() ;
+ fClu->SetEmcEnergyThreshold(0.01) ;
+ fClu->SetEmcClusteringThreshold(0.20) ;
+ fClu->SetPpsdEnergyThreshold (0.0000002) ;
+ fClu->SetPpsdClusteringThreshold(0.0000001) ;
+ fClu->SetLocalMaxCut(0.02) ;
+ fClu->SetCalibrationParameters(0., 0.00000001) ;
+ Int_t ievent;
+
+ for ( ievent=0; ievent<Nevents; ievent++)
+ {
+
+ //========== Event Number>
+ if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
+ cout << "AnalyzeResolutions > " << "Event is " << ievent << endl ;
+
+ //=========== Connects the various Tree's for evt
+ gAlice->GetEvent(ievent);
+
+ //=========== Gets the Kine TTree
+ gAlice->TreeK()->GetEvent(0) ;
+
+ //=========== Get the Digit Tree
+ gAlice->TreeD()->GetEvent(0) ;
+
+ //========== Creating branches ===================================
+ AliPHOSRecPoint::RecPointsList ** EmcRecPoints = fPHOS->EmcRecPoints() ;
+ gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP", EmcRecPoints ) ;
+
+ AliPHOSRecPoint::RecPointsList ** PpsdRecPoints = fPHOS->PpsdRecPoints() ;
+ gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", PpsdRecPoints ) ;
+
+ AliPHOSTrackSegment::TrackSegmentsList ** TrackSegmentsList = fPHOS->TrackSegments() ;
+ if( (*TrackSegmentsList) )
+ (*TrackSegmentsList)->Clear() ;
+ gAlice->TreeR()->SetBranchAddress( "PHOSTS", TrackSegmentsList ) ;
+
+ AliPHOSRecParticle::RecParticlesList ** RecParticleList = fPHOS->RecParticles() ;
+ if( (*RecParticleList) )
+ (*RecParticleList)->Clear() ;
+ gAlice->TreeR()->SetBranchAddress( "PHOSRP", RecParticleList ) ;
+
+
+ //=========== Gets the Reconstraction TTree
+ gAlice->TreeR()->GetEvent(0) ;
+
+ AliPHOSPpsdRecPoint * RecPoint ;
+ Int_t relid[4] ;
+ TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
+ while( ( RecPoint = (AliPHOSPpsdRecPoint *)nextRP() ) )
+ {
+ if(!(RecPoint->GetUp()) ) {
+ AliPHOSDigit *digit ;
+ Int_t iDigit ;
+ for(iDigit = 0; iDigit < fPHOS->Digits()->GetEntries(); iDigit++)
+ {
+ digit = (AliPHOSDigit *) fPHOS->Digits()->At(iDigit) ;
+ fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ if((relid[2]==1)&&(relid[3]==1)&&(relid[0]==RecPoint->GetPHOSMod())){
+ Float_t ConvertorEnergy = fClu->Calibrate(digit->GetAmp()) ;
+ fhEnergyCorrelations->Fill(ConvertorEnergy,RecPoint->GetTotalEnergy() );
+ break ;
+ }
+ }
+ break ;
+ }
+ }
+ }
+ SaveHistograms() ;
+ fhEnergyCorrelations->Draw("BOX") ;
}
+
//____________________________________________________________________________
- void AliPHOSAnalyze::AnalyzeManyEvents(Int_t Nevents, Int_t module) // analyzes many events
+void AliPHOSAnalyze::AnalyzeManyEvents(Int_t Nevents, Int_t module)
{
+ // analyzes Nevents events in a single PHOS module
+ // Events should be reconstructed by Reconstruct()
if ( fRootFile == 0 )
cout << "AnalyzeManyEvents > " << "Root File not openned" << endl ;
else
{
- //========== Get AliRun object from file
- gAlice = (AliRun*) fRootFile->Get("gAlice") ;
- //=========== Get the PHOS object and associated geometry from the file
- fPHOS = (AliPHOSv0 *)gAlice->GetDetector("PHOS") ;
- fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
//========== Booking Histograms
cout << "AnalyzeManyEvents > " << "Booking Histograms" << endl ;
BookingHistograms();
+
Int_t ievent;
- Int_t relid[4] ;
+ Int_t relid[4] ;
AliPHOSDigit * digit ;
AliPHOSEmcRecPoint * emc ;
AliPHOSPpsdRecPoint * ppsd ;
// AliPHOSTrackSegment * tracksegment ;
AliPHOSRecParticle * recparticle;
+
for ( ievent=0; ievent<Nevents; ievent++)
{
- if (ievent==0) cout << "AnalyzeManyEvents > " << "Starting Analyzing " << endl ;
- //========== Create the Clusterizer
- fClu = new AliPHOSClusterizerv1() ;
- fClu->SetEmcEnergyThreshold(0.025) ;
- fClu->SetEmcClusteringThreshold(0.75) ;
- fClu->SetPpsdEnergyThreshold (0.0000002) ;
- fClu->SetPpsdClusteringThreshold(0.0000001) ;
- fClu->SetLocalMaxCut(0.03) ;
- fClu->SetCalibrationParameters(0., 0.00000001) ;
- //========== Creates the track segment maker
- fTrs = new AliPHOSTrackSegmentMakerv1() ;
- //========== Creates the particle identifier
- fPID = new AliPHOSPIDv1() ;
- fPID->SetShowerProfileCuts(0.5, 1.5, 0.5, 1.5 ) ;
- fPID->Print() ;
- //========== Creates the Reconstructioner
- fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
- //========== Event Number
- if ( ( log10(ievent+1) - (Int_t)(log10(ievent+1)) ) == 0. ) cout << "AnalyzeManyEvents > " << "Event is " << ievent << endl ;
+ //========== Event Number>
+ if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
+ cout << "AnalyzeManyEvents > " << "Event is " << ievent << endl ;
+
//=========== Connects the various Tree's for evt
gAlice->GetEvent(ievent);
+
//=========== Gets the Digit TTree
gAlice->TreeD()->GetEvent(0) ;
+
//=========== Gets the number of entries in the Digits array
TIter nextdigit(fPHOS->Digits()) ;
while( ( digit = (AliPHOSDigit *)nextdigit() ) )
if (relid[1]>16) fhConvertorDigit->Fill(fClu->Calibrate(digit->GetAmp()));
}
}
- //=========== Do the reconstruction
- fPHOS->Reconstruction(fRec);
+
+
//=========== Cluster in module
- TIter nextEmc(fPHOS->EmcClusters() ) ;
+ TIter nextEmc(*fPHOS->EmcRecPoints() ) ;
while((emc = (AliPHOSEmcRecPoint *)nextEmc()))
{
if ( emc->GetPHOSMod() == module )
{
fhEmcCluster->Fill( emc->GetTotalEnergy() );
- TIter nextPpsd( fPHOS->PpsdClusters()) ;
+ TIter nextPpsd( *fPHOS->PpsdRecPoints()) ;
while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
{
if ( ppsd->GetPHOSMod() == module )
}
}
}
+
//=========== Cluster in module PPSD Down
- TIter nextPpsd(fPHOS->PpsdClusters() ) ;
+ TIter nextPpsd(*fPHOS->PpsdRecPoints() ) ;
while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
{
if ( ppsd->GetPHOSMod() == module )
if (ppsd->GetUp()) fhVetoCluster ->Fill(ppsd->GetTotalEnergy()) ;
}
}
+
//========== TRackSegments in the event
- TIter nextRecParticle(fPHOS->RecParticles() ) ;
+ TIter nextRecParticle(*fPHOS->RecParticles() ) ;
while((recparticle = (AliPHOSRecParticle *)nextRecParticle()))
{
if ( recparticle->GetPHOSTrackSegment()->GetPHOSMod() == module )
{
- cout << "Particle type is " << recparticle->GetType() << endl ;
- switch(recparticle->GetType())
- {
- case kGAMMA:
- fhPhotonEnergy->Fill(recparticle->Energy() ) ;
- //fhPhotonPositionX->Fill(recpart. ) ;
- //fhPhotonPositionY->Fill(recpart. ) ;
- cout << "PHOTON" << endl;
- break;
- case kELECTRON:
- fhElectronEnergy->Fill(recparticle->Energy() ) ;
- //fhElectronPositionX->Fill(recpart. ) ;
- //fhElectronPositionY->Fill(recpart. ) ;
- cout << "ELECTRON" << endl;
- break;
- case kNEUTRALHADRON:
- fhNeutralHadronEnergy->Fill(recparticle->Energy() ) ;
- //fhNeutralHadronPositionX->Fill(recpart. ) ;
- //fhNeutralHadronPositionY->Fill(recpart. ) ;
- cout << "NEUTRAl HADRON" << endl;
- break ;
- case kNEUTRALEM:
- fhNeutralEMEnergy->Fill(recparticle->Energy() ) ;
- //fhNeutralEMPositionX->Fill(recpart. ) ;
- //fhNeutralEMPositionY->Fill(recpart. ) ;
- //cout << "NEUTRAL EM" << endl;
- break ;
- case kCHARGEDHADRON:
- fhChargedHadronEnergy->Fill(recparticle->Energy() ) ;
- //fhChargedHadronPositionX->Fill(recpart. ) ;
- //fhChargedHadronPositionY->Fill(recpart. ) ;
- cout << "CHARGED HADRON" << endl;
- break ;
-
- }
+ cout << "Particle type is " << recparticle->GetType() << endl ;
+ Int_t numberofprimaries = 0 ;
+ Int_t * listofprimaries = recparticle->GetPrimaries(numberofprimaries) ;
+ cout << "Number of primaries = " << numberofprimaries << endl ;
+ Int_t index ;
+ for ( index = 0 ; index < numberofprimaries ; index++)
+ cout << " primary # " << index << " = " << listofprimaries[index] << endl ;
}
}
- // Deleting fClu, fTrs, fPID et fRec
- fClu->Delete();
- fTrs->Delete();
- fPID->Delete();
- fRec->Delete();
-
} // endfor
- SavingHistograms();
+ SaveHistograms();
} // endif
} // endfunction
+//____________________________________________________________________________
+ void AliPHOSAnalyze::Reconstruct(Int_t Nevents,Int_t FirstEvent )
+{
+ Int_t ievent ;
+ for ( ievent=FirstEvent; ievent<Nevents; ievent++)
+ {
+ if (ievent==FirstEvent)
+ {
+ cout << "Analyze > Starting Reconstructing " << endl ;
+ //========== Create the Clusterizer
+ fClu = new AliPHOSClusterizerv1() ;
+ fClu->SetEmcEnergyThreshold(0.03) ;
+ fClu->SetEmcClusteringThreshold(0.20) ;
+ fClu->SetPpsdEnergyThreshold (0.0000001) ;
+ fClu->SetPpsdClusteringThreshold(0.0000001) ;
+ fClu->SetLocalMaxCut(0.02) ;
+ fClu->SetCalibrationParameters(0., 0.00000001) ;
+
+ //========== Creates the track segment maker
+ fTrs = new AliPHOSTrackSegmentMakerv1() ;
+ // fTrs->UnsetUnfoldFlag() ;
+
+ //========== Creates the particle identifier
+ fPID = new AliPHOSPIDv1() ;
+ fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ;
+ fPID->SetDispersionCutOff(2.0) ;
+ fPID->SetRelativeDistanceCut(3.) ;
+
+ //========== Creates the Reconstructioner
+ fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
+ // fRec -> SetDebugReconstruction(kTRUE);
+
+ }
+
+ //========== Event Number>
+ // if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
+ cout << "Reconstruct > Event is " << ievent << 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);
+ }
+
+ fClu->Delete();
+ fClu=0 ;
+ fTrs->Delete();
+ fTrs = 0 ;
+ fPID->Delete();
+ fPID = 0 ;
+ fRec->Delete();
+ fRec = 0 ;
+
+}
+//____________________________________________________________________________
+ void AliPHOSAnalyze::InvariantMass(Int_t Nevents )
+{
+ // 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 ;
+
+ 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 ;
+ }
+
+
+ //Now calculate invariant mass:
+ Int_t irp1,irp2 ;
+ Int_t NCurEvent = 0 ;
+
+ for(irp1 = 0; irp1 < AllRecParticleList->GetEntries()-1; irp1++){
+ AliPHOSRecParticle * rp1 = (AliPHOSRecParticle *)AllRecParticleList->At(irp1) ;
+
+ 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
+
+
+ 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();
+
+}
+
+//____________________________________________________________________________
+ void AliPHOSAnalyze::AnalyzeResolutions(Int_t Nevents )
+{
+ // analyzes Nevents events and calculate Energy and Position resolution as well as
+ // probaility of correct indentifiing of the incident particle
+
+ //========== Booking Histograms
+ cout << "AnalyzeResolutions > " << "Booking Histograms" << endl ;
+ BookResolutionHistograms();
+
+ Int_t Counter[9][5] ;
+ Int_t i1,i2,TotalInd = 0 ;
+ for(i1 = 0; i1<9; i1++)
+ for(i2 = 0; i2<5; i2++)
+ Counter[i1][i2] = 0 ;
+
+ Int_t TotalPrimary = 0 ;
+ Int_t TotalRecPart = 0 ;
+ Int_t TotalRPwithPrim = 0 ;
+ Int_t ievent;
+
+ cout << "Start Analysing"<< endl ;
+ for ( ievent=0; ievent<Nevents; ievent++)
+ {
+
+ //========== Event Number>
+ // if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
+ cout << "AnalyzeResolutions > " << "Event is " << ievent << endl ;
+
+ //=========== Connects the various Tree's for evt
+ gAlice->GetEvent(ievent);
+
+ //=========== Gets the Kine TTree
+ gAlice->TreeK()->GetEvent(0) ;
+
+ //=========== Gets the list of Primari Particles
+ TClonesArray * PrimaryList = gAlice->Particles();
+
+ 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 ) {
+ 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() ;
+ gAlice->TreeR()->SetBranchAddress( "PHOSTS", TrackSegmentsList ) ;
+
+ AliPHOSRecParticle::RecParticlesList ** RecParticleList = fPHOS->RecParticles() ;
+ if( (*RecParticleList) )
+ (*RecParticleList)->Clear() ;
+ gAlice->TreeR()->SetBranchAddress( "PHOSRP", RecParticleList ) ;
+
+ //=========== Gets the Reconstraction TTree
+ gAlice->TreeR()->GetEvent(0) ;
+
+ 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 = 5. ;
+ Int_t ClosestPrimary = -1 ;
+
+ Int_t numberofprimaries ;
+ Int_t * listofprimaries = RecParticle->GetPrimaries(numberofprimaries) ;
+ Int_t index ;
+ TParticle * Primary ;
+ Double_t Distance = MinDistance ;
+ for ( index = 0 ; index < numberofprimaries ; index++){
+ Primary = (TParticle*)PrimaryList->UncheckedAt(listofprimaries[index]) ;
+ Int_t ModuleNumber ;
+ Double_t PrimX, PrimZ ;
+ fGeom->ImpactOnEmc(Primary->Theta(), Primary->Phi(), ModuleNumber, PrimX, PrimZ) ;
+ if(ModuleNumberRec == ModuleNumber)
+ Distance = TMath::Sqrt((RecX-PrimX)*(RecX-PrimX)+(RecZ-PrimZ)*(RecZ-PrimZ) ) ;
+ if(MinDistance > Distance)
+ {
+ MinDistance = Distance ;
+ ClosestPrimary = listofprimaries[index] ;
+ }
+ }
+ TotalRecPart++ ;
+
+ 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
+ fhNPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+ }
+
+ if(PrimaryType == -2112){ //neutron ~
+ fhNBarPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+ }
+ if(PrimaryCode == 2){
+ fhChargedPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+ }
+
+ fhAllPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+ fhVeto->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+ fhPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+
+ Counter[5][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();
+ cout << "Resolutions: Analyzed " << Nevents << " event(s)" << endl ;
+ cout << "Resolutions: Total primary " << TotalPrimary << endl ;
+ cout << "Resoluitons: Total reconstracted " << TotalRecPart << endl ;
+ cout << "TotalReconstructed with Primarie " << TotalRPwithPrim << endl ;
+ cout << " Primary: Photon Electron Ch. Hadr. Neutr. Hadr Kaons" << endl ;
+ cout << " Detected as photon " << Counter[0][0] << " " << Counter[0][1] << " " << Counter[0][2] << " " <<Counter[0][3] << " " << Counter[0][4] << endl ;
+ cout << " Detected as electron " << Counter[1][0] << " " << Counter[1][1] << " " << Counter[1][2] << " " <<Counter[1][3] << " " << Counter[1][4] << endl ;
+ cout << " Detected as neutral hadron " << Counter[2][0] << " " << Counter[2][1] << " " << Counter[2][2] << " " <<Counter[2][3] << " " << Counter[2][4] << endl ;
+ cout << " Detected as neutral EM " << Counter[3][0] << " " << Counter[3][1] << " " << Counter[3][2] << " " <<Counter[3][3] << " " << Counter[3][4] << endl ;
+ cout << " Detected as charged hadron " << Counter[4][0] << " " << Counter[4][1] << " " << Counter[4][2] << " " <<Counter[4][3] << " " << Counter[4][4] << endl ;
+ cout << " Detected as gamma-hadron " << Counter[5][0] << " " << Counter[5][1] << " " << Counter[5][2] << " " <<Counter[5][3] << " " << Counter[5][4] << endl ;
+ cout << " Detected as Absurd EM " << Counter[6][0] << " " << Counter[6][1] << " " << Counter[6][2] << " " <<Counter[6][3] << " " << Counter[6][4] << endl ;
+ cout << " Detected as absurd hadron " << Counter[7][0] << " " << Counter[7][1] << " " << Counter[7][2] << " " <<Counter[7][3] << " " << Counter[7][4] << endl ;
+ cout << " Detected as undefined " << Counter[8][0] << " " << Counter[8][1] << " " << Counter[8][2] << " " <<Counter[8][3] << " " << Counter[8][4] << endl ;
+
+ for(i1 = 0; i1<9; i1++)
+ for(i2 = 0; i2<5; i2++)
+ TotalInd+=Counter[i1][i2] ;
+ cout << "Indentified particles " << TotalInd << endl ;
+
+} // endfunction
+
//____________________________________________________________________________
void AliPHOSAnalyze::BookingHistograms()
{
- if (fhEmcDigit )
- delete fhEmcDigit ;
- if (fhVetoDigit )
- delete fhVetoDigit ;
- if (fhConvertorDigit )
- delete fhConvertorDigit ;
- if (fhEmcCluster )
- delete fhEmcCluster ;
- if (fhVetoCluster )
- delete fhVetoCluster ;
- if (fhConvertorCluster )
- delete fhConvertorCluster ;
- if (fhConvertorEmc )
- delete fhConvertorEmc ;
-
+ // Books the histograms where the results of the analysis are stored (to be changed)
+
+ delete fhEmcDigit ;
+ delete fhVetoDigit ;
+ delete fhConvertorDigit ;
+ delete fhEmcCluster ;
+ delete fhVetoCluster ;
+ delete fhConvertorCluster ;
+ delete fhConvertorEmc ;
+
fhEmcDigit = new TH1F("hEmcDigit", "hEmcDigit", 1000, 0. , 25.);
fhVetoDigit = new TH1F("hVetoDigit", "hVetoDigit", 500, 0. , 3.e-5);
fhConvertorDigit = new TH1F("hConvertorDigit","hConvertorDigit", 500, 0. , 3.e-5);
fhVetoCluster = new TH1F("hVetoCluster", "hVetoCluster", 500, 0. , 3.e-5);
fhConvertorCluster = new TH1F("hConvertorCluster","hConvertorCluster",500, 0. , 3.e-5);
fhConvertorEmc = new TH2F("hConvertorEmc", "hConvertorEmc", 200, 1. , 3., 200, 0., 3.e-5);
- fhPhotonEnergy = new TH1F("hPhotonEnergy", "hPhotonEnergy", 1000, 0. , 30.);
- fhElectronEnergy = new TH1F("hElectronEnergy","hElectronEnergy", 1000, 0. , 30.);
- fhNeutralHadronEnergy = new TH1F("hNeutralHadronEnergy", "hNeutralHadronEnergy", 1000, 0. , 30.);
- fhNeutralEMEnergy = new TH1F("hNeutralEMEnergy", "hNeutralEMEnergy", 1000, 0. , 30.);
- fhChargedHadronEnergy = new TH1F("hChargedHadronEnergy", "hChargedHadronEnergy", 1000, 0. , 30.);
- fhPhotonPositionX = new TH1F("hPhotonPositionX","hPhotonPositionX", 500,-80. , 80.);
- fhElectronPositionX = new TH1F("hElectronPositionX","hElectronPositionX",500,-80. , 80.);
- fhNeutralHadronPositionX = new TH1F("hNeutralHadronPositionX","hNeutralHadronPositionX",500,-80. , 80.);
- fhNeutralEMPositionX = new TH1F("hNeutralEMPositionX","hNeutralEMPositionX",500,-80. , 80.);
- fhChargedHadronPositionX = new TH1F("hChargedHadronPositionX","hChargedHadronPositionX",500,-80. , 80.);
- fhPhotonPositionY = new TH1F("hPhotonPositionY","hPhotonPositionY", 500,-80. , 80.);
- fhElectronPositionY = new TH1F("hElectronPositionY","hElectronPositionY",500,-80. , 80.);
- fhNeutralHadronPositionY = new TH1F("hNeutralHadronPositionY","hNeutralHadronPositionY",500,-80. , 80.);
- fhNeutralEMPositionY = new TH1F("hNeutralEMPositionY","hNeutralEMPositionY",500,-80. , 80.);
- fhChargedHadronPositionY = new TH1F("hChargedHadronPositionY","hChargedHadronPositionY",500,-80. , 80.);
+
+}
+//____________________________________________________________________________
+void AliPHOSAnalyze::BookResolutionHistograms()
+{
+ // Books the histograms where the results of the Resolution analysis are stored
+
+// 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 ;
+
+ 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(fhAllEM)
+// delete fhAllEM ;
+// if(fhPhotEM)
+// delete fhPhotEM ;
+// if(fhNEM)
+// delete fhNEM ;
+// if(fhNBarEM)
+// delete fhNBarEM ;
+// if(fhChargedEM)
+// delete fhChargedEM ;
+
+ 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 ;
+
+ 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
+
}
//____________________________________________________________________________
Bool_t AliPHOSAnalyze::Init(Int_t evt)
{
+ // Do a few initializations: open the root file
+ // get the AliRun object
+ // defines the clusterizer, tracksegment maker and particle identifier
+ // sets the associated parameters
Bool_t ok = kTRUE ;
//=========== Get the PHOS object and associated geometry from the file
- fPHOS = (AliPHOSv0 *)gAlice->GetDetector("PHOS") ;
- fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
+ fPHOS = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
+ fGeom = fPHOS->GetGeometry();
+ // fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
+
} // else !ok
} // if fRootFile
//========== Create the Clusterizer
- fClu = new AliPHOSClusterizerv1() ;
- fClu->SetEmcEnergyThreshold(0.025) ;
- fClu->SetEmcClusteringThreshold(0.75) ;
+ fClu = new AliPHOSClusterizerv1() ;
+ fClu->SetEmcEnergyThreshold(0.030) ;
+ fClu->SetEmcClusteringThreshold(0.20) ;
fClu->SetPpsdEnergyThreshold (0.0000002) ;
fClu->SetPpsdClusteringThreshold(0.0000001) ;
fClu->SetLocalMaxCut(0.03) ;
fTrs = new AliPHOSTrackSegmentMakerv1() ;
cout << "AnalyzeOneEvent > using tack segment maker " << fTrs->GetName() << endl ;
+ // fTrs->UnsetUnfoldFlag() ;
//========== Creates the particle identifier
fPID = new AliPHOSPIDv1() ;
cout << "AnalyzeOneEvent > using particle identifier " << fPID->GetName() << endl ;
-
+ //fPID->SetShowerProfileCuts(Float_t l1m, Float_t l1M, Float_t l2m, Float_t l2M) ;
+ fPID->SetShowerProfileCuts(0.7, 2.0 , 0.6 , 1.5) ;
+
//========== Creates the Reconstructioner
- fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
+ fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
+ fRec -> SetDebugReconstruction(kFALSE);
//=========== Connect the various Tree's for evt
//____________________________________________________________________________
void AliPHOSAnalyze::DisplayKineEvent(Int_t evt)
{
+ // Display particles from the Kine Tree in global Alice (theta, phi) coordinates.
+ // One PHOS module at the time.
+ // The particle type can be selected.
+
if (evt == -999)
evt = fEvt ;
sprintf(histoname,"Event %d: Incident particles in module %d", evt, module) ;
Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
- fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, kDegre) ;
+ fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, AliPHOSGeometry::Degre() ) ;
Double_t theta, phi ;
- fGeom->EmcXtalCoverage(theta, phi, kDegre) ;
+ fGeom->EmcXtalCoverage(theta, phi, AliPHOSGeometry::Degre() ) ;
Int_t tdim = (Int_t)( (tM - tm) / theta ) ;
Int_t pdim = (Int_t)( (pM - pm) / phi ) ;
fGeom->ImpactOnEmc(theta, phi, mod, z, x) ;
if ( mod == module ) {
nparticlein++ ;
- histoparticle->Fill(phi*kRADDEG, theta*kRADDEG, particle->Energy() ) ;
+ if (particle->Energy() > fClu->GetEmcClusteringThreshold() )
+ histoparticle->Fill(phi*kRADDEG, theta*kRADDEG, particle->Energy() ) ;
}
}
}
//____________________________________________________________________________
void AliPHOSAnalyze::DisplayRecParticles()
{
+ // Display reconstructed particles in global Alice(theta, phi) coordinates.
+ // One PHOS module at the time.
+ // Click on symbols indicate the reconstructed particle type.
+
if (fEvt == -999) {
- cout << "DisplayRecPoints > Analyze an event first ... (y/n) " ;
+ cout << "DisplayRecParticles > Analyze an event first ... (y/n) " ;
Text_t answer[1] ;
cin >> answer ; cout << answer ;
- if ( answer == "y" )
- AnalyzeOneEvent() ;
+// if ( answer == "y" )
+// AnalyzeOneEvent() ;
}
if (fEvt != -999) {
Int_t module ;
- cout << "DisplayRecPoints > which module (1-5, -1: all) ? " ;
+ cout << "DisplayRecParticles > which module (1-5, -1: all) ? " ;
cin >> module ; cout << module << endl ;
Text_t histoname[80] ;
sprintf(histoname,"Event %d: Reconstructed particles in module %d", fEvt, module) ;
Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
- fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, kDegre) ;
+ fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, AliPHOSGeometry::Degre() ) ;
Double_t theta, phi ;
- fGeom->EmcXtalCoverage(theta, phi, kDegre) ;
+ fGeom->EmcXtalCoverage(theta, phi, AliPHOSGeometry::Degre() ) ;
Int_t tdim = (Int_t)( (tM - tm) / theta ) ;
Int_t pdim = (Int_t)( (pM - pm) / phi ) ;
tm -= theta ;
Text_t canvasname[80] ;
sprintf(canvasname, "Reconstructed particles in PHOSmodule # %d", module) ;
TCanvas * rparticlecanvas = new TCanvas("RparticleCanvas", canvasname, 650, 500) ;
- RecParticlesList * rpl = fPHOS->RecParticles() ;
+ AliPHOSRecParticle::RecParticlesList * rpl = *fPHOS->RecParticles() ;
Int_t nRecParticles = rpl->GetEntries() ;
Int_t nRecParticlesInModule = 0 ;
TIter nextRecPart(rpl) ;
Double_t kRADDEG = 180. / TMath::Pi() ;
while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ;
- if ( ts->GetPHOSMod() == module ) {
+ if ( ts->GetPHOSMod() == module ) {
+ Int_t numberofprimaries = 0 ;
+ Int_t * listofprimaries = 0;
+ rp->GetPrimaries(numberofprimaries) ;
+ cout << "Number of primaries = " << numberofprimaries << endl ;
+ Int_t index ;
+ for ( index = 0 ; index < numberofprimaries ; index++)
+ cout << " primary # " << index << " = " << listofprimaries[index] << endl ;
+
nRecParticlesInModule++ ;
Double_t theta = rp->Theta() * kRADDEG ;
Double_t phi = rp->Phi() * kRADDEG ;
}
}
histoRparticle->Draw("color") ;
+
+ nextRecPart.Reset() ;
+ while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
+ AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ;
+ if ( ts->GetPHOSMod() == module )
+ rp->Draw("P") ;
+ }
+
Text_t text[80] ;
sprintf(text, "reconstructed particles: %d", nRecParticlesInModule) ;
TPaveText * pavetext = new TPaveText(292, 100, 300, 101);
//____________________________________________________________________________
void AliPHOSAnalyze::DisplayRecPoints()
{
+ // Display reconstructed points in local PHOS-module (x, z) coordinates.
+ // One PHOS module at the time.
+ // Click on symbols displays the EMC cluster, or PPSD information.
+
if (fEvt == -999) {
cout << "DisplayRecPoints > Analyze an event first ... (y/n) " ;
Text_t answer[1] ;
cin >> answer ; cout << answer ;
- if ( answer == "y" )
- AnalyzeOneEvent() ;
+// if ( answer == "y" )
+// AnalyzeOneEvent() ;
}
if (fEvt != -999) {
Float_t etot=0.;
Int_t relid[4]; Int_t nDigits = 0 ;
AliPHOSDigit * digit ;
+
+ // Making 2D histogram of the EMC module
while((digit = (AliPHOSDigit *)next()))
{
fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
- if (relid[0] == module)
+ if (relid[0] == module && relid[1] == 0)
{
- nDigits++ ;
energy = fClu->Calibrate(digit->GetAmp()) ;
- etot += energy ;
- fGeom->RelPosInModule(relid,y,z) ;
- if (energy > 0.01 )
+ cout << "Energy is " << energy << " and threshold is " << fClu->GetEmcEnergyThreshold() << endl;
+ if (energy > fClu->GetEmcEnergyThreshold() ){
+ nDigits++ ;
+ etot += energy ;
+ fGeom->RelPosInModule(relid,y,z) ;
hModule->Fill(y, z, energy) ;
+ }
}
}
cout <<"DrawRecPoints > Found in module "
//=========== Cluster in module
- TClonesArray * emcRP = fPHOS->EmcClusters() ;
+ // TClonesArray * emcRP = fPHOS->EmcClusters() ;
+ TObjArray * emcRP = *(fPHOS->EmcRecPoints()) ;
+
etot = 0.;
Int_t totalnClusters = 0 ;
Int_t nClusters = 0 ;
AliPHOSEmcRecPoint * emc ;
while((emc = (AliPHOSEmcRecPoint *)nextemc()))
{
- Int_t numberofprimaries ;
- Int_t * primariesarray = new Int_t[10] ;
- emc->GetPrimaries(numberofprimaries, primariesarray) ;
+ // Int_t numberofprimaries ;
+ // Int_t * primariesarray = new Int_t[10] ;
+ // emc->GetPrimaries(numberofprimaries, primariesarray) ;
totalnClusters++ ;
if ( emc->GetPHOSMod() == module )
{
//=========== Cluster in module PPSD Down
- TClonesArray * ppsdRP = fPHOS->PpsdClusters() ;
+ // TClonesArray * ppsdRP = fPHOS->PpsdClusters() ;
+ TObjArray * ppsdRP = *(fPHOS->PpsdRecPoints() );
+
etot = 0.;
TIter nextPpsd(ppsdRP) ;
AliPHOSPpsdRecPoint * ppsd ;
//=========== Cluster in module PPSD Up
- ppsdRP = fPHOS->PpsdClusters() ;
+ ppsdRP = *(fPHOS->PpsdRecPoints()) ;
+
etot = 0.;
TIter nextPpsdUp(ppsdRP) ;
while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsdUp()))
//____________________________________________________________________________
void AliPHOSAnalyze::DisplayTrackSegments()
{
+ // Display track segments in local PHOS-module (x, z) coordinates.
+ // One PHOS module at the time.
+ // One symbol per PHOS subsystem: EMC, upper PPSD, lower PPSD.
+
if (fEvt == -999) {
cout << "DisplayTrackSegments > Analyze an event first ... (y/n) " ;
Text_t answer[1] ;
cin >> answer ; cout << answer ;
- if ( answer == "y" )
- AnalyzeOneEvent() ;
+// if ( answer == "y" )
+// AnalyzeOneEvent() ;
}
if (fEvt != -999) {
TCanvas * trackcanvas = new TCanvas("TrackSegmentCanvas", canvasname, 650, 500) ;
histotrack->Draw() ;
- TrackSegmentsList * trsegl = fPHOS->TrackSegments() ;
+ AliPHOSTrackSegment::TrackSegmentsList * trsegl = *(fPHOS->TrackSegments()) ;
AliPHOSTrackSegment * trseg ;
Int_t nTrackSegments = trsegl->GetEntries() ;
//____________________________________________________________________________
Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
{
- fRootFile = new TFile(name) ;
+ // Open the root file named "name"
+
+ fRootFile = new TFile(name, "update") ;
return fRootFile->IsOpen() ;
}
//____________________________________________________________________________
-void AliPHOSAnalyze::SavingHistograms()
+void AliPHOSAnalyze::SaveHistograms()
{
- Text_t outputname[80] ;// = fRootFile->GetName();
+ // 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 (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 (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)
+ fhPhotElec->Write() ;
+ if(fhPhotNeuH)
+ fhPhotNeuH->Write() ;
+ if(fhPhotNuEM)
+ fhPhotNuEM->Write() ;
+ if(fhPhotNuEM)
+ fhPhotNuEM->Write() ;
+ if(fhPhotChHa)
+ fhPhotChHa->Write() ;
+ if(fhPhotGaHa)
+ fhPhotGaHa->Write() ;
+ if(fhEnergyCorrelations)
+ fhEnergyCorrelations->Write() ;
+
output.Write();
output.Close();
}
+//____________________________________________________________________________
+Float_t AliPHOSAnalyze::CorrectEnergy(Float_t ERecPart)
+{
+ return ERecPart/0.8783 ;
+}
+
+//____________________________________________________________________________
+void AliPHOSAnalyze::ResetHistograms()
+{
+ fhEnergyCorrelations = 0 ; //Energy correlations between Eloss in Convertor and PPSD(2)
+
+ fhEmcDigit = 0 ; // Histo of digit energies in the Emc
+ fhVetoDigit = 0 ; // Histo of digit energies in the Veto
+ fhConvertorDigit = 0 ; // Histo of digit energies in the Convertor
+ fhEmcCluster = 0 ; // Histo of Cluster energies in Emc
+ fhVetoCluster = 0 ; // Histo of Cluster energies in Veto
+ fhConvertorCluster = 0 ; // Histo of Cluster energies in Convertor
+ fhConvertorEmc = 0 ; // 2d Convertor versus Emc energies
+
+ 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 ;
+ fhPhotEM = 0 ;
+ fhAllEM = 0 ;
+ fhNEM = 0 ;
+ fhNBarEM = 0 ;
+ fhChargedEM = 0 ;
+ fhPhotPPSD = 0 ;
+ fhAllPPSD = 0 ;
+ fhNPPSD = 0 ;
+ fhNBarPPSD = 0 ;
+ fhChargedPPSD = 0 ;
+
+ fhPrimary = 0 ;
+
+ fhPhotPhot = 0 ;
+ fhPhotElec = 0 ;
+ fhPhotNeuH = 0 ;
+ fhPhotNuEM = 0 ;
+ fhPhotChHa = 0 ;
+ fhPhotGaHa = 0 ;
+
+
+}
+
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff