//////////////////////////////////////////////////////////////////////////////
// --- ROOT system ---
+#include "TTree.h"
// --- Standard library ---
{
// ctor
//We do not create objects, because these pointers will be overwritten durin reading from file.
+ fDigits = 0 ;
fEmcRecPoints = 0 ;
fPpsdRecPoints = 0 ;
fTrackSegments = 0 ;
AliPHOS::AliPHOS(const char* name, const char* title): AliDetector(name,title)
{
// ctor
-
- fEmcRecPoints = new AliPHOSRecPoint::RecPointsList(10) ;
- fPpsdRecPoints = new AliPHOSRecPoint::RecPointsList(10) ;
- fTrackSegments = new AliPHOSTrackSegment::TrackSegmentsList("AliPHOSTrackSegment", 10) ;
- fRecParticles = new AliPHOSRecParticle::RecParticlesList("AliPHOSRecParticle", 10) ;
- fDebugLevel = 0;
+
+ fDigits = 0 ;
+ fEmcRecPoints = 0 ;
+ fPpsdRecPoints = 0 ;
+ fTrackSegments = 0 ;
+ fRecParticles = 0 ;
}
//____________________________________________________________________________
gMC->Gstpar(idtmed[715], "STRA",2.) ;
}
-
//____________________________________________________________________________
void AliPHOS::SetTreeAddress()
{
// TBranch *branch;
AliDetector::SetTreeAddress();
- // //Branch address for TreeR: RecParticles
-// TTree *treeR = gAlice->TreeR();
-// if ( treeR && fRecParticles ) {
-// branch = treeR->GetBranch("PHOSRP");
-// if (branch) branch->SetAddress(&fRecParticles) ;
-// }
-// //Branch address for TreeR: TrackSegments
-// if ( treeR && fTrackSegments ) {
-// branch = treeR->GetBranch("PHOSTS");
-// if (branch) branch->SetAddress(&fTrackSegments) ;
-// }
-// //Branch address for TreeR: EmcRecPoint
-// if ( treeR && fEmcRecPoints ) {
-// branch = treeR->GetBranch("PHOSEmcRP");
-// if (branch) branch->SetAddress(&fEmcRecPoints) ;
- // }
+ TBranch * branch ;
+
+ if(fDigits)
+ fDigits->Clear();
+ else
+ fDigits = new TClonesArray("AliPHOSDigit",1) ;
+
+ if (gAlice->TreeD() && fDigits ) {
+ branch = gAlice->TreeD()->GetBranch("PHOS");
+ if (branch) branch->SetAddress(&fDigits) ;
+ }
+
+ TTree *treeR = gAlice->TreeR();
+
+ //Branch address for TreeR: EmcRecPoint
+
+ if(fEmcRecPoints)
+ fEmcRecPoints->Delete();
+ else
+ fEmcRecPoints = new AliPHOSRecPoint::RecPointsList(1) ;
+
+ if ( treeR && fEmcRecPoints ) {
+ branch = treeR->GetBranch("PHOSEmcRP");
+ if (branch) branch->SetAddress(&fEmcRecPoints) ;
+ }
+
+ //Branch address for TreeR: PPSDRecPoint
+ if(fPpsdRecPoints)
+ fPpsdRecPoints->Delete();
+ else
+ fPpsdRecPoints = new AliPHOSRecPoint::RecPointsList(1) ;
+
+ if ( treeR && fPpsdRecPoints ) {
+ branch = treeR->GetBranch("PHOSPpsdRP");
+ if (branch) branch->SetAddress(&fPpsdRecPoints) ;
+ }
+
+ //Branch address for TreeR: TrackSegments
+ if(fTrackSegments)
+ fTrackSegments->Clear() ;
+ else
+ fTrackSegments = new AliPHOSTrackSegment::TrackSegmentsList("AliPHOSTrackSegment", 1) ;
+
+ if ( treeR && fTrackSegments ) {
+ branch = treeR->GetBranch("PHOSTS");
+ if (branch) branch->SetAddress(&fTrackSegments) ;
+ }
+
+ //Branch address for TreeR: RecParticles
+ if(fRecParticles)
+ fRecParticles->Clear() ;
+ else
+ fRecParticles = new AliPHOSRecParticle::RecParticlesList("AliPHOSRecParticle", 1) ;
+
+ if ( treeR && fRecParticles ) {
+ branch = treeR->GetBranch("PHOSRP");
+ if (branch) branch->SetAddress(&fRecParticles) ;
+ }
+
}
return &fEmcRecPoints ;
}
virtual AliPHOSGeometry * GetGeometry() = 0 ;
+ virtual void Hit2Digit(){} // Will convert hits to digits in versions v1 etc
+
virtual Int_t IsVersion(void) const { return -1 ; }
virtual AliPHOSRecPoint::RecPointsList ** PpsdRecPoints() {
// to be redefined when ppsd is present
assert(0==1) ;
return *this ;
}
- void SetDebugLevel(Int_t deb) {fDebugLevel = deb;}
protected:
AliPHOSRecPoint::RecPointsList * fPpsdRecPoints ; // The RecPoints (clusters) list in PPSD (veto)
AliPHOSTrackSegment::TrackSegmentsList * fTrackSegments ;// The TrackSegment list in PHOS
AliPHOSRecParticle::RecParticlesList * fRecParticles ; // The reconstructed particles list in PHOS
- Int_t fDebugLevel ; // Debug level
-
ClassDef(AliPHOS,2) // Photon Spectrometer Detector (base class)
#include "AliPHOSTrackSegment.h"
#include "AliPHOSRecParticle.h"
#include "AliPHOSIndexToObject.h"
+#include "AliPHOSHit.h"
#include "AliPHOSCPVHit.h"
#include "AliPHOSCpvRecPoint.h"
if(fTrs) {delete fTrs ; fTrs =0 ;}
}
-
//____________________________________________________________________________
-void AliPHOSAnalyze::ActivePPSD(Int_t Nevents=1){
-
- fhEnergyCorrelations = new TH2F("hEnergyCorrelations","hEnergyCorrelations",40, 0., 0.15, 30, 0., 3.e-5);
+void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod){
+ //Draws pimary particles and reconstructed
+ //digits, RecPoints, RecPartices etc
+ //for event Nevent in the module Nmod.
+
+ TH2F * digitOccupancy = new TH2F("digitOccupancy","EMC digits", 64,-71.,71.,64,-71.,71.);
+ TH2F * emcOccupancy = new TH2F("emcOccupancy","EMC RecPoints",64,-71.,71.,64,-71.,71.);
+ TH2F * ppsdUp = new TH2F("ppsdUp","PPSD Up digits", 128,-71.,71.,128,-71.,71.) ;
+ TH2F * ppsdUpCl = new TH2F("ppsdUpCl","PPSD Up RecPoints",128,-71.,71.,128,-71.,71.) ;
+ TH2F * ppsdLow = new TH2F("ppsdLow","PPSD Low digits", 128,-71.,71.,128,-71.,71.) ;
+ TH2F * ppsdLowCl = new TH2F("ppsdLowCl","PPSD Low RecPoints",128,-71.,71.,128,-71.,71.) ;
+ TH2F * nbar = new TH2F("nbar","Primary nbar", 64,-71.,71.,64,-71.,71.);
+ TH2F * phot = new TH2F("phot","Primary Photon", 64,-71.,71.,64,-71.,71.);
+ TH2F * charg = new TH2F("charg","Primary charged",64,-71.,71.,64,-71.,71.);
+ TH2F * recPhot = new TH2F("recPhot","RecParticles with primary Photon",64,-71.,71.,64,-71.,71.);
+ TH2F * recNbar = new TH2F("recNbar","RecParticles with primary Nbar", 64,-71.,71.,64,-71.,71.);
+
//========== 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;
+ fClu->SetEmcEnergyThreshold(0.05) ;
+ fClu->SetEmcClusteringThreshold(0.20) ;
+ fClu->SetPpsdEnergyThreshold (0.0000002) ;
+ fClu->SetPpsdClusteringThreshold(0.0000001) ;
+ fClu->SetLocalMaxCut(0.03) ;
+ fClu->SetCalibrationParameters(0., 0.00000001) ;
- 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") ;
-}
-
+ gAlice->GetEvent(Nevent);
+
+ TClonesArray * primaryList = gAlice->Particles();
+
+ TParticle * primary ;
+ Int_t iPrimary ;
+ for ( iPrimary = 0 ; iPrimary < primaryList->GetEntries() ; iPrimary++)
+ {
+ primary = (TParticle*)primaryList->At(iPrimary) ;
+ Int_t primaryType = primary->GetPdgCode() ;
+ if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212) ) {
+ Int_t moduleNumber ;
+ Double_t primX, primZ ;
+ fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
+ if(moduleNumber==Nmod)
+ charg->Fill(primZ,primX,primary->Energy()) ;
+ }
+ if( primaryType == 22 ) {
+ Int_t moduleNumber ;
+ Double_t primX, primZ ;
+ fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
+ if(moduleNumber==Nmod)
+ phot->Fill(primZ,primX,primary->Energy()) ;
+ }
+ else{
+ if( primaryType == -2112 ) {
+ Int_t moduleNumber ;
+ Double_t primX, primZ ;
+ fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
+ if(moduleNumber==Nmod)
+ nbar->Fill(primZ,primX,primary->Energy()) ;
+ }
+ }
+ }
-//____________________________________________________________________________
-void AliPHOSAnalyze::AnalyzeManyEvents(Int_t Nevents, Int_t module)
-{
- // analyzes Nevents events in a single PHOS module
- // Events should be reconstructed by Reconstruct()
+ fPHOS->SetTreeAddress() ;
- if ( fRootFile == 0 )
- cout << "AnalyzeManyEvents > " << "Root File not openned" << endl ;
- else
+ gAlice->TreeD()->GetEvent(0) ;
+ gAlice->TreeR()->GetEvent(0) ;
+
+ TObjArray ** emcRecPoints = fPHOS->EmcRecPoints() ;
+ TObjArray ** ppsdRecPoints = fPHOS->PpsdRecPoints() ;
+ TClonesArray ** recParticleList = fPHOS->RecParticles() ;
+
+ Int_t iDigit ;
+ AliPHOSDigit * digit ;
+
+ for(iDigit = 0; iDigit < fPHOS->Digits()->GetEntries(); iDigit++)
{
- //========== Booking Histograms
- cout << "AnalyzeManyEvents > " << "Booking Histograms" << endl ;
- BookingHistograms();
-
- Int_t ievent;
- Int_t relid[4] ;
- AliPHOSDigit * digit ;
- AliPHOSEmcRecPoint * emc ;
- AliPHOSPpsdRecPoint * ppsd ;
- // AliPHOSTrackSegment * tracksegment ;
- AliPHOSRecParticle * recparticle;
-
- for ( ievent=0; ievent<Nevents; ievent++)
- {
- //========== 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() ) )
- {
- fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
- if (fClu->IsInEmc(digit)) fhEmcDigit->Fill(fClu->Calibrate(digit->GetAmp())) ;
- else
- {
- if (relid[1]<17) fhVetoDigit->Fill(fClu->Calibrate(digit->GetAmp()));
- if (relid[1]>16) fhConvertorDigit->Fill(fClu->Calibrate(digit->GetAmp()));
- }
- }
+ digit = (AliPHOSDigit *) fPHOS->Digits()->At(iDigit) ;
+ Int_t relid[4];
+ fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ Float_t x,z ;
+ fGeom->RelPosInModule(relid,x,z) ;
+ Float_t e = fClu->Calibrate(digit->GetAmp()) ;
+ if(relid[0]==Nmod){
+ if(relid[1]==0) //EMC
+ digitOccupancy->Fill(x,z,e) ;
+ if((relid[1]>0)&&(relid[1]<17))
+ ppsdUp->Fill(x,z,e) ;
+ if(relid[1]>16)
+ ppsdLow->Fill(x,z,e) ;
+ }
+ }
+
+ Int_t irecp ;
+ TVector3 pos ;
+
+ for(irecp = 0; irecp < (*emcRecPoints)->GetEntries() ; irecp ++){
+ AliPHOSEmcRecPoint * emc= (AliPHOSEmcRecPoint*)(*emcRecPoints)->At(irecp) ;
+ if(emc->GetPHOSMod()==Nmod){
+ emc->GetLocalPosition(pos) ;
+ emcOccupancy->Fill(pos.X(),pos.Z(),emc->GetEnergy());
+ }
+ }
+
+ for(irecp = 0; irecp < (*ppsdRecPoints)->GetEntries() ; irecp ++){
+ AliPHOSPpsdRecPoint * ppsd= (AliPHOSPpsdRecPoint *)(*ppsdRecPoints)->At(irecp) ;
+ if(ppsd->GetPHOSMod()==Nmod){
+ ppsd->GetLocalPosition(pos) ;
+ if(ppsd->GetUp())
+ ppsdUpCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy());
+ else
+ ppsdLowCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy());
+ }
+ }
+
+ AliPHOSRecParticle * recParticle ;
+ Int_t iRecParticle ;
+ for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ )
+ {
+ recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
+
+ Int_t moduleNumberRec ;
+ Double_t recX, recZ ;
+ fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
+ if(moduleNumberRec == Nmod){
+
+ 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 ;
-
- //=========== Cluster in module
- TIter nextEmc(*fPHOS->EmcRecPoints() ) ;
- while((emc = (AliPHOSEmcRecPoint *)nextEmc()))
+ for ( index = 0 ; index < numberofprimaries ; index++){
+ primary = (TParticle*)primaryList->At(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)
{
- if ( emc->GetPHOSMod() == module )
- {
- fhEmcCluster->Fill( emc->GetTotalEnergy() );
- TIter nextPpsd( *fPHOS->PpsdRecPoints()) ;
- while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
- {
- if ( ppsd->GetPHOSMod() == module )
- {
- if (!ppsd->GetUp()) fhConvertorEmc->Fill(emc->GetTotalEnergy(),ppsd->GetTotalEnergy()) ;
- }
- }
- }
- }
-
- //=========== Cluster in module PPSD Down
- TIter nextPpsd(*fPHOS->PpsdRecPoints() ) ;
- while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
- {
- if ( ppsd->GetPHOSMod() == module )
- {
- if (!ppsd->GetUp()) fhConvertorCluster->Fill(ppsd->GetTotalEnergy()) ;
- if (ppsd->GetUp()) fhVetoCluster ->Fill(ppsd->GetTotalEnergy()) ;
- }
+ minDistance = distance ;
+ closestPrimary = listofprimaries[index] ;
}
-
- //========== TRackSegments in the event
- TIter nextRecParticle(*fPHOS->RecParticles() ) ;
- while((recparticle = (AliPHOSRecParticle *)nextRecParticle()))
- {
- if ( recparticle->GetPHOSTrackSegment()->GetPHOSMod() == module )
- {
- 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 ;
- }
- }
- } // endfor
- SaveHistograms();
- } // endif
-} // endfunction
-
+ }
+
+ if(closestPrimary >=0 ){
+
+ Int_t primaryType = ((TParticle *)primaryList->At(closestPrimary))->GetPdgCode() ;
+
+ if(primaryType==22)
+ recPhot->Fill(recZ,recX,recParticle->Energy()) ;
+ else
+ if(primaryType==-2112)
+ recNbar->Fill(recZ,recX,recParticle->Energy()) ;
+ }
+ }
+ }
+
+
+ digitOccupancy->Draw("box") ;
+ emcOccupancy->SetLineColor(2) ;
+ emcOccupancy->Draw("boxsame") ;
+ ppsdUp->SetLineColor(3) ;
+ ppsdUp->Draw("boxsame") ;
+ ppsdLow->SetLineColor(4) ;
+ ppsdLow->Draw("boxsame") ;
+ phot->SetLineColor(8) ;
+ phot->Draw("boxsame") ;
+ nbar->SetLineColor(6) ;
+ nbar->Draw("boxsame") ;
+
+}
//____________________________________________________________________________
- void AliPHOSAnalyze::Reconstruct(Int_t Nevents,Int_t FirstEvent )
+ void AliPHOSAnalyze::Reconstruct(Int_t nevents,Int_t firstEvent )
{
- // Performs reconstruction of EMC and CPV (GPS2 or IHEP)
+ // Performs reconstruction of EMC and CPV (GPS2, IHEP or MIXT)
// for events from FirstEvent to Nevents
Int_t ievent ;
- for ( ievent=FirstEvent; ievent<Nevents; ievent++) {
- if (ievent==FirstEvent) {
+ for ( ievent=firstEvent; ievent<nevents; ievent++) {
+ if (ievent==firstEvent) {
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) {
+ if (strcmp(fGeom->GetName(),"GPS2") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0) {
fClu->SetPpsdEnergyThreshold (0.0000002) ;
fClu->SetPpsdClusteringThreshold(0.0000001) ;
}
- else if (strcmp(fGeom->GetName(),"IHEP") == 0) {
+ else if (strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0) {
fClu->SetLocalMaxCutCPV(0.03) ;
fClu->SetLogWeightCutCPV(4.0) ;
- fClu->SetPpsdEnergyThreshold (0.09) ;
+ fClu->SetCpvEnergyThreshold(0.09) ;
}
fClu->SetCalibrationParameters(0., 0.00000001) ;
// fTrs->UnsetUnfoldFlag() ;
//========== Creates the particle identifier for GPS2 only
- if (strcmp(fGeom->GetName(),"GPS2") == 0) {
+ if (strcmp(fGeom->GetName(),"GPS2") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0) {
fPID = new AliPHOSPIDv1() ;
fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ;
}
cout << "======= Analyze ======> Event " << ievent+1 << endl ;
//=========== Connects the various Tree's for evt
- gAlice->GetEvent(ievent);
-
- //=========== Gets the Digit TTree
- gAlice->TreeD()->GetEvent(0) ;
-
+ Int_t tracks = gAlice->GetEvent(ievent);
+
+ fPHOS->Hit2Digit(tracks) ;
+
//=========== Do the reconstruction
- fPHOS->Reconstruction(fRec);
- }
+ fPHOS->Reconstruction(fRec);
+ }
+
if(fClu) {delete fClu ; fClu =0 ;}
if(fPID) {delete fPID ; fPID =0 ;}
if(fRec) {delete fRec ; fRec =0 ;}
Int_t nOfModules = fGeom->GetNModules();
TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
- Int_t iModule ;
+ Int_t iModule = 0;
for (iModule=0; iModule < nOfModules; iModule++)
hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100);
// Get the Hits Tree for the Primary track itrack
gAlice->ResetHits();
gAlice->TreeH()->GetEvent(itrack);
- for (iModule=0; iModule < nOfModules; iModule++) {
+ for (Int_t iModule=0; iModule < nOfModules; iModule++) {
cpvModule = fPHOS->GetCPVModule(iModule);
cpvHits = cpvModule.Hits();
nCPVhits = cpvHits->GetEntriesFast();
void AliPHOSAnalyze::InvariantMass(Int_t Nevents )
{
// Calculates Real and Mixed invariant mass distributions
- const Int_t NMixedEvents = 4 ; //# of events used for calculation of 'mixed' distribution
- Int_t MixedLoops = (Int_t )TMath::Ceil(Nevents/NMixedEvents) ;
+
+ 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 * hMixedPhot= new TH2D("hMixedPhot","Mixed for kPhoton particles",250,0.,1.,40,0.,4.) ;
Int_t ievent;
- Int_t EventInMixedLoop ;
+ Int_t eventInMixedLoop ;
- Int_t NRecParticles[NMixedEvents] ;
+ Int_t nRecParticles[nMixedEvents] ;
- AliPHOSRecParticle::RecParticlesList * AllRecParticleList = new TClonesArray("AliPHOSRecParticle", NMixedEvents*1000) ;
+ AliPHOSRecParticle::RecParticlesList * allRecParticleList = new TClonesArray("AliPHOSRecParticle", nMixedEvents*1000) ;
- for(EventInMixedLoop = 0; EventInMixedLoop < MixedLoops; EventInMixedLoop++ ){
+ for(eventInMixedLoop = 0; eventInMixedLoop < mixedLoops; eventInMixedLoop++ ){
Int_t iRecPhot = 0 ;
- for ( ievent=0; ievent < NMixedEvents; ievent++){
+ for ( ievent=0; ievent < nMixedEvents; ievent++){
- Int_t AbsEventNumber = EventInMixedLoop*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) ;
-
+ gAlice->GetEvent(absEventNumber);
+
//========== Creating branches ===================================
+ fPHOS->SetTreeAddress() ;
- AliPHOSRecParticle::RecParticlesList ** RecParticleList = fPHOS->RecParticles() ;
- if( (*RecParticleList) )
- (*RecParticleList)->Clear() ;
- gAlice->TreeR()->SetBranchAddress( "PHOSRP", RecParticleList ) ;
-
- //=========== Gets the Reconstraction TTree
+ gAlice->TreeD()->GetEvent(0) ;
gAlice->TreeR()->GetEvent(0) ;
- AliPHOSRecParticle * RecParticle ;
+ TClonesArray ** recParticleList = fPHOS->RecParticles() ;
+
+
+ AliPHOSRecParticle * recParticle ;
Int_t iRecParticle ;
- for(iRecParticle = 0; iRecParticle < (*RecParticleList)->GetEntries() ;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) ;
+ recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
+ if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
+ (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)){
+ new( (*allRecParticleList)[iRecPhot] ) AliPHOSRecParticle(*recParticle) ;
iRecPhot++;
}
}
- NRecParticles[ievent] = iRecPhot-1 ;
+ nRecParticles[ievent] = iRecPhot-1 ;
}
//Now calculate invariant mass:
Int_t irp1,irp2 ;
- Int_t NCurEvent = 0 ;
+ Int_t nCurEvent = 0 ;
- for(irp1 = 0; irp1 < AllRecParticleList->GetEntries()-1; irp1++){
- AliPHOSRecParticle * rp1 = (AliPHOSRecParticle *)AllRecParticleList->At(irp1) ;
+ 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) ;
+ 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())-
+ 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);
+ 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()));
+ 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(irp1 > nRecParticles[nCurEvent])
+ nCurEvent++;
- if(irp2 <= NRecParticles[NCurEvent]){ //'Real' event
- hRealEM->Fill(InvMass,Pt);
+ if(irp2 <= nRecParticles[nCurEvent]){ //'Real' event
+ hRealEM->Fill(invMass,pt);
if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
- hRealPhot->Fill(InvMass,Pt);
+ hRealPhot->Fill(invMass,pt);
}
else{
- hMixedEM->Fill(InvMass,Pt);
+ hMixedEM->Fill(invMass,pt);
if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
- hMixedPhot->Fill(InvMass,Pt);
+ hMixedPhot->Fill(invMass,pt);
} //real-mixed
} //loop over second rp
}//loop over first rp
- AllRecParticleList->Delete() ;
+ allRecParticleList->Delete() ;
} //Loop over events
- delete AllRecParticleList ;
+ delete allRecParticleList ;
//writing output
TFile output("invmass.root","RECREATE");
}
+//____________________________________________________________________________
+ void AliPHOSAnalyze::ReadAndPrintEMC(Int_t EvFirst=0, Int_t EvLast=0)
+{
+ //
+ // Read and print generated and reconstructed hits in EMC
+ // for events from EvFirst to Nevent.
+ // If only EvFirst is defined, print only this one event.
+ // Author: Yuri Kharlov
+ // 24 November 2000
+ //
+
+ if (EvFirst!=0 && EvLast==0) EvLast=EvFirst;
+ Int_t ievent;
+ for (ievent=EvFirst; ievent<=EvLast; ievent++) {
+
+ //========== Event Number>
+ cout << endl << "==== ReadAndPrintEMC ====> Event is " << ievent+1 << endl ;
+
+ //=========== Connects the various Tree's for evt
+ Int_t ntracks = gAlice->GetEvent(ievent);
+ fPHOS->SetTreeAddress() ;
+
+ gAlice->TreeD()->GetEvent(0) ;
+ gAlice->TreeR()->GetEvent(0) ;
+
+ // Loop over reconstructed particles
+
+ TClonesArray ** recParticleList = fPHOS->RecParticles() ;
+ AliPHOSRecParticle * recParticle ;
+ Int_t iRecParticle ;
+ Int_t *primList;
+ Int_t nPrimary;
+ for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ ) {
+ recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
+ Float_t recE = recParticle->Energy();
+ primList = recParticle->GetPrimaries(nPrimary);
+ Int_t moduleNumberRec ;
+ Double_t recX, recZ ;
+ fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
+ printf("Rec point: module %d, (X,Z) = (%8.4f,%8.4f) cm, E = %.3f GeV, primary = %d\n",
+ moduleNumberRec,recX,recZ,recE,*primList);
+ }
+
+ // Read and print EMC hits from EMCn branches
+
+ AliPHOSCPVModule emcModule;
+ TClonesArray *emcHits;
+ Int_t nEMChits;
+ AliPHOSCPVHit *emcHit;
+ TLorentzVector p;
+ Float_t xgen, zgen;
+ Int_t ipart, primary;
+ Int_t nGenHits = 0;
+ for (Int_t itrack=0; itrack<ntracks; itrack++) {
+ //=========== Get the Hits Tree for the Primary track itrack
+ gAlice->ResetHits();
+ gAlice->TreeH()->GetEvent(itrack);
+ Int_t iModule = 0 ;
+ for (iModule=0; iModule < fGeom->GetNModules(); iModule++) {
+ emcModule = fPHOS->GetEMCModule(iModule);
+ emcHits = emcModule.Hits();
+ nEMChits = emcHits->GetEntriesFast();
+ for (Int_t ihit=0; ihit<nEMChits; ihit++) {
+ nGenHits++;
+ emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
+ p = emcHit->GetMomentum();
+ xgen = emcHit->X();
+ zgen = emcHit->Y();
+ ipart = emcHit->GetIpart();
+ primary= emcHit->GetTrack();
+ printf("EMC hit A: module %d, ",iModule+1);
+ printf(" p = (%f .4, %f .4, %f .4, %f .4) GeV,\n",
+ p.Px(),p.Py(),p.Pz(),p.Energy());
+ printf(" (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
+ xgen,zgen,ipart,primary);
+ }
+ }
+ }
+
+// // Read and print EMC hits from PHOS branch
+
+// for (Int_t itrack=0; itrack<ntracks; itrack++) {
+// //=========== Get the Hits Tree for the Primary track itrack
+// gAlice->ResetHits();
+// gAlice->TreeH()->GetEvent(itrack);
+// TClonesArray *hits = fPHOS->Hits();
+// AliPHOSHit *hit ;
+// Int_t ihit;
+// for ( ihit = 0 ; ihit < hits->GetEntries() ; ihit++ ) {
+// hit = (AliPHOSHit*)hits->At(ihit) ;
+// Float_t hitXYZ[3];
+// hitXYZ[0] = hit->X();
+// hitXYZ[1] = hit->Y();
+// hitXYZ[2] = hit->Z();
+// ipart = hit->GetPid();
+// primary = hit->GetPrimary();
+// Int_t absId = hit->GetId();
+// Int_t relId[4];
+// fGeom->AbsToRelNumbering(absId, relId) ;
+// Int_t module = relId[0];
+// if (relId[1]==0 && !(hitXYZ[0]==0 && hitXYZ[2]==0))
+// printf("EMC hit B: module %d, (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
+// module,hitXYZ[0],hitXYZ[2],ipart,primary);
+// }
+// }
+
+ }
+}
+
+//____________________________________________________________________________
+ void AliPHOSAnalyze::AnalyzeEMC(Int_t Nevents)
+{
+ //
+ // Read generated and reconstructed hits in EMC for Nevents events.
+ // Plots the coordinate and energy resolution histograms.
+ // Coordinate resolution is a difference between the reconstructed
+ // coordinate and the exact coordinate on the face of the PHOS
+ // Author: Yuri Kharlov
+ // 27 November 2000
+ //
+
+ // Book histograms
+
+ TH1F *hDx1 = new TH1F("hDx1" ,"EMC x-resolution", 100,-5. , 5.);
+ TH1F *hDz1 = new TH1F("hDz1" ,"EMC z-resolution", 100,-5. , 5.);
+ TH1F *hDE1 = new TH1F("hDE1" ,"EMC E-resolution", 100,-2. , 2.);
+
+ TH2F *hDx2 = new TH2F("hDx2" ,"EMC x-resolution", 100, 0., 10., 100,-5. , 5.);
+ TH2F *hDz2 = new TH2F("hDz2" ,"EMC z-resolution", 100, 0., 10., 100,-5. , 5.);
+ TH2F *hDE2 = new TH2F("hDE2" ,"EMC E-resolution", 100, 0., 10., 100, 0. , 5.);
+
+ cout << "Start EMC 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 << "==== AnalyzeEMC ====> Event is " << ievent+1 << endl ;
+
+ //=========== Connects the various Tree's for evt
+ Int_t ntracks = gAlice->GetEvent(ievent);
+
+ fPHOS->SetTreeAddress() ;
+
+ gAlice->TreeD()->GetEvent(0) ;
+ gAlice->TreeR()->GetEvent(0) ;
+
+ // Create and fill arrays of hits for each EMC module
+
+ Int_t nOfModules = fGeom->GetNModules();
+ TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
+ Int_t iModule;
+ for (iModule=0; iModule < nOfModules; iModule++)
+ hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100);
+
+ AliPHOSCPVModule emcModule;
+ TClonesArray *emcHits;
+ Int_t nEMChits;
+ AliPHOSCPVHit *emcHit;
+
+ // First go through all primary tracks and fill the arrays
+ // of hits per each EMC module
+
+ for (Int_t itrack=0; itrack<ntracks; itrack++) {
+ // Get the Hits Tree for the Primary track itrack
+ gAlice->ResetHits();
+ gAlice->TreeH()->GetEvent(itrack);
+ for (Int_t iModule=0; iModule < nOfModules; iModule++) {
+ emcModule = fPHOS->GetEMCModule(iModule);
+ emcHits = emcModule.Hits();
+ nEMChits = emcHits->GetEntriesFast();
+ for (Int_t ihit=0; ihit<nEMChits; ihit++) {
+ emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
+ TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
+ new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*emcHit);
+ }
+ emcModule.Clear();
+ }
+ }
+
+ // Loop over reconstructed particles
+
+ TClonesArray ** recParticleList = fPHOS->RecParticles() ;
+ AliPHOSRecParticle * recParticle ;
+ Int_t nEMCrecs = (*recParticleList)->GetEntries();
+ if (nEMCrecs == 1) {
+ recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(0) ;
+ Float_t recE = recParticle->Energy();
+ Int_t phosModule;
+ Double_t recX, recZ ;
+ fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), phosModule, recX, recZ) ;
+
+ // for this rec.point take the hit list in the same PHOS module
+
+ emcHits = hitsPerModule[phosModule-1];
+ Int_t nEMChits = emcHits->GetEntriesFast();
+ if (nEMChits == 1) {
+ Float_t genX, genZ, genE;
+ for (Int_t ihit=0; ihit<nEMChits; ihit++) {
+ emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
+ genX = emcHit->X();
+ genZ = emcHit->Y();
+ genE = emcHit->GetMomentum().E();
+ }
+ Float_t dx = recX - genX;
+ Float_t dz = recZ - genZ;
+ Float_t de = recE - genE;
+ hDx1 ->Fill(dx);
+ hDz1 ->Fill(dz);
+ hDE1 ->Fill(de);
+ hDx2 ->Fill(genE,dx);
+ hDz2 ->Fill(genE,dz);
+ hDE2 ->Fill(genE,recE);
+ }
+ }
+ delete [] hitsPerModule;
+ }
+ // Save histograms
+
+ Text_t outputname[80] ;
+ sprintf(outputname,"%s.analyzed",fRootFile->GetName());
+ TFile output(outputname,"RECREATE");
+ output.cd();
+
+ hDx1 ->Write() ;
+ hDz1 ->Write() ;
+ hDE1 ->Write() ;
+ hDx2 ->Write() ;
+ hDz2 ->Write() ;
+ hDE2 ->Write() ;
+
+ // Plot histograms
+
+ TCanvas *emcCanvas = new TCanvas("EMC","EMC analysis",20,20,700,300);
+ gStyle->SetOptStat(111111);
+ gStyle->SetOptFit(1);
+ gStyle->SetOptDate(1);
+ emcCanvas->Divide(3,1);
+
+ emcCanvas->cd(1);
+ gPad->SetFillColor(10);
+ hDx1->SetFillColor(16);
+ hDx1->Draw();
+
+ emcCanvas->cd(2);
+ gPad->SetFillColor(10);
+ hDz1->SetFillColor(16);
+ hDz1->Draw();
+
+ emcCanvas->cd(3);
+ gPad->SetFillColor(10);
+ hDE1->SetFillColor(16);
+ hDE1->Draw();
+
+ emcCanvas->Print("EMC.ps");
+
+}
+
//____________________________________________________________________________
void AliPHOSAnalyze::AnalyzeResolutions(Int_t Nevents )
{
cout << "AnalyzeResolutions > " << "Booking Histograms" << endl ;
BookResolutionHistograms();
- Int_t Counter[9][5] ;
- Int_t i1,i2,TotalInd = 0 ;
+ 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 ;
+ counter[i1][i2] = 0 ;
- Int_t TotalPrimary = 0 ;
- Int_t TotalRecPart = 0 ;
- Int_t TotalRPwithPrim = 0 ;
+ Int_t totalPrimary = 0 ;
+ Int_t totalRecPart = 0 ;
+ Int_t totalRPwithPrim = 0 ;
Int_t ievent;
cout << "Start Analysing"<< endl ;
gAlice->TreeK()->GetEvent(0) ;
//=========== Gets the list of Primari Particles
- TClonesArray * PrimaryList = gAlice->Particles();
+ TClonesArray * primaryList = gAlice->Particles();
- TParticle * Primary ;
+ TParticle * primary ;
Int_t iPrimary ;
- for ( iPrimary = 0 ; iPrimary < PrimaryList->GetEntries() ; 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++ ;
+ 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 ) ;
+ fPHOS->SetTreeAddress() ;
- //=========== Gets the Reconstraction TTree
+ gAlice->TreeD()->GetEvent(0) ;
gAlice->TreeR()->GetEvent(0) ;
- AliPHOSRecParticle * RecParticle ;
+ TClonesArray ** recParticleList = fPHOS->RecParticles() ;
+
+ AliPHOSRecParticle * recParticle ;
Int_t iRecParticle ;
- for(iRecParticle = 0; iRecParticle < (*RecParticleList)->GetEntries() ;iRecParticle++ )
+ for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ )
{
- RecParticle = (AliPHOSRecParticle *) (*RecParticleList)->At(iRecParticle) ;
- fhAllRP->Fill(CorrectEnergy(RecParticle->Energy())) ;
+ 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) ;
+ Int_t moduleNumberRec ;
+ Double_t recX, recZ ;
+ fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
- Double_t MinDistance = 2. ;
- Int_t ClosestPrimary = -1 ;
+ Double_t minDistance = 100. ;
+ Int_t closestPrimary = -1 ;
Int_t numberofprimaries ;
- Int_t * listofprimaries = RecParticle->GetPrimaries(numberofprimaries) ;
+ Int_t * listofprimaries = recParticle->GetPrimaries(numberofprimaries) ;
Int_t index ;
- TParticle * Primary ;
- Double_t Distance = MinDistance ;
+ TParticle * primary ;
+ Double_t distance = minDistance ;
+ Double_t dX, dZ;
+ Double dXmin = 0.;
+ Double dZmin = 0. ;
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] ;
+ 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) {
+ dX = recX - primX;
+ dZ = recZ - primZ;
+ distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
+ if(minDistance > distance) {
+ minDistance = distance ;
+ dXmin = dX;
+ dZmin = dZ;
+ closestPrimary = listofprimaries[index] ;
}
+ }
}
- TotalRecPart++ ;
+ totalRecPart++ ;
- if(ClosestPrimary >=0 ){
- TotalRPwithPrim++;
+ if(closestPrimary >=0 ){
+ totalRPwithPrim++;
- Int_t PrimaryType = ((TParticle *)PrimaryList->At(ClosestPrimary))->GetPdgCode() ;
-// TParticlePDG* PDGparticle = ((TParticle *)PrimaryList->At(ClosestPrimary))->GetPDG();
+ Int_t primaryType = ((TParticle *)primaryList->At(closestPrimary))->GetPdgCode() ;
+// TParticlePDG* pDGparticle = ((TParticle *)primaryList->At(closestPrimary))->GetPDG();
// Double_t charge = PDGparticle->Charge() ;
// if(charge)
-// cout <<"Primary " <<PrimaryType << " E " << ((TParticle *)PrimaryList->At(ClosestPrimary))->Energy() << endl ;
- Int_t PrimaryCode ;
- switch(PrimaryType)
+// cout <<"Primary " <<primaryType << " E " << ((TParticle *)primaryList->At(closestPrimary))->Energy() << endl ;
+ 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) ;
+ primaryCode = 0; //Photon
+ fhAllEnergy ->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), recParticle->Energy()) ;
+ fhAllPosition ->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), minDistance) ;
+ fhAllPositionX->Fill(dXmin);
+ fhAllPositionZ->Fill(dZmin);
break;
case 11 :
- PrimaryCode = 1; //Electron
+ primaryCode = 1; //Electron
break;
case -11 :
- PrimaryCode = 1; //positron
+ primaryCode = 1; //positron
break;
case 321 :
- PrimaryCode = 4; //K+
+ primaryCode = 4; //K+
break;
case -321 :
- PrimaryCode = 4; //K-
+ primaryCode = 4; //K-
break;
case 310 :
- PrimaryCode = 4; //K0s
+ primaryCode = 4; //K0s
break;
case 130 :
- PrimaryCode = 4; //K0l
+ primaryCode = 4; //K0l
break;
case 211 :
- PrimaryCode = 2; //K0l
+ primaryCode = 2; //K0l
break;
case -211 :
- PrimaryCode = 2; //K0l
+ primaryCode = 2; //K0l
break;
case 2212 :
- PrimaryCode = 2; //K0l
+ primaryCode = 2; //K0l
break;
case -2212 :
- PrimaryCode = 2; //K0l
+ primaryCode = 2; //K0l
break;
default:
- PrimaryCode = 3; //ELSE
+ primaryCode = 3; //ELSE
break;
}
- switch(RecParticle->GetType())
+ 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() ) ;
+ 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) ;
+ 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()) ) ;
+ fhPhotReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+ fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+ fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
- fhPhotPhot->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
+ 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(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()) ) ;
+ 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]++;
+ 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 == 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
+ 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(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()) ) ;
+ 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]++;
+ 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()) ) ;
+ if(primaryType == 22)
+ fhPhotNeuH->Fill(CorrectEnergy(recParticle->Energy()) ) ;
- fhVeto->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
- Counter[2][PrimaryCode]++;
+ 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 ) ;
+ 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()) ) ;
+ fhPhotNuEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+ fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
}
- if(PrimaryType == 2112) //neutron
- fhNEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+ if(primaryType == 2112) //neutron
+ fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
- if(PrimaryType == -2112) //neutron ~
- fhNBarEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+ if(primaryType == -2112) //neutron ~
+ fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
- if(PrimaryCode == 2)
- fhChargedEM->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()) ) ;
+ fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+ fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
+ fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
- Counter[3][PrimaryCode]++;
+ counter[3][primaryCode]++;
break ;
case AliPHOSFastRecParticle::kCHARGEDHA:
- if(PrimaryType == 22) //photon
- fhPhotChHa->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+ if(primaryType == 22) //photon
+ fhPhotChHa->Fill(CorrectEnergy(recParticle->Energy()) ) ;
- Counter[4][PrimaryCode]++ ;
+ 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 == 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
+ fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
}
- if(PrimaryType == -2112){ //neutron ~
- fhNBarPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
+ if(primaryType == -2112){ //neutron ~
+ fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
}
- if(PrimaryCode == 2){
- fhChargedPPSD->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]++ ;
+ 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()) ) ;
+ counter[6][primaryCode]++ ;
+ fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
break;
case AliPHOSFastRecParticle::kABSURDHA:
- Counter[7][PrimaryCode]++ ;
+ counter[7][primaryCode]++ ;
break;
default:
- Counter[8][PrimaryCode]++ ;
+ 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 << "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 ;
+ 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 ;
+ totalInd+=counter[i1][i2] ;
+ cout << "Indentified particles " << totalInd << endl ;
} // endfunction
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.);
+ fhAllPositionX = new TH1F("hAllPositionX", "#Delta X of any RP with primary photon",100, -2., 2.);
+ fhAllPositionZ = new TH1F("hAllPositionZ", "#Delta X of any RP with primary photon",100, -2., 2.);
+
// if(fhAllReg)
// delete fhAllReg ;
// if(fhPhotReg)
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 ;
-
- //========== Open galice root file
-
- if ( fRootFile == 0 ) {
- Text_t * name = new Text_t[80] ;
- cout << "AnalyzeOneEvent > Enter file root file name : " ;
- cin >> name ;
- Bool_t ok = OpenRootFile(name) ;
- if ( !ok )
- cout << " AliPHOSAnalyze > Error opening " << name << endl ;
- else {
- //========== 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 = fPHOS->GetGeometry();
- // fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
-
- } // else !ok
- } // if fRootFile
-
- if ( ok ) {
-
- //========== Create the Clusterizer
-
- fClu = new AliPHOSClusterizerv1() ;
- fClu->SetEmcEnergyThreshold(0.030) ;
- fClu->SetEmcClusteringThreshold(0.20) ;
- fClu->SetPpsdEnergyThreshold (0.0000002) ;
- fClu->SetPpsdClusteringThreshold(0.0000001) ;
- fClu->SetLocalMaxCut(0.03) ;
- fClu->SetCalibrationParameters(0., 0.00000001) ;
- cout << "AnalyzeOneEvent > using clusterizer " << fClu->GetName() << endl ;
- fClu->PrintParameters() ;
-
- //========== Creates the track segment maker
-
- 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 -> SetDebugReconstruction(kFALSE);
- fRec -> SetDebugReconstruction(kTRUE);
-
- //=========== Connect the various Tree's for evt
-
- if ( evt == -999 ) {
- cout << "AnalyzeOneEvent > Enter event number : " ;
- cin >> evt ;
- cout << evt << endl ;
- }
- fEvt = evt ;
- gAlice->GetEvent(evt);
-
- //=========== Get the Digit TTree
-
- gAlice->TreeD()->GetEvent(0) ;
-
- } // ok
-
- return ok ;
-}
-
-
-//____________________________________________________________________________
-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 ;
-
- Int_t module ;
- cout << "DisplayKineEvent > which module (1-5, -1: all) ? " ;
- cin >> module ; cout << module << endl ;
-
- Int_t testparticle ;
- cout << " 22 : PHOTON " << endl
- << " (-)11 : (POSITRON)ELECTRON " << endl
- << " (-)2112 : (ANTI)NEUTRON " << endl
- << " -999 : Everything else " << endl ;
- cout << "DisplayKineEvent > enter PDG particle code to display " ;
- cin >> testparticle ; cout << testparticle << endl ;
-
- Text_t histoname[80] ;
- 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, AliPHOSGeometry::Degre() ) ;
-
- Double_t theta, phi ;
- fGeom->EmcXtalCoverage(theta, phi, AliPHOSGeometry::Degre() ) ;
-
- Int_t tdim = (Int_t)( (tM - tm) / theta ) ;
- Int_t pdim = (Int_t)( (pM - pm) / phi ) ;
-
- tm -= theta ;
- tM += theta ;
- pm -= phi ;
- pM += phi ;
-
- TH2F * histoparticle = new TH2F("histoparticle", histoname,
- pdim, pm, pM, tdim, tm, tM) ;
- histoparticle->SetStats(kFALSE) ;
-
- // Get pointers to Alice Particle TClonesArray
-
- TParticle * particle;
- TClonesArray * particlearray = gAlice->Particles();
-
- Text_t canvasname[80];
- sprintf(canvasname,"Particles incident in PHOS/EMC module # %d",module) ;
- TCanvas * kinecanvas = new TCanvas("kinecanvas", canvasname, 650, 500) ;
-
- // get the KINE Tree
-
- TTree * kine = gAlice->TreeK() ;
- Stat_t nParticles = kine->GetEntries() ;
- cout << "DisplayKineEvent > events in kine " << nParticles << endl ;
-
- // loop over particles
-
- Double_t kRADDEG = 180. / TMath::Pi() ;
- Int_t index1 ;
- Int_t nparticlein = 0 ;
- for (index1 = 0 ; index1 < nParticles ; index1++){
- Int_t nparticle = particlearray->GetEntriesFast() ;
- Int_t index2 ;
- for( index2 = 0 ; index2 < nparticle ; index2++) {
- particle = (TParticle*)particlearray->UncheckedAt(index2) ;
- Int_t particletype = particle->GetPdgCode() ;
- if (testparticle == -999 || testparticle == particletype) {
- Double_t phi = particle->Phi() ;
- Double_t theta = particle->Theta() ;
- Int_t mod ;
- Double_t x, z ;
- fGeom->ImpactOnEmc(theta, phi, mod, z, x) ;
- if ( mod == module ) {
- nparticlein++ ;
- if (particle->Energy() > fClu->GetEmcClusteringThreshold() )
- histoparticle->Fill(phi*kRADDEG, theta*kRADDEG, particle->Energy() ) ;
- }
- }
- }
- }
- kinecanvas->Draw() ;
- histoparticle->Draw("color") ;
- TPaveText * pavetext = new TPaveText(294, 100, 300, 101);
- Text_t text[40] ;
- sprintf(text, "Particles: %d ", nparticlein) ;
- pavetext->AddText(text) ;
- pavetext->Draw() ;
- kinecanvas->Update();
-
-}
-//____________________________________________________________________________
-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 << "DisplayRecParticles > Analyze an event first ... (y/n) " ;
- Text_t answer[1] ;
- cin >> answer ; cout << answer ;
-// if ( answer == "y" )
-// AnalyzeOneEvent() ;
- }
- if (fEvt != -999) {
-
- Int_t module ;
- 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, AliPHOSGeometry::Degre() ) ;
- Double_t theta, phi ;
- fGeom->EmcXtalCoverage(theta, phi, AliPHOSGeometry::Degre() ) ;
- Int_t tdim = (Int_t)( (tM - tm) / theta ) ;
- Int_t pdim = (Int_t)( (pM - pm) / phi ) ;
- tm -= theta ;
- tM += theta ;
- pm -= phi ;
- TH2F * histoRparticle = new TH2F("histoRparticle", histoname,
- pdim, pm, pM, tdim, tm, tM) ;
- histoRparticle->SetStats(kFALSE) ;
- Text_t canvasname[80] ;
- sprintf(canvasname, "Reconstructed particles in PHOSmodule # %d", module) ;
- TCanvas * rparticlecanvas = new TCanvas("RparticleCanvas", canvasname, 650, 500) ;
- AliPHOSRecParticle::RecParticlesList * rpl = *fPHOS->RecParticles() ;
- Int_t nRecParticles = rpl->GetEntries() ;
- Int_t nRecParticlesInModule = 0 ;
- TIter nextRecPart(rpl) ;
- AliPHOSRecParticle * rp ;
- cout << "DisplayRecParticles > " << nRecParticles << " reconstructed particles " << endl ;
- Double_t kRADDEG = 180. / TMath::Pi() ;
- while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
- AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ;
- 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 ;
- Double_t energy = rp->Energy() ;
- histoRparticle->Fill(phi, theta, energy) ;
- }
- }
- 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);
- pavetext->AddText(text) ;
- pavetext->Draw() ;
- rparticlecanvas->Update() ;
- }
-}
-
-//____________________________________________________________________________
-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 (fEvt != -999) {
-
- Int_t module ;
- cout << "DisplayRecPoints > which module (1-5, -1: all) ? " ;
- cin >> module ; cout << module << endl ;
-
- Text_t canvasname[80];
- sprintf(canvasname,"Digits in PHOS/EMC module # %d",module) ;
- TCanvas * modulecanvas = new TCanvas("module", canvasname, 650, 500) ;
- modulecanvas->Draw() ;
-
- //=========== Creating 2d-histogram of the PHOS module
- // a little bit junkie but is used to test Geom functinalities
-
- Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
-
- fGeom->EmcModuleCoverage(module, tm, tM, pm, pM);
- // convert angles into coordinates local to the EMC module of interest
-
- Int_t emcModuleNumber ;
- Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
- Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
- fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
- fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
- Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;
- Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
- Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ;
- Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ;
- Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ;
- Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;
- Text_t histoname[80];
- sprintf(histoname,"Event %d: Digits and RecPoints in module %d", fEvt, module) ;
- TH2F * hModule = new TH2F("HistoReconstructed", histoname,
- xdim, xmin, xMax, zdim, zmin, zMax) ;
- hModule->SetMaximum(2.0);
- hModule->SetMinimum(0.0);
- hModule->SetStats(kFALSE);
-
- TIter next(fPHOS->Digits()) ;
- Float_t energy, y, z;
- 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 && relid[1] == 0)
- {
- energy = fClu->Calibrate(digit->GetAmp()) ;
- 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 "
- << module << " " << nDigits << " digits with total energy " << etot << endl ;
- hModule->Draw("col2") ;
-
- //=========== Cluster in module
-
- // TClonesArray * emcRP = fPHOS->EmcClusters() ;
- TObjArray * emcRP = *(fPHOS->EmcRecPoints()) ;
-
- etot = 0.;
- Int_t totalnClusters = 0 ;
- Int_t nClusters = 0 ;
- TIter nextemc(emcRP) ;
- AliPHOSEmcRecPoint * emc ;
- while((emc = (AliPHOSEmcRecPoint *)nextemc()))
- {
- // Int_t numberofprimaries ;
- // Int_t * primariesarray = new Int_t[10] ;
- // emc->GetPrimaries(numberofprimaries, primariesarray) ;
- totalnClusters++ ;
- if ( emc->GetPHOSMod() == module )
- {
- nClusters++ ;
- energy = emc->GetTotalEnergy() ;
- etot+= energy ;
- emc->Draw("M") ;
- }
- }
- cout << "DrawRecPoints > Found " << totalnClusters << " EMC Clusters in PHOS" << endl ;
- cout << "DrawRecPoints > Found in module " << module << " " << nClusters << " EMC Clusters " << endl ;
- cout << "DrawRecPoints > total energy " << etot << endl ;
-
- TPaveText * pavetext = new TPaveText(22, 80, 83, 90);
- Text_t text[40] ;
- sprintf(text, "digits: %d; clusters: %d", nDigits, nClusters) ;
- pavetext->AddText(text) ;
- pavetext->Draw() ;
- modulecanvas->Update();
-
- //=========== Cluster in module PPSD Down
-
- // TClonesArray * ppsdRP = fPHOS->PpsdClusters() ;
- TObjArray * ppsdRP = *(fPHOS->PpsdRecPoints() );
-
- etot = 0.;
- TIter nextPpsd(ppsdRP) ;
- AliPHOSPpsdRecPoint * ppsd ;
- while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
- {
- totalnClusters++ ;
- if ( ppsd->GetPHOSMod() == module )
- {
- nClusters++ ;
- energy = ppsd->GetEnergy() ;
- etot+=energy ;
- if (!ppsd->GetUp()) ppsd->Draw("P") ;
- }
- }
- cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Down Clusters in PHOS" << endl ;
- cout << "DrawRecPoints > Found in module " << module << " " << nClusters << " Ppsd Down Clusters " << endl ;
- cout << "DrawRecPoints > total energy " << etot << endl ;
-
- //=========== Cluster in module PPSD Up
-
- ppsdRP = *(fPHOS->PpsdRecPoints()) ;
-
- etot = 0.;
- TIter nextPpsdUp(ppsdRP) ;
- while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsdUp()))
- {
- totalnClusters++ ;
- if ( ppsd->GetPHOSMod() == module )
- {
- nClusters++ ;
- energy = ppsd->GetEnergy() ;
- etot+=energy ;
- if (ppsd->GetUp()) ppsd->Draw("P") ;
- }
- }
- cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Up Clusters in PHOS" << endl ;
- cout << "DrawRecPoints > Found in module " << module << " " << nClusters << " Ppsd Up Clusters " << endl ;
- cout << "DrawRecPoints > total energy " << etot << endl ;
-
- } // if !-999
-}
-
-//____________________________________________________________________________
-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 (fEvt != -999) {
-
- Int_t module ;
- cout << "DisplayTrackSegments > which module (1-5, -1: all) ? " ;
- cin >> module ; cout << module << endl ;
- //=========== Creating 2d-histogram of the PHOS module
- // a little bit junkie but is used to test Geom functinalities
-
- Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
-
- fGeom->EmcModuleCoverage(module, tm, tM, pm, pM);
- // convert angles into coordinates local to the EMC module of interest
-
- Int_t emcModuleNumber ;
- Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
- Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
- fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
- fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
- Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;
- Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
- Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ;
- Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ;
- Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ;
- Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;
- Text_t histoname[80];
- sprintf(histoname,"Event %d: Track Segments in module %d", fEvt, module) ;
- TH2F * histotrack = new TH2F("histotrack", histoname,
- xdim, xmin, xMax, zdim, zmin, zMax) ;
- histotrack->SetStats(kFALSE);
- Text_t canvasname[80];
- sprintf(canvasname,"Track segments in PHOS/EMC-PPSD module # %d", module) ;
- TCanvas * trackcanvas = new TCanvas("TrackSegmentCanvas", canvasname, 650, 500) ;
- histotrack->Draw() ;
-
- AliPHOSTrackSegment::TrackSegmentsList * trsegl = *(fPHOS->TrackSegments()) ;
- AliPHOSTrackSegment * trseg ;
-
- Int_t nTrackSegments = trsegl->GetEntries() ;
- Int_t index ;
- Float_t etot = 0 ;
- Int_t nTrackSegmentsInModule = 0 ;
- for(index = 0; index < nTrackSegments ; index++){
- trseg = (AliPHOSTrackSegment * )trsegl->At(index) ;
- etot+= trseg->GetEnergy() ;
- if ( trseg->GetPHOSMod() == module ) {
- nTrackSegmentsInModule++ ;
- trseg->Draw("P");
- }
- }
- Text_t text[80] ;
- sprintf(text, "track segments: %d", nTrackSegmentsInModule) ;
- TPaveText * pavetext = new TPaveText(22, 80, 83, 90);
- pavetext->AddText(text) ;
- pavetext->Draw() ;
- trackcanvas->Update() ;
- cout << "DisplayTrackSegments > Found " << trsegl->GetEntries() << " Track segments with total energy "<< etot << endl ;
-
- }
-}
//____________________________________________________________________________
Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
{
fRootFile = new TFile(name, "update") ;
return fRootFile->IsOpen() ;
}
+
//____________________________________________________________________________
void AliPHOSAnalyze::SaveHistograms()
{
fhPPSDEnergy->Write() ;
if(fhAllPosition)
fhAllPosition->Write() ;
+ if(fhAllPositionX)
+ fhAllPositionX->Write() ;
+ if(fhAllPositionZ)
+ fhAllPositionZ->Write() ;
if(fhPhotPosition)
fhPhotPosition->Write() ;
if(fhEMPosition)
fhEMEnergy = 0 ; // Spectrum of detected electrons with electron primary
fhPPSDEnergy = 0 ;
fhAllPosition = 0 ;
+ fhAllPositionX = 0 ;
+ fhAllPositionZ = 0 ;
fhPhotPosition = 0 ;
fhEMPosition = 0 ;
fhPPSDPosition = 0 ;
fhPhotChHa = 0 ;
fhPhotGaHa = 0 ;
-
}
AliPHOSAnalyze(const AliPHOSAnalyze & ana) ; // cpy ctor
virtual ~AliPHOSAnalyze() ; // dtor
- void ActivePPSD(Int_t Nevents) ;
- void AnalyzeManyEvents(Int_t Nevtents = 100, Int_t Module=0) ; // analyzes many events ;
- void AnalyzeResolutions(Int_t Nevtents) ; // analyzes Energy and Position resolutions ;
+ void DrawRecon(Int_t Nevent= 0,Int_t Nmod = 1) ; // draws positions of entering and
+ //reconstructed points in PHOS plain
+ void AnalyzeResolutions (Int_t Nevents) ; // analyzes Energy and Position resolutions ;
+ void ReadAndPrintEMC(Int_t EvFirst=0, Int_t EvLast=0); // Read & print generated and reconstructed hits in EMC
void ReadAndPrintCPV(Int_t EvFirst=0, Int_t EvLast=0); // Read & print generated and reconstructed hits in CPV
void AnalyzeCPV(Int_t Nevents); // analyzes various CPV characteristics
+ void AnalyzeEMC(Int_t Nevents); // analyzes EMC resolution
void InvariantMass(Int_t Nevents = 100) ;
void Reconstruct(Int_t Nevtents = 100,Int_t FirstEvent = 0) ;
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
- void DisplayRecPoints() ; // displays RecPoints in module coordinate
- void DisplayTrackSegments() ; // displays TrackSegments in module coordinate
Bool_t OpenRootFile(Text_t * name) ; // opens the root file
void SaveHistograms() ; // Save histograms in a root file
void ResetHistograms() ; //
TH1F * fhConvertorCluster ; // Histo of Cluster energies in Convertor
TH2F * fhConvertorEmc ; // 2d Convertor versus Emc energies
- TH2F * fhAllEnergy ; // Spectrum of detected photons with photon primary
- TH2F * fhPhotEnergy ; // Total spectrum of detected photons
+ 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 * 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
+ 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 * fhAllPositionX ; // X-Position Resolution of photons with photon primary
+ TH1F * fhAllPositionZ ; // Z-Position Resolution of photons with photon primary
TH1F * fhPhotonPositionY ; // Y distribution of detected photons
TH1F * fhElectronPositionY ; // Y distribution of detected electrons
virtual Float_t GetPhiDisplacement(void) = 0 ;
virtual Float_t GetPPSDModuleSize(Int_t index) = 0 ;
virtual Float_t GetZDisplacement(void) = 0 ;
+ virtual void SetLeadConverterThickness(Float_t x) = 0 ;
// Return CPV geometrical parameters
virtual Int_t GetNumberOfPadsZ(void) { AssertPPSD("GetNumberOfPadsZ"); return 0; }
virtual Int_t GetNumberOfModulesPhi(void) { AssertPPSD("GetNumberOfModulesPhi"); return 0; }
virtual Int_t GetNumberOfModulesZ(void) { AssertPPSD("GetNumberOfModulesZ"); return 0; }
+ virtual void SetLeadConverterThickness(Float_t x) { AssertPPSD("SetLeadConverterThickness"); }
private:
/*
$Log$
+ Revision 1.4 2000/11/20 09:50:03 schutz
+ AliPHOSCPVHit inherits from AliHit
+
Revision 1.3 2000/11/13 11:45:36 schutz
DEC and HP compilers sufisfied
//______________________________________________________________________________
-void AliPHOSCPVModule::MakeBranch(Int_t i)
+void AliPHOSCPVModule::MakeBranch(char *title,Int_t i)
{
//
- // Create a new branch for a CPV module #i in the current Root Tree
+ // Create a new branch for a EMC or CPV module #i in the current Root Tree
//
char branchname[10];
- sprintf(branchname,"CPV%d",i);
+ sprintf(branchname,"%s%d",title,i);
gAlice->TreeH()->Branch(branchname,&fHits, 1000);
}
//_____________________________________________________________________________
-void AliPHOSCPVModule::SetTreeAddress(Int_t i)
+void AliPHOSCPVModule::SetTreeAddress(char *title, Int_t i)
{
//
- // Set branch address for the Hits Tree for a CPV module #i
+ // Set branch address for the Hits Tree for a CPV or EMC module #i
//
TBranch *branch;
char branchname[20];
TTree *treeH = gAlice->TreeH();
if (treeH){
- sprintf(branchname,"CPV%d",i);
+ sprintf(branchname,"%s%d",title,i);
branch = treeH->GetBranch(branchname);
if (branch) branch->SetAddress(&fHits);
}
void Clear(Option_t *opt="");
void Print(Option_t *opt="");
void AddHit(Int_t shunt, Int_t track, TLorentzVector p, Float_t *xy, Int_t ipart);
- void MakeBranch(Int_t i);
- void SetTreeAddress(Int_t i);
+ void MakeBranch (char *title, Int_t i);
+ void SetTreeAddress(char *title, Int_t i);
TClonesArray *Hits (void) {return fHits;}
virtual void SetLogWeightCutCPV(Float_t w) = 0 ;
virtual void SetPpsdClusteringThreshold(Float_t cluth) = 0 ;
virtual void SetPpsdEnergyThreshold(Float_t enth) = 0 ;
+ virtual void SetCpvClusteringThreshold(Float_t cluth) = 0 ;
+ virtual void SetCpvEnergyThreshold(Float_t enth) = 0 ;
ClassDef(AliPHOSClusterizer,1) // Clusterization algorithm class (abstract base class)
//*-- Author: Yves Schutz (SUBATECH)
//////////////////////////////////////////////////////////////////////////////
+#include <assert.h>
+
// --- ROOT system ---
#include "TMath.h"
fEmcEnergyThreshold = 0.01;
fPpsdClusteringThreshold = 0.0;
fPpsdEnergyThreshold = 0.1;
+ fCpvClusteringThreshold = 0.0;
+ fCpvEnergyThreshold = 0.1;
fW0 = 4.5 ;
fLocMaxCut = 0.06 ;
fW0CPV = 4.5 ;
}
//Do NOT clusterize upper PPSD
- if(fGeom->GetNumberOfCPVLayers() > 1 &&
- relid1[1] > 0 &&
+ if( IsInPpsd(d1) && IsInPpsd(d2) &&
+ relid1[1] > 0 &&
relid1[1] < fGeom->GetNumberOfPadsPhi()*fGeom->GetNumberOfPadsPhi() ) rv = 2 ;
return rv ;
Int_t relid[4] ;
fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
- if ( relid[1] == 0 )
- rv = kTRUE;
+ if ( relid[1] == 0 ) rv = kTRUE;
+
+ return rv ;
+}
+
+//____________________________________________________________________________
+Bool_t AliPHOSClusterizerv1::IsInPpsd(AliPHOSDigit * digit)
+{
+ // Tells if (true) or not (false) the digit is in a PHOS-EMC module
+
+ Bool_t rv = kFALSE ;
+
+ Int_t relid[4] ;
+ fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
+
+ if ( relid[1] > 0 && relid[0] > fGeom->GetNCPVModules() ) rv = kTRUE;
+
+ return rv ;
+}
+
+//____________________________________________________________________________
+Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit)
+{
+ // Tells if (true) or not (false) the digit is in a PHOS-EMC module
+
+ Bool_t rv = kFALSE ;
+
+ Int_t relid[4] ;
+ fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
+
+ if ( relid[1] > 0 && relid[0] <= fGeom->GetNCPVModules() ) rv = kTRUE;
return rv ;
}
AliPHOSDigit ** clusterdigitslist = new AliPHOSDigit*[dl->GetEntries()] ;
Int_t index ;
- if (( ( IsInEmc(digit) ) && ( Calibrate(digit->GetAmp() ) > fEmcClusteringThreshold ) ) ||
- ( ( !IsInEmc(digit) ) && ( Calibrate(digit->GetAmp() ) > fPpsdClusteringThreshold ) ) ) {
+ if (( ( IsInEmc (digit) ) && ( Calibrate(digit->GetAmp() ) > fEmcClusteringThreshold ) ) ||
+ ( ( IsInPpsd(digit) ) && ( Calibrate(digit->GetAmp() ) > fPpsdClusteringThreshold ) ) ||
+ ( ( IsInCpv (digit) ) && ( Calibrate(digit->GetAmp() ) > fCpvClusteringThreshold ) ) ) {
Int_t iDigitInCluster = 0 ;
// start a new PPSD cluster
if(fNumberOfPpsdClusters >= ppsdl->GetSize()) ppsdl->Expand(2*fNumberOfPpsdClusters+1);
- if (strcmp(fGeom->GetName(),"GPS2") == 0) {
+ if (IsInPpsd(digit)) {
(*ppsdl)[fNumberOfPpsdClusters] = new AliPHOSPpsdRecPoint() ;
clu = (AliPHOSPpsdRecPoint *) ppsdl->At(fNumberOfPpsdClusters) ;
}
- else if (strcmp(fGeom->GetName(),"IHEP") == 0) {
+ else if (IsInCpv(digit) ) {
(*ppsdl)[fNumberOfPpsdClusters] = new AliPHOSCpvRecPoint(fW0CPV, fLocMaxCutCPV) ;
clu = (AliPHOSCpvRecPoint *) ppsdl->At(fNumberOfPpsdClusters) ;
}
+ else {
+ cout << "AliPHOSClusterizerv1::MakeClusters: unknown configuration " << fGeom->GetName() << endl;
+ assert(0==1);
+ }
fNumberOfPpsdClusters++ ;
clu->AddDigit(*digit, Calibrate(digit->GetAmp()) ) ;
clusterdigitslist[iDigitInCluster] = digit ;
<< " EMC Clustering threshold = " << fEmcClusteringThreshold << endl
<< " EMC Energy threshold = " << fEmcEnergyThreshold << endl
<< " PPSD Clustering threshold = " << fPpsdClusteringThreshold << endl
- << " PPSD Energy threshold = " << fPpsdEnergyThreshold << endl ;
+ << " PPSD Energy threshold = " << fPpsdEnergyThreshold << endl
+ << " CPV Clustering threshold = " << fCpvClusteringThreshold << endl
+ << " CPV Energy threshold = " << fCpvEnergyThreshold << endl ;
}
virtual Float_t GetLogWeightCutCPV() { return fW0CPV;}
virtual Float_t GetPpsdClusteringThreshold() { return fPpsdClusteringThreshold; }
virtual Float_t GetPpsdEnergyThreshold() { return fPpsdEnergyThreshold; }
+ virtual Float_t GetCpvClusteringThreshold() { return fCpvClusteringThreshold; }
+ virtual Float_t GetCpvEnergyThreshold() { return fCpvEnergyThreshold; }
- virtual Bool_t IsInEmc(AliPHOSDigit * digit) ; // Tells if id digit is in EMC
+ virtual Bool_t IsInEmc (AliPHOSDigit * digit) ; // Tells if id digit is in EMC
+ virtual Bool_t IsInPpsd(AliPHOSDigit * digit) ; // Tells if id digit is in PPSD
+ virtual Bool_t IsInCpv (AliPHOSDigit * digit) ; // Tells if id digit is in CPV
virtual void MakeClusters(const DigitsList * dl,
AliPHOSRecPoint::RecPointsList * emcl,
AliPHOSRecPoint::RecPointsList * ppsdl) ; // does the job
virtual void SetLogWeightCutCPV(Float_t w) { fW0CPV = w ; }
virtual void SetPpsdClusteringThreshold(Float_t cluth) { fPpsdClusteringThreshold = cluth ; }
virtual void SetPpsdEnergyThreshold(Float_t enth) { fPpsdEnergyThreshold = enth ; }
+ virtual void SetCpvClusteringThreshold(Float_t cluth) { fCpvClusteringThreshold = cluth ; }
+ virtual void SetCpvEnergyThreshold(Float_t enth) { fCpvEnergyThreshold = enth ; }
private:
Float_t fA ; // offset of the energy calibration
Float_t fB ; // gain of the energy calibration
AliPHOSGeometry * fGeom ; // pointer to geometry
+ Int_t fNumberOfEmcClusters ; // number of EMC clusters found
Float_t fEmcClusteringThreshold ; // minimum energy to include a EMC digit in a cluster
Float_t fEmcEnergyThreshold ; // minimum energy of EMC digit to be considered
- Int_t fNumberOfEmcClusters ; // number of EMC clusters found
Int_t fNumberOfPpsdClusters ; // number of PPSD clusters found
Float_t fPpsdClusteringThreshold ; // minimum energy to include a PPSD digit in a cluster
Float_t fPpsdEnergyThreshold ; // minimum energy of PPSD digit to be considered
+ Int_t fNumberOfCpvClusters ; // number of CPV clusters found
+ Float_t fCpvClusteringThreshold ; // minimum energy to include a CPV digit in a cluster
+ Float_t fCpvEnergyThreshold ; // minimum energy of CPV digit to be considered
Float_t fLocMaxCut ; // minimum energy difference to distinguish local maxima in a cluster
Float_t fW0 ; // logarithmic weight for the cluster center of gravity calculation
-
Float_t fLocMaxCutCPV ; // minimum energy difference to distinguish local maxima in a CPV cluster
Float_t fW0CPV ; // logarithmic weight for the CPV cluster center of gravity calculation
- ClassDef(AliPHOSClusterizerv1,1) // Clusterizer implementation version 1
+ ClassDef(AliPHOSClusterizerv1,1) // Clusterizer implementation version 1
};
void GetClusterLengths(Int_t &lengX, Int_t &lengZ); // cluster lengths along x and z
Bool_t IsCPV(void) {
// true if the recpoint is in CPV
- return kTRUE ; }
+ return (fPHOSMod <= fGeom->GetNCPVModules()) ; }
Bool_t IsSortable() const {
// says that emcrecpoints are sortable objects
return kTRUE ; }
fNZ = 64 ;
fXtlSize[0] = 2.2 ;
- fXtlSize[1] = 18.0 ;
+ fXtlSize[1] = 22.0 ;
fXtlSize[2] = 2.2 ;
// all these numbers coming next are subject to changes
- fOuterBoxThickness[0] = 2.8 ;
+ fOuterBoxThickness[0] = 2.5 ;
fOuterBoxThickness[1] = 5.0 ;
fOuterBoxThickness[2] = 5.0 ;
fTextolitBoxThickness[1] = 0.0 ;
fTextolitBoxThickness[2] = 3.0 ;
- fAirThickness[0] = 1.56 ;
+ fAirThickness[0] = 0.4 ;
fAirThickness[1] = 20.5175 ;
fAirThickness[2] = 2.48 ;
if (fRotMatrixArray) fRotMatrixArray->Delete() ;
if (fRotMatrixArray) delete fRotMatrixArray ;
if (fPHOSAngle ) delete fPHOSAngle ;
-// if (fGeometryEMCA ) detete fGeometryEMCA;
-// if (fGeometryCPV ) detete fGeometryCPV ;
}
//____________________________________________________________________________
if ( ((strcmp( fName, "default" )) == 0) ||
((strcmp( fName, "GPS2" )) == 0) ||
- ((strcmp( fName, "IHEP" )) == 0) ) {
+ ((strcmp( fName, "IHEP" )) == 0) ||
+ ((strcmp( fName, "MIXT" )) == 0) ) {
fgInit = kTRUE ;
- fGeometryEMCA = new AliPHOSEMCAGeometry();
- if ( ((strcmp( fName, "GPS2" )) == 0) ) fGeometryCPV = new AliPHOSPPSDGeometry();
- if ( ((strcmp( fName, "IHEP" )) == 0) ) fGeometryCPV = new AliPHOSCPVGeometry ();
- fGeometrySUPP = new AliPHOSSupportGeometry();
- fNModules = 5;
+
+ fNModules = 5;
+ fNPPSDModules = 0;
+ fAngle = 20;
+
+ fGeometryEMCA = new AliPHOSEMCAGeometry();
+ if ( ((strcmp( fName, "GPS2" )) == 0) ) {
+ fGeometryPPSD = new AliPHOSPPSDGeometry();
+ fGeometryCPV = 0;
+ fNPPSDModules = fNModules;
+ }
+ else if ( ((strcmp( fName, "IHEP" )) == 0) ) {
+ fGeometryCPV = new AliPHOSCPVGeometry ();
+ fGeometryPPSD = 0;
+ fNPPSDModules = 0;
+ }
+ else if ( ((strcmp( fName, "MIXT" )) == 0) ) {
+ fGeometryCPV = new AliPHOSCPVGeometry ();
+ fGeometryPPSD = new AliPHOSPPSDGeometry();
+ fNPPSDModules = 1;
+ }
+ fGeometrySUPP = new AliPHOSSupportGeometry();
+
fPHOSAngle = new Float_t[fNModules] ;
Int_t index ;
for ( index = 0; index < fNModules; index++ )
this->SetPHOSAngles() ;
fRotMatrixArray = new TObjArray(fNModules) ;
}
- else {
- fgInit = kFALSE ;
- cout << "PHOS Geometry setup: option not defined " << fName << endl ;
- }
+ else {
+ fgInit = kFALSE ;
+ cout << "PHOS Geometry setup: option not defined " << fName << endl ;
+ }
}
//____________________________________________________________________________
// Calculates the position in ALICE of the PHOS modules
Double_t const kRADDEG = 180.0 / kPI ;
- Float_t pphi = TMath::ATan( GetOuterBoxSize(0) / ( 2.0 * GetIPtoOuterCoverDistance() ) ) ;
+ Float_t pphi = 2 * TMath::ATan( GetOuterBoxSize(0) / ( 2.0 * GetIPtoOuterCoverDistance() ) ) ;
pphi *= kRADDEG ;
+ if (pphi > fAngle) cout << "AliPHOSGeometry: PHOS modules overlap!\n";
+ pphi = fAngle;
for( Int_t i = 1; i <= fNModules ; i++ ) {
- Float_t angle = pphi * 2 * ( i - fNModules / 2.0 - 0.5 ) ;
+ Float_t angle = pphi * ( i - fNModules / 2.0 - 0.5 ) ;
fPHOSAngle[i-1] = - angle ;
}
}
}
else if ( strcmp(fName,"IHEP") == 0 ) {
id -= GetNPhi() * GetNZ() * GetNModules() ;
- relid[0] = (Int_t) TMath::Ceil( id / ( GetNumberOfPadsPhi() * GetNumberOfPadsZ() ) ) ;
+ Float_t nCPV = GetNumberOfCPVPadsPhi() * GetNumberOfCPVPadsZ() ;
+ relid[0] = (Int_t) TMath::Ceil( id / nCPV ) ;
relid[1] = 1 ;
- id -= ( relid[0] - 1 ) * GetNumberOfPadsPhi() * GetNumberOfPadsZ() ;
- relid[2] = (Int_t) TMath::Ceil( id / GetNumberOfPadsZ() ) ;
- relid[3] = (Int_t) ( id - ( relid[2] - 1 ) * GetNumberOfPadsZ() ) ;
+ id -= ( relid[0] - 1 ) * nCPV ;
+ relid[2] = (Int_t) TMath::Ceil( id / GetNumberOfCPVPadsZ() ) ;
+ relid[3] = (Int_t) ( id - ( relid[2] - 1 ) * GetNumberOfCPVPadsZ() ) ;
+ }
+ else if ( strcmp(fName,"MIXT") == 0 ) {
+ id -= GetNPhi() * GetNZ() * GetNModules() ;
+ Float_t nPPSD = 2 * GetNumberOfModulesPhi() * GetNumberOfModulesZ() * GetNumberOfPadsPhi() * GetNumberOfPadsZ() ;
+ Float_t nCPV = GetNumberOfCPVPadsPhi() * GetNumberOfCPVPadsZ() ;
+ if (id <= nCPV*GetNCPVModules()) { // this pad belons to CPV
+ relid[0] = (Int_t) TMath::Ceil( id / nCPV ) ;
+ relid[1] = 1 ;
+ id -= ( relid[0] - 1 ) * nCPV ;
+ relid[2] = (Int_t) TMath::Ceil( id / GetNumberOfCPVPadsZ() ) ;
+ relid[3] = (Int_t) ( id - ( relid[2] - 1 ) * GetNumberOfCPVPadsZ() ) ;
+ }
+ else { // this pad belons to PPSD
+ id -= nCPV*GetNCPVModules();
+ relid[0] = (Int_t)TMath::Ceil( id / nPPSD );
+ id -= ( relid[0] - 1 ) * nPPSD ;
+ relid[0] += GetNCPVModules();
+ relid[1] = (Int_t)TMath::Ceil( id / ( GetNumberOfPadsPhi() * GetNumberOfPadsZ() ) ) ;
+ id -= ( relid[1] - 1 ) * GetNumberOfPadsPhi() * GetNumberOfPadsZ() ;
+ relid[2] = (Int_t)TMath::Ceil( id / GetNumberOfPadsPhi() ) ;
+ relid[3] = (Int_t) ( id - ( relid[2] - 1 ) * GetNumberOfPadsPhi() ) ;
+ }
}
}
else { // its a PW04 crystal
Bool_t AliPHOSGeometry::RelToAbsNumbering(const Int_t * relid, Int_t & AbsId)
{
// Converts the relative numbering into the absolute numbering
- // AbsId = 1 to fNModules * fNPhi * fNZ -> PbWO4
- // AbsId = N(total PHOS crystals) +
- // 1 to fNModules * 2 * (fNumberOfModulesPhi * fNumberOfModulesZ) * fNumberOfPadsPhi * fNumberOfPadsZ -> PPSD
- // AbsId = N(total PHOS crystals) +
- // 1:fNModules * fNumberOfCPVPadsPhi * fNumberOfCPVPadsZ -> CPV
+ // EMCA crystals:
+ // AbsId = from 1 to fNModules * fNPhi * fNZ
+ // PPSD gas cell:
+ // AbsId = from N(total EMCA crystals) + 1
+ // to NCPVModules * fNumberOfCPVPadsPhi * fNumberOfCPVPadsZ +
+ // fNModules * 2 * (fNumberOfModulesPhi * fNumberOfModulesZ) * fNumberOfPadsPhi * fNumberOfPadsZ
+ // CPV pad:
+ // AbsId = from N(total PHOS crystals) + 1
+ // to NCPVModules * fNumberOfCPVPadsPhi * fNumberOfCPVPadsZ
Bool_t rv = kTRUE ;
- if ( relid[1] > 0 ) { // it is a PPSD pad
- AbsId = GetNPhi() * GetNZ() * GetNModules() // the offset to separate EMCA crystals from PPSD pads
- + ( relid[0] - 1 ) * GetNumberOfModulesPhi() * GetNumberOfModulesZ() // the pads offset of PHOS modules
+ if ( relid[1] > 0 && strcmp(fName,"GPS2")==0) { // it is a PPSD pad
+ AbsId = GetNPhi() * GetNZ() * GetNModules() // the offset to separate EMCA crystals from PPSD pads
+ + ( relid[0] - 1 ) * GetNumberOfModulesPhi() * GetNumberOfModulesZ() // the pads offset of PPSD modules
* GetNumberOfPadsPhi() * GetNumberOfPadsZ() * 2
- + ( relid[1] - 1 ) * GetNumberOfPadsPhi() * GetNumberOfPadsZ() // the pads offset of PPSD modules
- + ( relid[2] - 1 ) * GetNumberOfPadsPhi() // the pads offset of a PPSD row
- + relid[3] ; // the column number
+ + ( relid[1] - 1 ) * GetNumberOfPadsPhi() * GetNumberOfPadsZ() // the pads offset of PPSD modules
+ + ( relid[2] - 1 ) * GetNumberOfPadsPhi() // the pads offset of a PPSD row
+ + relid[3] ; // the column number
+ }
+
+ else if ( relid[1] > 0 && strcmp(fName,"MIXT")==0) { // it is a PPSD pad
+ AbsId = GetNPhi() * GetNZ() * GetNModules() // the offset to separate EMCA crystals from PPSD pads
+ + GetNCPVModules() * GetNumberOfCPVPadsPhi() * GetNumberOfCPVPadsZ() // the pads offset of CPV modules if any
+ + ( relid[0] - 1 - GetNCPVModules())
+ * GetNumberOfModulesPhi() * GetNumberOfModulesZ() // the pads offset of PPSD modules
+ * GetNumberOfPadsPhi() * GetNumberOfPadsZ() * 2
+ + ( relid[1] - 1 ) * GetNumberOfPadsPhi() * GetNumberOfPadsZ() // the pads offset of PPSD modules
+ + ( relid[2] - 1 ) * GetNumberOfPadsPhi() // the pads offset of a PPSD row
+ + relid[3] ; // the column number
}
else if ( relid[1] == 0 ) { // it is a Phos crystal
else if ( relid[1] == -1 ) { // it is a CPV pad
AbsId = GetNPhi() * GetNZ() * GetNModules() // the offset to separate EMCA crystals from CPV pads
- + ( relid[0] - 1 ) * GetNumberOfPadsPhi() * GetNumberOfPadsZ() // the pads offset of PHOS modules
- + ( relid[2] - 1 ) * GetNumberOfPadsZ() // the pads offset of a CPV row
+ + ( relid[0] - 1 ) * GetNumberOfCPVPadsPhi() * GetNumberOfCPVPadsZ() // the pads offset of PHOS modules
+ + ( relid[2] - 1 ) * GetNumberOfCPVPadsZ() // the pads offset of a CPV row
+ relid[3] ; // the column number
}
void AliPHOSGeometry::RelPosInAlice(const Int_t id, TVector3 & pos )
{
// Converts the absolute numbering into the global ALICE coordinates
+ // It works only for the GPS2 geometry
- if (id > 0) {
-
- Int_t relid[4] ;
-
- AbsToRelNumbering(id , relid) ;
-
- Int_t phosmodule = relid[0] ;
-
- Float_t y0 = 0 ;
-
- if ( relid[1] == 0 ) // it is a PbW04 crystal
- { y0 = -(GetIPtoOuterCoverDistance() + GetUpperPlateThickness()
- + GetSecondUpperPlateThickness() + GetUpperCoolingPlateThickness()) ;
- }
- if ( relid[1] > 0 ) { // its a PPSD pad
- if ( relid[1] > GetNumberOfModulesPhi() * GetNumberOfModulesZ() ) // its an bottom module
- {
- y0 = -( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() / 2.0) ;
- }
- else // its an upper module
- y0 = -( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() - GetLeadToMicro2Gap()
- - GetLeadConverterThickness() - GetMicro1ToLeadGap() - GetMicromegas1Thickness() / 2.0) ;
- }
-
- Float_t x, z ;
- RelPosInModule(relid, x, z) ;
-
- pos.SetX(x) ;
- pos.SetZ(z) ;
- pos.SetY( TMath::Sqrt(x*x + z*z + y0*y0) ) ;
-
-
-
- Float_t phi = GetPHOSAngle( phosmodule) ;
- Double_t const kRADDEG = 180.0 / kPI ;
- Float_t rphi = phi / kRADDEG ;
-
- TRotation rot ;
- rot.RotateZ(-rphi) ; // a rotation around Z by angle
-
- TRotation dummy = rot.Invert() ; // to transform from original frame to rotate frame
-
- pos.Transform(rot) ; // rotate the baby
+ if (id > 0 && strcmp(fName,"GPS2")==0) {
+
+ Int_t relid[4] ;
+
+ AbsToRelNumbering(id , relid) ;
+
+ Int_t phosmodule = relid[0] ;
+
+ Float_t y0 = 0 ;
+
+ if ( relid[1] == 0 ) { // it is a PbW04 crystal
+ y0 = -(GetIPtoOuterCoverDistance() + GetUpperPlateThickness()
+ + GetSecondUpperPlateThickness() + GetUpperCoolingPlateThickness()) ;
+ }
+ if ( relid[1] > 0 ) { // its a PPSD pad
+ if ( relid[1] > GetNumberOfModulesPhi() * GetNumberOfModulesZ() ) { // its an bottom module
+ y0 = -( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() / 2.0) ;
+ }
+ else // its an upper module
+ y0 = -( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() - GetLeadToMicro2Gap()
+ - GetLeadConverterThickness() - GetMicro1ToLeadGap() - GetMicromegas1Thickness() / 2.0) ;
+ }
+
+ Float_t x, z ;
+ RelPosInModule(relid, x, z) ;
+
+ pos.SetX(x) ;
+ pos.SetZ(z) ;
+ pos.SetY( TMath::Sqrt(x*x + z*z + y0*y0) ) ;
+
+
+
+ Float_t phi = GetPHOSAngle( phosmodule) ;
+ Double_t const kRADDEG = 180.0 / kPI ;
+ Float_t rphi = phi / kRADDEG ;
+
+ TRotation rot ;
+ rot.RotateZ(-rphi) ; // a rotation around Z by angle
+
+ TRotation dummy = rot.Invert() ; // to transform from original frame to rotate frame
+
+ pos.Transform(rot) ; // rotate the baby
}
else {
- pos.SetX(0.);
- pos.SetY(0.);
- pos.SetZ(0.);
- }
+ pos.SetX(0.);
+ pos.SetY(0.);
+ pos.SetZ(0.);
+ }
}
//____________________________________________________________________________
{
// Converts the relative numbering into the local PHOS-module (x, z) coordinates
// Note: sign of z differs from that in the previous version (Yu.Kharlov, 12 Oct 2000)
+
+ Bool_t padOfCPV = (strcmp(fName,"IHEP")==0) ||
+ ((strcmp(fName,"MIXT")==0) && relid[0]<=GetNCPVModules()) ;
+ Bool_t padOfPPSD = (strcmp(fName,"GPS2")==0) ||
+ ((strcmp(fName,"MIXT")==0) && relid[0]> GetNCPVModules()) ;
Int_t ppsdmodule ;
Float_t x0,z0;
Int_t row = relid[2] ; //offset along x axiz
Int_t column = relid[3] ; //offset along z axiz
- Float_t padsizeZ = GetPadSizeZ();
- Float_t padsizeX = GetPadSizePhi();
-
+ Float_t padsizeZ = 0;
+ Float_t padsizeX = 0;
+ Int_t nOfPadsPhi = 0;
+ Int_t nOfPadsZ = 0;
+ if ( padOfPPSD ) {
+ padsizeZ = GetPPSDModuleSize(2) / GetNumberOfPadsZ();
+ padsizeX = GetPPSDModuleSize(0) / GetNumberOfPadsPhi();
+ nOfPadsPhi = GetNumberOfPadsPhi();
+ nOfPadsZ = GetNumberOfPadsZ();
+ }
+ else if ( padOfCPV ) {
+ padsizeZ = GetPadSizeZ();
+ padsizeX = GetPadSizePhi();
+ nOfPadsPhi = GetNumberOfCPVPadsPhi();
+ nOfPadsZ = GetNumberOfCPVPadsZ();
+ }
+
if ( relid[1] == 0 ) { // its a PbW04 crystal
x = - ( GetNPhi()/2. - row + 0.5 ) * GetCrystalSize(0) ; // position ox Xtal with respect
z = ( GetNZ() /2. - column + 0.5 ) * GetCrystalSize(2) ; // of center of PHOS module
}
else {
- if ( relid[1] > GetNumberOfModulesPhi() * GetNumberOfModulesZ() )
- ppsdmodule = relid[1]-GetNumberOfModulesPhi() * GetNumberOfModulesZ();
- else
- ppsdmodule = relid[1] ;
- Int_t modrow = 1+(Int_t)TMath::Ceil( (Float_t)ppsdmodule / GetNumberOfModulesPhi()-1. ) ;
- Int_t modcol = ppsdmodule - ( modrow - 1 ) * GetNumberOfModulesPhi() ;
- if ( ((strcmp( fName, "GPS2" )) == 0) ) {
+ if ( padOfPPSD ) {
+ if ( relid[1] > GetNumberOfModulesPhi() * GetNumberOfModulesZ() )
+ ppsdmodule = relid[1]-GetNumberOfModulesPhi() * GetNumberOfModulesZ();
+ else
+ ppsdmodule = relid[1] ;
+ Int_t modrow = 1+(Int_t)TMath::Ceil( (Float_t)ppsdmodule / GetNumberOfModulesPhi()-1. ) ;
+ Int_t modcol = ppsdmodule - ( modrow - 1 ) * GetNumberOfModulesPhi() ;
x0 = ( GetNumberOfModulesPhi() / 2. - modrow + 0.5 ) * GetPPSDModuleSize(0) ;
z0 = ( GetNumberOfModulesZ() / 2. - modcol + 0.5 ) * GetPPSDModuleSize(2) ;
} else {
x0 = 0;
z0 = 0;
}
- x = - ( GetNumberOfPadsPhi()/2. - row - 0.5 ) * padsizeX + x0 ; // position of pad with respect
- z = ( GetNumberOfPadsZ() /2. - column - 0.5 ) * padsizeZ - z0 ; // of center of PHOS module
+ x = - ( nOfPadsPhi/2. - row - 0.5 ) * padsizeX + x0 ; // position of pad with respect
+ z = ( nOfPadsZ /2. - column - 0.5 ) * padsizeZ - z0 ; // of center of PHOS module
}
}
// Return general PHOS parameters
Int_t GetNModules(void) const { return fNModules ; }
+ Int_t GetNPPSDModules(void) const { return fNPPSDModules ; }
+ Int_t GetNCPVModules(void) const { return fNModules - fNPPSDModules ; }
Float_t GetPHOSAngle(Int_t index) const { return fPHOSAngle[index-1] ; }
// Return EMCA geometrical parameters
// Return PPSD geometrical parameters
- Float_t GetAnodeThickness(void) const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetAnodeThickness(); }
- Float_t GetAvalancheGap(void) const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetAvalancheGap(); }
- Float_t GetCathodeThickness(void) const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetCathodeThickness(); }
- Float_t GetCompositeThickness(void) const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetCompositeThickness(); }
- Float_t GetConversionGap(void) const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetConversionGap(); }
- Float_t GetLeadConverterThickness(void) const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetLeadConverterThickness(); }
- Float_t GetLeadToMicro2Gap(void) const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetLeadToMicro2Gap(); }
- Float_t GetLidThickness(void) const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetLidThickness(); }
- Float_t GetMicromegas1Thickness(void) const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetMicromegas1Thickness(); }
- Float_t GetMicromegas2Thickness(void) const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetMicromegas2Thickness(); }
- Float_t GetMicromegasWallThickness(void) const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetMicromegasWallThickness();}
- Float_t GetMicro1ToLeadGap(void) const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetMicro1ToLeadGap(); }
- Float_t GetPCThickness(void) const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetPCThickness(); }
- Float_t GetPhiDisplacement(void) const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetPhiDisplacement(); }
- Float_t GetPPSDModuleSize(Int_t index) const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetPPSDModuleSize(index); }
- Float_t GetZDisplacement(void) const { return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetZDisplacement(); }
+ Float_t GetAnodeThickness(void) const { return fGeometryPPSD->GetAnodeThickness(); }
+ Float_t GetAvalancheGap(void) const { return fGeometryPPSD->GetAvalancheGap(); }
+ Float_t GetCathodeThickness(void) const { return fGeometryPPSD->GetCathodeThickness(); }
+ Float_t GetCompositeThickness(void) const { return fGeometryPPSD->GetCompositeThickness(); }
+ Float_t GetConversionGap(void) const { return fGeometryPPSD->GetConversionGap(); }
+ Float_t GetLeadConverterThickness(void) const { return fGeometryPPSD->GetLeadConverterThickness(); }
+ Float_t GetLeadToMicro2Gap(void) const { return fGeometryPPSD->GetLeadToMicro2Gap(); }
+ Float_t GetLidThickness(void) const { return fGeometryPPSD->GetLidThickness(); }
+ Float_t GetMicromegas1Thickness(void) const { return fGeometryPPSD->GetMicromegas1Thickness(); }
+ Float_t GetMicromegas2Thickness(void) const { return fGeometryPPSD->GetMicromegas2Thickness(); }
+ Float_t GetMicromegasWallThickness(void) const { return fGeometryPPSD->GetMicromegasWallThickness();}
+ Float_t GetMicro1ToLeadGap(void) const { return fGeometryPPSD->GetMicro1ToLeadGap(); }
+ Int_t GetNumberOfModulesPhi(void) const { return fGeometryPPSD->GetNumberOfModulesPhi(); }
+ Int_t GetNumberOfModulesZ(void) const { return fGeometryPPSD->GetNumberOfModulesZ(); }
+ Int_t GetNumberOfPadsPhi(void) const { return fGeometryPPSD->GetNumberOfPadsPhi(); }
+ Int_t GetNumberOfPadsZ(void) const { return fGeometryPPSD->GetNumberOfPadsZ(); }
+ Float_t GetPCThickness(void) const { return fGeometryPPSD->GetPCThickness(); }
+ Float_t GetPhiDisplacement(void) const { return fGeometryPPSD->GetPhiDisplacement(); }
+ Float_t GetPPSDModuleSize(Int_t index) const { return fGeometryPPSD->GetPPSDModuleSize(index); }
+ Float_t GetZDisplacement(void) const { return fGeometryPPSD->GetZDisplacement(); }
+ void SetLeadConverterThickness(Float_t x) const { fGeometryPPSD->SetLeadConverterThickness(x);}
// Return CPV geometrical parameters
- Bool_t IsLeadConverterExists(void) const { return ((AliPHOSCPVGeometry*) fGeometryCPV)->IsLeadConverterExists(); }
- Float_t GetCPVActiveSize(Int_t index) const { return ((AliPHOSCPVGeometry*) fGeometryCPV)->GetCPVActiveSize(index); }
- Int_t GetNumberOfCPVChipsPhi(void) const { return ((AliPHOSCPVGeometry*) fGeometryCPV)->GetNumberOfCPVChipsPhi(); }
- Int_t GetNumberOfCPVChipsZ(void) const { return ((AliPHOSCPVGeometry*) fGeometryCPV)->GetNumberOfCPVChipsZ(); }
- Float_t GetGassiplexChipSize(Int_t index) const { return ((AliPHOSCPVGeometry*) fGeometryCPV)->GetGassiplexChipSize(index); }
- Float_t GetCPVGasThickness(void) const { return ((AliPHOSCPVGeometry*) fGeometryCPV)->GetCPVGasThickness(); }
- Float_t GetCPVTextoliteThickness(void) const { return ((AliPHOSCPVGeometry*) fGeometryCPV)->GetCPVTextoliteThickness(); }
- Float_t GetCPVCuNiFoilThickness(void) const { return ((AliPHOSCPVGeometry*) fGeometryCPV)->GetCPVCuNiFoilThickness(); }
- Float_t GetFTPosition(Int_t index) const { return ((AliPHOSCPVGeometry*) fGeometryCPV)->GetFTPosition(index); }
- Float_t GetCPVFrameSize(Int_t index) const { return ((AliPHOSCPVGeometry*) fGeometryCPV)->GetCPVFrameSize(index); }
+ Int_t GetNumberOfCPVLayers(void) const { return fGeometryCPV ->GetNumberOfCPVLayers(); }
+ Bool_t IsLeadConverterExists(void) const { return fGeometryCPV->IsLeadConverterExists(); }
+ Float_t GetCPVActiveSize(Int_t index) const { return fGeometryCPV->GetCPVActiveSize(index); }
+ Int_t GetNumberOfCPVChipsPhi(void) const { return fGeometryCPV->GetNumberOfCPVChipsPhi(); }
+ Int_t GetNumberOfCPVChipsZ(void) const { return fGeometryCPV->GetNumberOfCPVChipsZ(); }
+ Int_t GetNumberOfCPVPadsPhi(void) const { return fGeometryCPV->GetNumberOfCPVPadsPhi(); }
+ Int_t GetNumberOfCPVPadsZ(void) const { return fGeometryCPV->GetNumberOfCPVPadsZ(); }
+ Float_t GetPadSizePhi(void) const { return fGeometryCPV->GetCPVPadSizePhi(); }
+ Float_t GetPadSizeZ(void) const { return fGeometryCPV->GetCPVPadSizeZ(); }
+ Float_t GetGassiplexChipSize(Int_t index) const { return fGeometryCPV->GetGassiplexChipSize(index); }
+ Float_t GetCPVGasThickness(void) const { return fGeometryCPV->GetCPVGasThickness(); }
+ Float_t GetCPVTextoliteThickness(void) const { return fGeometryCPV->GetCPVTextoliteThickness(); }
+ Float_t GetCPVCuNiFoilThickness(void) const { return fGeometryCPV->GetCPVCuNiFoilThickness(); }
+ Float_t GetFTPosition(Int_t index) const { return fGeometryCPV->GetFTPosition(index); }
+ Float_t GetCPVFrameSize(Int_t index) const { return fGeometryCPV->GetCPVFrameSize(index); }
// Common PPSD and CPV parameters
- Int_t GetNumberOfCPVLayers(void) const {
- if (strcmp(fName,"GPS2")==0) return 2;
- else if (strcmp(fName,"IHEP")==0) return ((AliPHOSCPVGeometry*) fGeometryCPV)->GetNumberOfCPVLayers();
- else return 0;
- }
-
Float_t GetCPVBoxSize(Int_t index) const {
- if (strcmp(fName,"GPS2")==0) return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetCPVBoxSize(index);
- else if (strcmp(fName,"IHEP")==0) return ((AliPHOSCPVGeometry* ) fGeometryCPV)->GetCPVBoxSize(index);
- else return 0;
- }
-
- Int_t GetNumberOfModulesPhi(void) const {
- if (strcmp(fName,"GPS2")==0) return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetNumberOfModulesPhi();
- else if (strcmp(fName,"IHEP")==0) return 1;
- else return 0;
- }
-
- Int_t GetNumberOfModulesZ(void) const {
- if (strcmp(fName,"GPS2")==0) return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetNumberOfModulesZ();
- else if (strcmp(fName,"IHEP")==0) return 1;
- else return 0;
- }
-
- Int_t GetNumberOfPadsPhi(void) const {
- if (strcmp(fName,"GPS2")==0) return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetNumberOfPadsPhi();
- else if (strcmp(fName,"IHEP")==0) return ((AliPHOSCPVGeometry* ) fGeometryCPV)->GetNumberOfCPVPadsPhi();
- else return 0;
- }
-
- Int_t GetNumberOfPadsZ(void) const {
- if (strcmp(fName,"GPS2")==0) return ((AliPHOSPPSDGeometry*) fGeometryCPV)->GetNumberOfPadsZ();
- else if (strcmp(fName,"IHEP")==0) return ((AliPHOSCPVGeometry* ) fGeometryCPV)->GetNumberOfCPVPadsZ();
- else return 0;
- }
-
- Float_t GetPadSizePhi(void) const {
- if (strcmp(fName,"GPS2")==0) return GetPPSDModuleSize(0) / GetNumberOfPadsPhi();
- else if (strcmp(fName,"IHEP")==0) return ((AliPHOSCPVGeometry*) fGeometryCPV)->GetCPVPadSizePhi();
- else return 0;
- }
-
- Float_t GetPadSizeZ(void) const {
- if (strcmp(fName,"GPS2")==0) return GetPPSDModuleSize(2) / GetNumberOfPadsZ();
- else if (strcmp(fName,"IHEP")==0) return ((AliPHOSCPVGeometry*) fGeometryCPV)->GetCPVPadSizeZ();
+ if (strcmp(fName,"GPS2")==0) return fGeometryPPSD->GetCPVBoxSize(index);
+ else if (strcmp(fName,"IHEP")==0) return fGeometryCPV ->GetCPVBoxSize(index);
+ else if (strcmp(fName,"MIXT")==0) return TMath::Max(fGeometryCPV ->GetCPVBoxSize(index),
+ fGeometryPPSD->GetCPVBoxSize(index));
else return 0;
}
private:
Int_t fNModules ; // Number of modules constituing PHOS
+ Int_t fNPPSDModules ; // Number of PPSD modules
+ Float_t fAngle ; // Position angles between modules
Float_t *fPHOSAngle ; //[fNModules] Position angles of modules
TObjArray *fRotMatrixArray ; // Liste of rotation matrices (one per phos module)
AliPHOSEMCAGeometry *fGeometryEMCA ; // Geometry object for Electromagnetic calorimeter
- AliPHOSCPVBaseGeometry *fGeometryCPV ; // Geometry object for CPV (either GPS2 or IHEP)
+ AliPHOSCPVGeometry *fGeometryCPV ; // Geometry object for CPV (IHEP)
+ AliPHOSPPSDGeometry *fGeometryPPSD ; // Geometry object for PPSD (GPS2)
AliPHOSSupportGeometry *fGeometrySUPP ; // Geometry object for PHOS support
- void SetPHOSAngles(); // calculates the PHOS modules PHI angle
+ void SetPHOSAngles(); // calculates the PHOS modules PHI angle
static AliPHOSGeometry * fgGeom ; // pointer to the unique instance of the singleton
static Bool_t fgInit ; // Tells if geometry has been succesfully set up
Bool_t rv = kFALSE ;
- if ( fId == rValue.GetId() && fPrimary == rValue.GetPrimary() )
+ if ( fId == rValue.GetId() && fPrimary == rValue.GetPrimary() )
rv = kTRUE;
return rv;
// return the identificator of this his
return fId ;
}
+ Int_t GetPid(void) const {
+ // return the particle PDG code which initiates this hit
+ return fPid ;
+ }
Int_t GetPrimary(void) const {
// returns the primary particle id at the origine of this hit
return fPrimary ;
//*-- Author: Yves Schutz (SUBATECH)
//////////////////////////////////////////////////////////////////////////////
+#include <iostream.h>
+
// --- ROOT system ---
#include "TTree.h"
{
// Calculates the distance between the EMC RecPoint and the PPSD RecPoint
- Float_t r;
+ Float_t r = 1.e+10;
TVector3 vecEmc ;
TVector3 vecPpsd ;
virtual Float_t GetCPVBoxSize(Int_t index) { return fPPSDBoxSize[index] ; }
virtual Float_t GetPPSDModuleSize(Int_t index) { return fPPSDModuleSize[index] ; }
virtual Float_t GetZDisplacement(void) { return fZDisplacement ; }
+ virtual void SetLeadConverterThickness(Float_t x) { fLeadConverterThickness = x; }
// CPV functions cannot be used for PPSD
// --- Standard library ---
+#include <assert.h>
+
// --- AliRoot header files ---
#include "AliPHOSRecParticle.h"
//
//*-- Author: Gines Martinez (SUBATECH)
+#include <assert.h>
+
// --- ROOT system ---
#include "TMarker.h"
#include "AliRecPoint.h"
#include "AliPHOSDigit.h"
+#include "AliPHOSGeometry.h"
+class AliPHOSGeometry;
class AliPHOSRecPoint : public AliRecPoint {
return *this ;
}
- protected:
+protected:
Int_t fPHOSMod ; // PHOS Module number in which the RecPoint is found
-
-
+ AliPHOSGeometry * fGeom ; // pointer to the PHOS geometry class
+
ClassDef(AliPHOSRecPoint,1) // RecPoint for PHOS (Base Class)
};
}
//____________________________________________________________________________
- void AliPHOSReconstructioner::Make(DigitsList * dl,
- AliPHOSRecPoint::RecPointsList * emccl,
- AliPHOSRecPoint::RecPointsList * ppsdl,
- AliPHOSTrackSegment::TrackSegmentsList * trsl,
- AliPHOSRecParticle::RecParticlesList * rpl)
+ void AliPHOSReconstructioner::MakePPSD(DigitsList * dl,
+ AliPHOSRecPoint::RecPointsList * emccl,
+ AliPHOSRecPoint::RecPointsList * ppsdl,
+ AliPHOSTrackSegment::TrackSegmentsList * trsl,
+ AliPHOSRecParticle::RecParticlesList * rpl)
{
// Launches the Reconstruction process in the sequence: Make the reconstructed poins (clusterize)
// Make the track segments
}
-
-
// Making Clusters
if (fDebugReconstruction) cout << "DebugReconstruction>>> Start making reconstructed points (clusterizing)" << endl;
}
//____________________________________________________________________________
- void AliPHOSReconstructioner::Make(DigitsList * dl,
- AliPHOSRecPoint::RecPointsList * emccl,
- AliPHOSRecPoint::RecPointsList * cpvcl)
+ void AliPHOSReconstructioner::MakeCPV(DigitsList * dl,
+ AliPHOSRecPoint::RecPointsList * emccl,
+ AliPHOSRecPoint::RecPointsList * cpvcl)
{
// Launches the Reconstruction process of EMC and CPV in the sequence:
fClusterizer->MakeClusters(dl, emccl, cpvcl);
if (fDebugReconstruction){
- // Digit Debuging
+ // Digit Debuging
cout << "AliPHOSReconstructioner: Digit list entries are " << dl->GetEntries() << endl ;
cout << "AliPHOSReconstructioner: EMC list entries are " << emccl->GetEntries() << endl ;
cout << "AliPHOSReconstructioner: CPV list entries are " << cpvcl->GetEntries() << endl ;
AliPHOSCpvRecPoint * cpvrp ;
for (index = 0 ; index < cpvcl->GetEntries() ; index++) {
cpvrp = (AliPHOSCpvRecPoint * )cpvcl->At(index) ;
+ if (cpvrp->IsCPV()) break;
cpvrp ->SetIndexInList(index) ;
TVector3 locpos; cpvrp->GetLocalPosition(locpos);
Int_t lengX,lengZ; cpvrp->GetClusterLengths(lengX,lengZ);
AliPHOSClusterizer * GetClusterizer() { return fClusterizer ; }
void Init(AliPHOSClusterizer * Clusterizer, AliPHOSTrackSegmentMaker * Tracker,
AliPHOSPID * Identifier) ;
- void Make(TClonesArray * DL,
- AliPHOSRecPoint::RecPointsList * emccl,
- AliPHOSRecPoint::RecPointsList * ppsdl,
- AliPHOSTrackSegment::TrackSegmentsList * trsl,
- AliPHOSRecParticle::RecParticlesList * rpl) ; // does the job for EMC+PPSD
- void Make(TClonesArray * DL,
- AliPHOSRecPoint::RecPointsList * emccl,
- AliPHOSRecPoint::RecPointsList * ppsdl) ; // does the job for EMC+CPV
+ void MakePPSD(TClonesArray * DL,
+ AliPHOSRecPoint::RecPointsList * emccl,
+ AliPHOSRecPoint::RecPointsList * ppsdl,
+ AliPHOSTrackSegment::TrackSegmentsList * trsl,
+ AliPHOSRecParticle::RecParticlesList * rpl) ; // does the job for EMC+PPSD
+ void MakeCPV (TClonesArray * DL,
+ AliPHOSRecPoint::RecPointsList * emccl,
+ AliPHOSRecPoint::RecPointsList * ppsdl) ; // does the job for EMC+CPV
void SetDebugReconstruction(Bool_t deb) { fDebugReconstruction = deb; }
{
// Makes the track segments out of the list of EMC and PPSD Recpoints and stores them in a list
- Int_t phosmod = 1 ;
+ Int_t phosmod = fGeom->GetNCPVModules() + 1 ;
Int_t emcStopedAt = 0 ;
Int_t ppsdStopedAt = 0 ;
AliPHOS(name,title)
{
// ctor : title is used to identify the layout
- // GPS2 = 5 modules (EMC + PPSD)
+ // GPS2 = 5 modules (EMC + PPSD)
+ // IHEP = 5 modules (EMC + CPV)
+ // MIXT = 4 modules (EMC + CPV) and 1 module (EMC + PPSD)
// gets an instance of the geometry parameters class
//END_HTML
this->BuildGeometryforPHOS() ;
- if ( ( strcmp(fGeom->GetName(), "GPS2" ) == 0 ) )
+ if (strcmp(fGeom->GetName(),"GPS2") == 0)
this->BuildGeometryforPPSD() ;
- else if ( ( strcmp(fGeom->GetName(), "IHEP" ) == 0 ) )
+ else if (strcmp(fGeom->GetName(),"IHEP") == 0)
this->BuildGeometryforCPV() ;
+ else if (strcmp(fGeom->GetName(),"MIXT") == 0) {
+ this->BuildGeometryforPPSD() ;
+ this->BuildGeometryforCPV() ;
+ }
else
- cout << "AliPHOSv0::BuildGeometry : no charged particle identification system installed" << endl;
+ cout << "AliPHOSv0::BuildGeometry : no charged particle identification system installed: "
+ << "Geometry name = " << fGeom->GetName() << endl;
}
Int_t number = 988 ;
TNode * top = gAlice->GetGeometry()->GetNode("alice") ;
- for( Int_t i = 1; i <= fGeom->GetNModules(); i++ ) { // the number of PHOS modules
+ Int_t firstModule = 0 ;
+ if (strcmp(fGeom->GetName(),"GPS2") == 0)
+ firstModule = 1;
+ else if (strcmp(fGeom->GetName(),"MIXT") == 0)
+ firstModule = fGeom->GetNModules() - fGeom->GetNPPSDModules() + 1;
+
+ for( Int_t i = firstModule; i <= fGeom->GetNModules(); i++ ) { // the number of PHOS modules
Float_t angle = fGeom->GetPHOSAngle(i) ;
- sprintf(rotname, "%s%d", "rotg", number++) ;
+ sprintf(rotname, "%s%d", "rotg", number+i) ;
new TRotMatrix(rotname, rotname, 90, angle, 90, 90 + angle, 0, 0);
top->cd();
sprintf(nodename, "%s%d", "Moduleg", i) ;
Float_t r = fGeom->GetIPtoCPVDistance() + fGeom->GetCPVBoxSize(1) / 2.0 ;
Int_t number = 988 ;
TNode * top = gAlice->GetGeometry()->GetNode("alice") ;
- for( Int_t i = 1; i <= fGeom->GetNModules(); i++ ) { // the number of PHOS modules
+
+ Int_t lastModule = 0 ;
+ if (strcmp(fGeom->GetName(),"IHEP") == 0)
+ lastModule = fGeom->GetNModules();
+ else if (strcmp(fGeom->GetName(),"MIXT") == 0)
+ lastModule = fGeom->GetNModules() - fGeom->GetNPPSDModules();
+
+ for( Int_t i = 1; i <= lastModule; i++ ) { // the number of PHOS modules
// One CPV module
Float_t angle = fGeom->GetPHOSAngle(i) ;
- sprintf(rotname, "%s%d", "rotg", number++) ;
+ sprintf(rotname, "%s%d", "rotg", number+i) ;
new TRotMatrix(rotname, rotname, 90, angle, 90, 90 + angle, 0, 0);
top->cd();
sprintf(nodename, "%s%d", "CPVModule", i) ;
// Get pointer to the array containing media indeces
Int_t *idtmed = fIdtmed->GetArray() - 699 ;
+ // Create a box a PHOS module.
+ // In case of MIXT geometry 2 different boxes are needed
+
Float_t bigbox[3] ;
bigbox[0] = fGeom->GetOuterBoxSize(0) / 2.0 ;
bigbox[1] = ( fGeom->GetOuterBoxSize(1) + fGeom->GetCPVBoxSize(1) ) / 2.0 ;
bigbox[2] = fGeom->GetOuterBoxSize(2) / 2.0 ;
- gMC->Gsvolu("PHOS", "BOX ", idtmed[798], bigbox, 3) ;
+ gMC->Gsvolu("PHOS", "BOX ", idtmed[798], bigbox, 3) ;
+
+ if ( strcmp( fGeom->GetName(),"MIXT") == 0 && fGeom->GetNPPSDModules() > 0)
+ gMC->Gsvolu("PHO1", "BOX ", idtmed[798], bigbox, 3) ;
- this->CreateGeometryforPHOS() ;
+ this->CreateGeometryforPHOS() ;
if ( strcmp( fGeom->GetName(), "GPS2") == 0 )
this->CreateGeometryforPPSD() ;
else if ( strcmp( fGeom->GetName(), "IHEP") == 0 )
this->CreateGeometryforCPV() ;
+ else if ( strcmp( fGeom->GetName(), "MIXT") == 0 ) {
+ this->CreateGeometryforPPSD() ;
+ this->CreateGeometryforCPV() ;
+ }
else
cout << "AliPHOSv0::CreateGeometry : no charged particle identification system installed" << endl;
- this->CreateGeometryforSupport() ;
+
+// this->CreateGeometryforSupport() ;
// --- Position PHOS mdules in ALICE setup ---
Int_t idrotm[99] ;
Double_t const kRADDEG = 180.0 / kPI ;
- for( Int_t i = 1; i <= fGeom->GetNModules(); i++ ) {
+ Int_t lastModule;
+ if (strcmp(fGeom->GetName(),"MIXT") == 0)
+ lastModule = fGeom->GetNModules() - fGeom->GetNPPSDModules();
+ else
+ lastModule = fGeom->GetNModules();
+
+ Int_t i;
+ for( i = 1; i <= lastModule ; i++ ) {
Float_t angle = fGeom->GetPHOSAngle(i) ;
AliMatrix(idrotm[i-1], 90.0, angle, 90.0, 90.0+angle, 0.0, 0.0) ;
Float_t r = fGeom->GetIPtoOuterCoverDistance() + ( fGeom->GetOuterBoxSize(1) + fGeom->GetCPVBoxSize(1) ) / 2.0 ;
- Float_t xP1 = r * TMath::Sin( angle / kRADDEG ) ;
+ Float_t xP1 = r * TMath::Sin( angle / kRADDEG ) ;
Float_t yP1 = -r * TMath::Cos( angle / kRADDEG ) ;
gMC->Gspos("PHOS", i, "ALIC", xP1, yP1, 0.0, idrotm[i-1], "ONLY") ;
} // for GetNModules
+ for( i = lastModule+1; i <= fGeom->GetNModules(); i++ ) {
+
+ Float_t angle = fGeom->GetPHOSAngle(i) ;
+ AliMatrix(idrotm[i-1], 90.0, angle, 90.0, 90.0+angle, 0.0, 0.0) ;
+
+ Float_t r = fGeom->GetIPtoOuterCoverDistance() + ( fGeom->GetOuterBoxSize(1) + fGeom->GetCPVBoxSize(1) ) / 2.0 ;
+
+ Float_t xP1 = r * TMath::Sin( angle / kRADDEG ) ;
+ Float_t yP1 = -r * TMath::Cos( angle / kRADDEG ) ;
+
+ gMC->Gspos("PHO1", i-lastModule, "ALIC", xP1, yP1, 0.0, idrotm[i-1], "ONLY") ;
+
+ } // for GetNModules
+
}
//____________________________________________________________________________
Float_t yO = - fGeom->GetCPVBoxSize(1) / 2.0 ;
- gMC->Gspos("EMCA", 1, "PHOS", 0.0, yO, 0.0, 0, "ONLY") ;
+ gMC->Gspos("EMCA", 1, "PHOS", 0.0, yO, 0.0, 0, "ONLY") ;
+ if ( strcmp( fGeom->GetName(),"MIXT") == 0 && fGeom->GetNPPSDModules() > 0)
+ gMC->Gspos("EMCA", 1, "PHO1", 0.0, yO, 0.0, 0, "ONLY") ;
// ---
// --- Define Textolit Wall box, position inside EMCA ---
Float_t yO = fGeom->GetOuterBoxSize(1) / 2.0 ;
- gMC->Gspos("PPSD", 1, "PHOS", 0.0, yO, 0.0, 0, "ONLY") ;
+ if ( strcmp( fGeom->GetName(),"MIXT") == 0 && fGeom->GetNPPSDModules() > 0)
+ gMC->Gspos("PPSD", 1, "PHO1", 0.0, yO, 0.0, 0, "ONLY") ;
+ else
+ gMC->Gspos("PPSD", 1, "PHOS", 0.0, yO, 0.0, 0, "ONLY") ;
// Now we build a micromegas module
// The box containing the whole module filled with epoxy (FR4)
// --- Divide GGPP in X (phi) and Z directions --
gMC->Gsdvn("GROW", "GGPS", fGeom->GetNumberOfPadsPhi(), 1) ;
- gMC->Gsdvn("GCEL", "GROW", fGeom->GetNumberOfPadsZ() , 3) ;
+ gMC->Gsdvn("GCEL", "GROW", fGeom->GetNumberOfPadsZ() , 3) ;
y0 = y0 - fGeom->GetAnodeThickness() / 2. - ( fGeom->GetConversionGap() + fGeom->GetAvalancheGap() ) / 2. ;
}
+//____________________________________________________________________________
+Float_t AliPHOSv0::ZMin(void) const
+{
+ // Overall dimension of the PHOS (min)
+ // Take it twice more than the PHOS module size
+ return -fGeom->GetOuterBoxSize(2);
+}
+
+//____________________________________________________________________________
+Float_t AliPHOSv0::ZMax(void) const
+{
+ // Overall dimension of the PHOS (max)
+ // Take it twice more than the PHOS module size
+ return fGeom->GetOuterBoxSize(2);
+}
+
//____________________________________________________________________________
void AliPHOSv0::Init(void)
{
// useless since there are no hits
assert(0==1) ;
}
- virtual void BuildGeometry(void) ; // creates the geometry for the ROOT display
- void BuildGeometryforPHOS(void) ; // creates the PHOS geometry for the ROOT display
- void BuildGeometryforPPSD(void) ; // creates the PPSD geometry for the ROOT display
- void BuildGeometryforCPV(void) ; // creates the CPV geometry for the ROOT display
- virtual void CreateGeometry(void) ; // creates the geometry for GEANT
- void CreateGeometryforPHOS(void) ; // creates the PHOS geometry for GEANT
- void CreateGeometryforPPSD(void) ; // creates the PPSD geometry for GEANT
- void CreateGeometryforCPV(void) ; // creates the CPV geometry for GEANT
- void CreateGeometryforSupport(void) ; // creates the Support geometry for GEANT
+ virtual void BuildGeometry(void) ; // creates the geometry for the ROOT display
+ void BuildGeometryforPHOS(void) ; // creates the PHOS geometry for the ROOT display
+ void BuildGeometryforPPSD(void) ; // creates the PPSD geometry for the ROOT display
+ void BuildGeometryforCPV(void) ; // creates the CPV geometry for the ROOT display
+ virtual void CreateGeometry(void) ; // creates the geometry for GEANT
+ void CreateGeometryforPHOS(void) ; // creates the PHOS geometry for GEANT
+ void CreateGeometryforPPSD(void) ; // creates the PPSD geometry for GEANT
+ void CreateGeometryforCPV(void) ; // creates the CPV geometry for GEANT
+ void CreateGeometryforSupport(void) ; // creates the Support geometry for GEANT
+ virtual Float_t ZMin() const; // overall dimension of the module (min)
+ virtual Float_t ZMax() const; // overall dimension of the module (max)
+
virtual AliPHOSGeometry * GetGeometry() {
// gets the pointer to the AliPHOSGeometry unique instance
return fGeom ;
// Implementation version v1 of PHOS Manager class
//---
// Layout EMC + PPSD has name GPS2:
-// Produces cumulated hits (no hits) and digits
+// Produces cumulated hits
//---
// Layout EMC + CPV has name IHEP:
-// Produces hits for CPV, cumulated hits and digits
+// Produces hits for CPV, cumulated hits
+//---
+// Layout EMC + CPV + PPSD has name GPS:
+// Produces hits for CPV, cumulated hits
+//---
//*-- Author: Yves Schutz (SUBATECH)
AliPHOSv1::AliPHOSv1()
{
// ctor
- fNTmpHits = 0 ;
- fTmpHits = 0 ;
// Create an empty array of AliPHOSCPVModule to satisfy
// AliPHOSv1::Streamer when reading root file
+
+ fReconstructioner = 0;
+ fTrackSegmentMaker = 0;
+
+ if ( 0==(fEMCModules=new TClonesArray("AliPHOSCPVModule",0)) ) {
+ Error("AliPHOSv1","Can not create array of EMC modules");
+ exit(1);
+ }
- if ( NULL==(fCPVModules=new TClonesArray("AliPHOSCPVModule",0)) ) {
+ if ( 0==(fCPVModules=new TClonesArray("AliPHOSCPVModule",0)) ) {
Error("AliPHOSv1","Can not create array of CPV modules");
exit(1);
}
// ctor : title is used to identify the layout
// GPS2 = 5 modules (EMC + PPSD)
// IHEP = 5 modules (EMC + CPV )
- // We use 2 arrays of hits :
+ // MIXT = 4 modules (EMC + CPV ) and 1 module (EMC + PPSD)
//
+ // We store hits :
// - fHits (the "normal" one), which retains the hits associated with
// the current primary particle being tracked
// (this array is reset after each primary has been tracked).
//
- // - fTmpHits, which retains all the hits of the current event. It
- // is used for the digitization part.
fPinElectronicNoise = 0.010 ;
fDigitThreshold = 0.1 ; // 1 GeV
// It is put in the Digit Tree because the TreeH is filled after each primary
// and the TreeD at the end of the event (branch is set in FinishEvent() ).
- fTmpHits= new TClonesArray("AliPHOSHit",1000) ;
+ fHits= new TClonesArray("AliPHOSHit",1000) ;
- fNTmpHits = fNhits = 0 ;
+ fNhits = 0 ;
- fDigits = new TClonesArray("AliPHOSDigit",1000) ;
+ fReconstructioner = 0;
+ fTrackSegmentMaker = 0;
fIshunt = 1 ; // All hits are associated with primary particles
+ // Create array of EMC modules of the size of PHOS modules number
+
+ if ( 0==(fEMCModules=new TClonesArray("AliPHOSCPVModule",fGeom->GetNModules())) ) {
+ Error("AliPHOSv1","Can not create array of EMC modules");
+ exit(1);
+ }
+ TClonesArray &lemcmodule = *fEMCModules;
+ for (Int_t i=0; i<fGeom->GetNModules(); i++) new(lemcmodule[i]) AliPHOSCPVModule();
+
// Create array of CPV modules for the IHEP's version of CPV
- if ( strcmp(fGeom->GetName(),"IHEP") == 0 ) {
- // Create array of AliPHOSCPVmodule of the size of PHOS modules number
+ if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
+ // Create array of CPV modules of the size of PHOS modules number
- if ( NULL==(fCPVModules=new TClonesArray("AliPHOSCPVModule",fGeom->GetNModules())) ) {
+ if ( 0==(fCPVModules=new TClonesArray("AliPHOSCPVModule",fGeom->GetNCPVModules())) ) {
Error("AliPHOSv1","Can not create array of CPV modules");
exit(1);
}
TClonesArray &lcpvmodule = *fCPVModules;
- for (Int_t i=0; i<fGeom->GetNModules(); i++) new(lcpvmodule[i]) AliPHOSCPVModule();
+ for (Int_t i=0; i<fGeom->GetNCPVModules(); i++) new(lcpvmodule[i]) AliPHOSCPVModule();
}
else {
// Create an empty array of AliPHOSCPVModule to satisfy
// We do not want to save in TreeH the raw hits
//fHits = new TClonesArray("AliPHOSHit",100) ;
- fDigits = new TClonesArray("AliPHOSDigit",1000) ;
- fTmpHits= new TClonesArray("AliPHOSHit",1000) ;
+ fDigits = 0 ;
+ fHits= new TClonesArray("AliPHOSHit",1000) ;
- fNTmpHits = fNhits = 0 ;
+ fNhits = 0 ;
fIshunt = 1 ; // All hits are associated with primary particles
{
// dtor
- if ( fTmpHits) {
- fTmpHits->Delete() ;
- delete fTmpHits ;
- fTmpHits = 0 ;
+ if ( fHits) {
+ fHits->Delete() ;
+ delete fHits ;
+ fHits = 0 ;
}
-// if ( fEmcRecPoints ) {
-// fEmcRecPoints->Delete() ;
-// delete fEmcRecPoints ;
-// fEmcRecPoints = 0 ;
-// }
-
-// if ( fPpsdRecPoints ) {
-// fPpsdRecPoints->Delete() ;
-// delete fPpsdRecPoints ;
-// fPpsdRecPoints = 0 ;
-// }
-
-// if ( fTrackSegments ) {
-// fTrackSegments->Delete() ;
-// delete fTrackSegments ;
-// fTrackSegments = 0 ;
-// }
+ if ( fDigits) {
+ fDigits->Delete() ;
+ delete fDigits ;
+ fDigits = 0 ;
+ }
+ if ( fEmcRecPoints ) {
+ fEmcRecPoints->Delete() ;
+ delete fEmcRecPoints ;
+ fEmcRecPoints = 0 ;
+ }
+
+ if ( fPpsdRecPoints ) {
+ fPpsdRecPoints->Delete() ;
+ delete fPpsdRecPoints ;
+ fPpsdRecPoints = 0 ;
+ }
+
+ if ( fTrackSegments ) {
+ fTrackSegments->Delete() ;
+ delete fTrackSegments ;
+ fTrackSegments = 0 ;
+ }
+
}
//____________________________________________________________________________
// or in a single PPSD gas cell
Int_t hitCounter ;
- TClonesArray <mphits = *fTmpHits ;
AliPHOSHit *newHit ;
AliPHOSHit *curHit ;
Bool_t deja = kFALSE ;
- // In any case, fills the fTmpHit TClonesArray (with "accumulated hits")
-
newHit = new AliPHOSHit(shunt, primary, tracknumber, Id, hits, trackpid) ;
- // We do not want to save in TreeH the raw hits
- // TClonesArray &lhits = *fHits;
-
- for ( hitCounter = 0 ; hitCounter < fNTmpHits && !deja ; hitCounter++ ) {
- curHit = (AliPHOSHit*) ltmphits[hitCounter] ;
- if( *curHit == *newHit ) {
- *curHit = *curHit + *newHit ;
- deja = kTRUE ;
+ for ( hitCounter = 0 ; hitCounter < fNhits && !deja ; hitCounter++ ) {
+ curHit = (AliPHOSHit*) (*fHits)[hitCounter] ;
+ if( *curHit == *newHit ) {
+ *curHit = *curHit + *newHit ;
+ deja = kTRUE ;
}
}
if ( !deja ) {
- new(ltmphits[fNTmpHits]) AliPHOSHit(*newHit) ;
- fNTmpHits++ ;
+ new((*fHits)[fNhits]) AliPHOSHit(*newHit) ;
+ fNhits++ ;
}
- // We do not want to save in TreeH the raw hits
- // new(lhits[fNhits]) AliPHOSHit(*newHit) ;
- // fNhits++ ;
-
- // Please note that the fTmpHits array must survive up to the
- // end of the events, so it does not appear e.g. in ResetHits() (
- // which is called at the end of each primary).
-
delete newHit;
-
}
//___________________________________________________________________________
return chan ;
}
-//___________________________________________________________________________
-void AliPHOSv1::FinishEvent()
-{
- // Makes the digits from the sum of summed hit in a single crystal or PPSD gas cell
- // Adds to the energy the electronic noise
- // Keeps digits with energy above fDigitThreshold
-
- // Save the cumulated hits instead of raw hits (need to create the branch myself)
- // It is put in the Digit Tree because the TreeH is filled after each primary
- // and the TreeD at the end of the event.
+//____________________________________________________________________________
+void AliPHOSv1::Hit2Digit(Int_t ntracks){
+ //Collects all hits in the same active volume into digits
+
+ if(fDigits!= 0)
+ fDigits->Clear() ;
+ else
+ fDigits = new TClonesArray("AliPHOSDigit",1000) ;
+
+ // Branch address for digit tree
+ char branchname[20];
+ sprintf(branchname,"%s",GetName());
+ gAlice->TreeD()->Branch(branchname,&fDigits,fBufferSize);
+ gAlice->TreeD()->GetEvent(0);
+
Int_t i ;
Int_t relid[4];
Int_t j ;
- TClonesArray &lDigits = *fDigits ;
AliPHOSHit * hit ;
AliPHOSDigit * newdigit ;
AliPHOSDigit * curdigit ;
Bool_t deja = kFALSE ;
- for ( i = 0 ; i < fNTmpHits ; i++ ) {
- hit = (AliPHOSHit*)fTmpHits->At(i) ;
-
- // Assign primary number only if contribution is significant
- if( hit->GetEnergy() > fDigitThreshold)
- newdigit = new AliPHOSDigit( hit->GetPrimary(), hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
- else
- newdigit = new AliPHOSDigit( -1 , hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
- deja =kFALSE ;
- for ( j = 0 ; j < fNdigits ; j++) {
- curdigit = (AliPHOSDigit*) lDigits[j] ;
- if ( *curdigit == *newdigit) {
- *curdigit = *curdigit + *newdigit ;
- deja = kTRUE ;
+ Int_t itrack ;
+ for (itrack=0; itrack<ntracks; itrack++){
+
+ //=========== Get the Hits Tree for the Primary track itrack
+ gAlice->ResetHits();
+ gAlice->TreeH()->GetEvent(itrack);
+
+ for ( i = 0 ; i < fHits->GetEntries() ; i++ ) {
+ hit = (AliPHOSHit*)fHits->At(i) ;
+
+ // Assign primary number only if contribution is significant
+ if( hit->GetEnergy() > fDigitThreshold)
+ newdigit = new AliPHOSDigit( hit->GetPrimary(), hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
+ else
+ newdigit = new AliPHOSDigit( -1 , hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
+ deja =kFALSE ;
+
+
+ for ( j = 0 ; j < fNdigits ; j++) {
+ curdigit = (AliPHOSDigit*) fDigits->At(j) ;
+ if ( *curdigit == *newdigit) {
+ *curdigit = *curdigit + *newdigit ;
+ deja = kTRUE ;
+ }
}
- }
- if ( !deja ) {
- new(lDigits[fNdigits]) AliPHOSDigit(* newdigit) ;
- fNdigits++ ;
- }
-
- delete newdigit ;
- }
-
- // Noise induced by the PIN diode of the PbWO crystals
+ if ( !deja ) {
+ new((*fDigits)[fNdigits]) AliPHOSDigit(* newdigit) ;
+ fNdigits++ ;
+ }
+
+ delete newdigit ;
+ }
+
+ } // loop over tracks
+
+ // Noise induced by the PIN diode of the PbWO crystals
+
Float_t energyandnoise ;
for ( i = 0 ; i < fNdigits ; i++ ) {
newdigit = (AliPHOSDigit * ) fDigits->At(i) ;
-
+
fGeom->AbsToRelNumbering(newdigit->GetId(), relid) ;
-
+
if (relid[1]==0){ // Digits belong to EMC (PbW0_4 crystals)
energyandnoise = newdigit->GetAmp() + Digitize(gRandom->Gaus(0., fPinElectronicNoise)) ;
-
+
if (energyandnoise < 0 )
energyandnoise = 0 ;
-
+
if ( newdigit->GetAmp() < fDigitThreshold ) // if threshold not surpassed, remove digit from list
fDigits->RemoveAt(i) ;
}
}
-
+
fDigits->Compress() ;
+
+ fNdigits = fDigits->GetEntries() ;
+ fDigits->Expand(fNdigits) ;
- fNdigits = fDigits->GetEntries() ;
for (i = 0 ; i < fNdigits ; i++) {
newdigit = (AliPHOSDigit *) fDigits->At(i) ;
- newdigit->SetIndexInList(i) ;
-
-// fGeom->AbsToRelNumbering(newdigit->GetId(), relid) ;
-// printf("FinishEvent(): relid=(%d,%d,%d,%d) Amp=%d\n",
-// relid[0],relid[1],relid[2],relid[3], newdigit->GetAmp());
+ newdigit->SetIndexInList(i) ;
}
-}
+ gAlice->TreeD()->Fill() ;
+ gAlice->TreeD()->Write(0,TObject::kOverwrite) ;
+
+}
//___________________________________________________________________________
void AliPHOSv1::MakeBranch(Option_t* opt)
{
// Create new branche in the current Root Tree in the digit Tree
AliDetector::MakeBranch(opt) ;
- char branchname[10];
-
- // Create new branche PHOSCH in the current Root Tree in the Hit Tree for accumulated Hits
- if ( ! (gAlice->IsLegoRun()) ) { // only when not in lego plot mode
- char *cdH = strstr(opt,"H");
- if ( fTmpHits && gAlice->TreeH() && cdH) {
- sprintf(branchname, "%sCH", GetName()) ;
- gAlice->TreeH()->Branch(branchname, &fTmpHits, fBufferSize) ;
- }
- }
- sprintf(branchname,"%s",GetName());
- char *cdD = strstr(opt,"D");
- if (fDigits && gAlice->TreeD() && cdD) {
- gAlice->TreeD()->Branch(branchname, &fDigits, fBufferSize);
- }
-
-
+ // Create new branches EMC<i> for hits in EMC modules
+
+ for( Int_t i=0; i<fGeom->GetNModules(); i++ ) GetEMCModule(i).MakeBranch("EMC",i+1);
+
// Create new branches CPV<i> for hits in CPV modules for IHEP geometry
- // Yuri Kharlov, 28 September 2000.
-
- if ( strcmp(fGeom->GetName(),"IHEP") == 0 ) {
- for( Int_t i=0; i<fGeom->GetNModules(); i++ ) GetCPVModule(i).MakeBranch(i+1);
+
+ if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
+ for( Int_t i=0; i<fGeom->GetNCPVModules(); i++ ) GetCPVModule(i).MakeBranch("CPV",i+1);
}
-
+
}
//_____________________________________________________________________________
fRecParticles->Delete() ;
- if (strcmp(fGeom->GetName(),"GPS2") == 0) {
+ if (strcmp(fGeom->GetName(),"GPS2") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0) {
if ( fRecParticles && gAlice->TreeR() ) {
sprintf(branchname,"%sRP",GetName()) ;
gAlice->TreeR()->Branch(branchname, &fRecParticles, fBufferSize) ;
}
// 3.
- if (strcmp(fGeom->GetName(),"GPS2") == 0)
- fReconstructioner->Make(fDigits, fEmcRecPoints, fPpsdRecPoints, fTrackSegments, fRecParticles);
- else if (strcmp(fGeom->GetName(),"IHEP") == 0)
- fReconstructioner->Make(fDigits, fEmcRecPoints, fPpsdRecPoints);
+ if (strcmp(fGeom->GetName(),"GPS2") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0)
+ fReconstructioner->MakePPSD(fDigits, fEmcRecPoints, fPpsdRecPoints, fTrackSegments, fRecParticles);
+ if (strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0)
+ fReconstructioner->MakeCPV (fDigits, fEmcRecPoints, fPpsdRecPoints);
// 4. Expand or Shrink the arrays to the proper size
size = fEmcRecPoints->GetEntries() ;
fEmcRecPoints->Expand(size) ;
-
+
size = fPpsdRecPoints->GetEntries() ;
fPpsdRecPoints->Expand(size) ;
//____________________________________________________________________________
void AliPHOSv1::ResetHits()
{
- // Reset hit tree for CPV in IHEP geometry
+ // Reset hit tree for EMC and CPV
// Yuri Kharlov, 28 September 2000
-
- if ( fTmpHits ) {
- fTmpHits->Clear() ;
- fNTmpHits = 0;
- }
-
AliDetector::ResetHits();
- if ( strcmp(fGeom->GetName(),"IHEP") == 0 ) {
- for (Int_t i=0; i<fGeom->GetNModules(); i++) ((AliPHOSCPVModule*)(*fCPVModules)[i]) -> Clear();
- }
-}
-//____________________________________________________________________________
-void AliPHOSv1::ResetDigits()
-{
- // May sound strange, but cumulative hits are store in digits Tree
- AliDetector::ResetDigits();
- if( fTmpHits ) {
- fTmpHits->Delete();
- fNTmpHits = 0 ;
+ for (Int_t i=0; i<fGeom->GetNModules(); i++) ((AliPHOSCPVModule*)(*fEMCModules)[i]) -> Clear();
+ if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
+ for (Int_t i=0; i<fGeom->GetNCPVModules(); i++) ((AliPHOSCPVModule*)(*fCPVModules)[i]) -> Clear();
}
+
}
//____________________________________________________________________________
void AliPHOSv1::ResetReconstruction()
// TBranch *branch;
AliPHOS::SetTreeAddress();
- // //Branch address for TreeR: RecPpsdRecPoint
+// //Branch address for TreeR: RecPpsdRecPoint
// TTree *treeR = gAlice->TreeR();
// if ( treeR && fPpsdRecPoints ) {
// branch = treeR->GetBranch("PHOSPpsdRP");
// if (branch) branch->SetAddress(&fPpsdRecPoints) ;
// }
+ // Set branch address for the Hits Tree for hits in EMC modules
+ // Yuri Kharlov, 23 November 2000.
+
+ for( Int_t i=0; i<fGeom->GetNModules(); i++ ) GetEMCModule(i).SetTreeAddress("EMC",i+1);
+
// Set branch address for the Hits Tree for hits in CPV modules for IHEP geometry
// Yuri Kharlov, 28 September 2000.
- if ( strcmp(fGeom->GetName(),"IHEP") == 0 ) {
- for( Int_t i=0; i<fGeom->GetNModules(); i++ ) GetCPVModule(i).SetTreeAddress(i+1);
+ if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
+ for( Int_t i=0; i<fGeom->GetNCPVModules(); i++ ) GetCPVModule(i).SetTreeAddress("CPV",i+1);
}
}
{
// Accumulates hits as long as the track stays in a single crystal or PPSD gas Cell
+// if (gMC->IsTrackEntering())
+// cout << "Track enters the volume " << gMC->CurrentVolName() << endl;
+// if (gMC->IsTrackExiting())
+// cout << "Track leaves the volume " << gMC->CurrentVolName() << endl;
+
Int_t relid[4] ; // (box, layer, row, column) indices
Int_t absid ; // absolute cell ID number
Float_t xyze[4] ; // position wrt MRS and energy deposited
TLorentzVector pos ; // Lorentz vector of the track current position
Int_t copy ;
+ Int_t i ;
Int_t tracknumber = gAlice->CurrentTrack() ;
Int_t primary = gAlice->GetPrimary( gAlice->CurrentTrack() );
TString name = fGeom->GetName() ;
Int_t trackpid = gMC->TrackPid() ;
- if ( name == "GPS2" ) { // ======> CPV is a GPS' PPSD
+ if ( name == "GPS2" || name == "MIXT" ) { // ======> CPV is a GPS' PPSD
if( gMC->CurrentVolID(copy) == gMC->VolId("GCEL") ) // We are inside a gas cell
{
if ( xyze[3] != 0 ) { // there is deposited energy
gMC->CurrentVolOffID(5, relid[0]) ; // get the PHOS Module number
+ if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(5),"PHO1") == 0 )
+ relid[0] += fGeom->GetNModules() - fGeom->GetNPPSDModules();
gMC->CurrentVolOffID(3, relid[1]) ; // get the Micromegas Module number
- // 1-> Geom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() upper
- // > fGeom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() lower
+ // 1-> fGeom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() upper
+ // > fGeom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() lower
gMC->CurrentVolOffID(1, relid[2]) ; // get the row number of the cell
gMC->CurrentVolID(relid[3]) ; // get the column number
} // We are inside the gas of the CPV
} // GPS2 configuration
- else if ( name == "IHEP" ) { // ======> CPV is a IHEP's one
+ if ( name == "IHEP" || name == "MIXT" ) { // ======> CPV is a IHEP's one
// Yuri Kharlov, 28 September 2000
// Add the current particle in the list of the CPV hits.
- phos.GetCPVModule(moduleNumber).AddHit(fIshunt,tracknumber,pmom,xyd,ipart);
+ phos.GetCPVModule(moduleNumber).AddHit(fIshunt,primary,pmom,xyd,ipart);
- if (fDebugLevel == 1) {
+ if (fDebug == 1) {
printf("CPV hit added to module #%2d: p = (% .4f, % .4f, % .4f, % .4f) GeV,\n",
moduleNumber+1,pmom.Px(),pmom.Py(),pmom.Pz(),pmom.E());
- printf( " xy = (%8.4f, %8.4f) cm, ipart = %d\n",
- xyd[0],xyd[1],ipart);
+ printf( " xy = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
+ xyd[0],xyd[1],ipart,primary);
}
// Digitize the current CPV hit:
xyze[1] = pos[1] ;
xyze[2] = pos[2] ;
xyze[3] = gMC->Edep() ;
-
+
+ // Track enters to the crystal from the top edge
+
+ if (gMC->IsTrackEntering()) {
+ Float_t posloc[3];
+ gMC -> Gmtod (xyze, posloc, 1);
+ if (posloc[1] > fGeom->GetCrystalSize(1)/2-0.01) {
+ Int_t row,cel;
+ Float_t xyd[2];
+ AliPHOSv1 &phos = *(AliPHOSv1*)gAlice->GetModule("PHOS");
+
+ Int_t moduleNumber;
+ gMC->CurrentVolOffID(10,moduleNumber);
+ if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(10),"PHO1") == 0 )
+ moduleNumber += fGeom->GetNModules() - fGeom->GetNPPSDModules();
+ moduleNumber--;
+
+ gMC->CurrentVolOffID(4, row) ; // get the row number inside the module
+ gMC->CurrentVolOffID(3, cel) ; // get the cell number inside the module
+ xyd[0] = -(posloc[2] + (cel-0.5-fGeom->GetNZ() /2) *
+ (fGeom->GetCrystalSize(2) + 2 * fGeom->GetGapBetweenCrystals()));
+ xyd[1] = posloc[0] + (row-0.5-fGeom->GetNPhi()/2) *
+ (fGeom->GetCrystalSize(0) + 2 * fGeom->GetGapBetweenCrystals());
+
+ // Current momentum of the hit's track in the CPV module ref. system
+
+ TLorentzVector pmom;
+ gMC -> TrackMomentum(pmom);
+ Float_t pm[3], pd[3];
+ for (i=0; i<3; i++) pm[i] = pmom[i];
+ gMC -> Gmtod (pm, pd, 2); // transform 3-momentum from master to daughter system
+ pmom[0] = pd[0];
+ pmom[1] =-pd[1];
+ pmom[2] =-pd[2];
+
+ // Current particle type of the hit's track
+
+ Int_t ipart = gMC->TrackPid();
+
+ // Add the current particle in the list of the EMC hits.
+
+ phos.GetEMCModule(moduleNumber).AddHit(fIshunt,primary,pmom,xyd,ipart);
+
+ if (fDebug == 1) {
+ printf("EMC hit added to module #%2d: p = (% .4f, % .4f, % .4f, % .4f) GeV,\n",
+ moduleNumber+1,pmom.Px(),pmom.Py(),pmom.Pz(),pmom.E());
+ printf( " xy = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
+ xyd[0],xyd[1],ipart,primary);
+ }
+ }
+ }
+
+ // Track is inside the crystal and deposits some energy
+
if ( xyze[3] != 0 ) {
gMC->CurrentVolOffID(10, relid[0]) ; // get the PHOS module number ;
+ if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(5),"PHO1") == 0 )
+ relid[0] += fGeom->GetNModules() - fGeom->GetNPPSDModules();
relid[1] = 0 ; // means PBW04
gMC->CurrentVolOffID(4, relid[2]) ; // get the row number inside the module
gMC->CurrentVolOffID(3, relid[3]) ; // get the cell number inside the module
// get the absolute Id number
fGeom->RelToAbsNumbering(relid, absid) ;
-
+
// add current hit to the hit list
AddHit(fIshunt, primary,tracknumber, absid, xyze, trackpid);
-
+
} // there is deposited energy
} // we are inside a PHOS Xtal
}
// Finite size of ionization region
- Int_t nCellZ = fGeom->GetNumberOfPadsZ();
- Int_t nCellX = fGeom->GetNumberOfPadsPhi();
+ Int_t nCellZ = fGeom->GetNumberOfCPVPadsZ();
+ Int_t nCellX = fGeom->GetNumberOfCPVPadsPhi();
Int_t nz3 = (kNgamz+1)/2;
Int_t nx3 = (kNgamx+1)/2;
cpvDigits->Expand(nIter*kNgamx*kNgamz);
virtual void AddHit( Int_t shunt, Int_t primary, Int_t track, Int_t id, Float_t *hits, Int_t pid) ;
Int_t Digitize(Float_t Energy);
- virtual void FinishEvent(void) ;
+ virtual void Hit2Digit(Int_t event) ;
virtual Int_t IsVersion(void) const {
// Gives the version number
return 1 ;
virtual void MakeBranch(Option_t* opt) ;
void Reconstruction(AliPHOSReconstructioner * Reconstructioner) ;
void ResetClusters(){} ;
- virtual void ResetDigits() ;
virtual void ResetHits() ;
virtual void ResetReconstruction() ; // Reset reconstructed objects
void SetReconstructioner(AliPHOSReconstructioner& Reconstructioner) {
return *this ;
}
- // IHEP's CPV specific functions
-
+ AliPHOSCPVModule &GetEMCModule(int n) { return *(AliPHOSCPVModule*)fEMCModules->operator[](n); }
AliPHOSCPVModule &GetCPVModule(int n) { return *(AliPHOSCPVModule*)fCPVModules->operator[](n); }
+
void CPVDigitize (TLorentzVector p, Float_t *xy, Int_t moduleNumber, TClonesArray *digits) ;
Float_t CPVPadResponseFunction(Float_t qhit, Float_t zg, Float_t xg) ;
Double_t CPVCumulPadResponse(Double_t x, Double_t y) ;
protected:
Float_t fDigitThreshold ; // Threshold for the digit registration
- Int_t fNTmpHits ; //! Used internally for digitalization
Float_t fPinElectronicNoise ; // Electronic Noise in the PIN
- AliPHOSReconstructioner * fReconstructioner ; // Reconstrutioner of the PHOS event: Clusterization and subtracking procedures
- TClonesArray * fTmpHits ; //! Used internally for digitalization
+ AliPHOSReconstructioner * fReconstructioner ; // Reconstrutioner of the PHOS event: Clusterization and subtracking procedures
AliPHOSTrackSegmentMaker * fTrackSegmentMaker ; // Reconstructioner of the PHOS track segment: 2 x PPSD + 1 x EMC
-
- TClonesArray * fCPVModules; // Array of CPV modules for the IHEP's version of CPV
+ TClonesArray * fEMCModules; // Array of EMC modules
+ TClonesArray * fCPVModules; // Array of CPV modules for the IHEP's version of CPV
ClassDef(AliPHOSv1,1) // Implementation of PHOS manager class for layout EMC+PPSD
// AddHit, StepManager,and FinishEvent are redefined
//
//*-- Author: Gines MARTINEZ (SUBATECH)
-
+//*-- Modified Nov. 22 2000 by Dmitri Peressounko
+// All hits are stored.
+// Note, that primaries will not be assigned to digits:
+// becouse of tiny energy deposition at each step.
+//
// --- ROOT system ---
AliPHOSv2::AliPHOSv2()
{
// default ctor
- fNTmpHits = 0 ;
- fTmpHits = 0 ;
+
}
//____________________________________________________________________________
AliPHOSv1(name,title)
{
// ctor
-
- fHits= new TClonesArray("AliPHOSHit",1000) ;
}
//____________________________________________________________________________
AliPHOSv2::~AliPHOSv2()
{
// dtor
- if ( fTmpHits) {
- fTmpHits->Delete() ;
- delete fTmpHits ;
- fTmpHits = 0 ;
- }
-
- if ( fEmcRecPoints ) {
- fEmcRecPoints->Delete() ;
- delete fEmcRecPoints ;
- fEmcRecPoints = 0 ;
- }
-
- if ( fPpsdRecPoints ) {
- fPpsdRecPoints->Delete() ;
- delete fPpsdRecPoints ;
- fPpsdRecPoints = 0 ;
- }
-
- if ( fTrackSegments ) {
- fTrackSegments->Delete() ;
- delete fTrackSegments ;
- fTrackSegments = 0 ;
- }
}
//____________________________________________________________________________
void AliPHOSv2::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id, Float_t * hits, Int_t pid)
{
// Add a hit to the hit list.
- // A PHOS hit is the sum of all hits in a single crystal
- // or in a single PPSD gas cell
- Int_t hitCounter ;
- TClonesArray <mphits = *fTmpHits ;
AliPHOSHit *newHit ;
- AliPHOSHit *curHit ;
- Bool_t deja = kFALSE ;
-
- // In any case, fills the fTmpHit TClonesArray (with "accumulated hits")
newHit = new AliPHOSHit(shunt, primary, tracknumber, Id, hits, pid) ;
- // We do want to save in TreeH the raw hits
- TClonesArray &lhits = *fHits;
-
- for ( hitCounter = 0 ; hitCounter < fNTmpHits && !deja ; hitCounter++ ) {
- curHit = (AliPHOSHit*) ltmphits[hitCounter] ;
- if( *curHit == *newHit ) {
- *curHit = *curHit + *newHit ;
- deja = kTRUE ;
- }
- }
-
- if ( !deja ) {
- new(ltmphits[fNTmpHits]) AliPHOSHit(*newHit) ;
- fNTmpHits++ ;
- }
-
- // We do want to save in TreeH the raw hits
- new(lhits[fNhits]) AliPHOSHit(*newHit) ;
+ new((*fHits)[fNhits]) AliPHOSHit(*newHit) ;
fNhits++ ;
- // Please note that the fTmpHits array must survive up to the
- // end of the events, so it does not appear e.g. in ResetHits() (
- // which is called at the end of each primary).
-
delete newHit;
}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff