#include "AliPHOSTrackSegment.h"
#include "AliPHOSRecParticle.h"
#include "AliPHOSCpvRecPoint.h"
-#include "AliPHOSPpsdRecPoint.h"
#include "AliPHOSGetter.h"
//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 * sdigitOccupancy = new TH2F("sdigitOccupancy","EMC sdigits", 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 ObjectGetter
AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ;
const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ;
gime->Event(Nevent);
+
+ Int_t nx = phosgeom->GetNPhi() ;
+ Int_t nz = phosgeom->GetNZ() ;
+ Float_t * cri= phosgeom->GetEMCAGeometry()->GetCrystalHalfSize() ;
+ Float_t x = nx*cri[0] ;
+ Float_t z = nz*cri[2] ;
+ Int_t nxCPV = (Int_t) (nx*phosgeom->GetPadSizePhi()/(2.*cri[0])) ;
+ Int_t nzCPV = (Int_t) (nz*phosgeom->GetPadSizeZ()/(2.*cri[2])) ;
+
+ TH2F * emcDigits = new TH2F("emcDigits","EMC digits", nx,-x,x,nz,-z,z);
+ TH2F * emcSdigits = new TH2F("emcSdigits","EMC sdigits", nx,-x,x,nz,-z,z);
+ TH2F * emcRecPoints = new TH2F("emcRecPoints","EMC RecPoints",nx,-x,x,nz,-z,z);
+ TH2F * cpvSdigits = new TH2F("cpvSdigits","CPV sdigits", nx,-x,x,nz,-z,z);
+ TH2F * cpvDigits = new TH2F("cpvDigits","CPV digits", nxCPV,-x,x,nzCPV,-z,z) ;
+ TH2F * cpvRecPoints = new TH2F("cpvCl","CPV RecPoints", nxCPV,-x,x,nzCPV,-z,z) ;
+ TH2F * nbar = new TH2F("nbar","Primary nbar", nx,-x,x,nz,-z,z);
+ TH2F * phot = new TH2F("phot","Primary Photon", nx,-x,x,nz,-z,z);
+ TH2F * charg = new TH2F("charg","Primary charged",nx,-x,x,nz,-z,z);
+ TH2F * recPhot = new TH2F("recPhot","RecParticles with primary Photon",nx,-x,x,nz,-z,z);
+ TH2F * recNbar = new TH2F("recNbar","RecParticles with primary Nbar", nx,-x,x,nz,-z,z);
+
//Plot Primary Particles
const TParticle * primary ;
{
primary = gime->Primary(iPrimary) ;
Int_t primaryType = primary->GetPdgCode() ;
- if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212) ) {
+ if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212)
+ ||(primaryType == 11)||(primaryType == -11) ) {
Int_t moduleNumber ;
Double_t primX, primZ ;
phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
}
}
}
+
Int_t iSDigit ;
AliPHOSDigit * sdigit ;
-
- TClonesArray * sdigits = gime->SDigits(branchTitle,ffileName) ;
- if(sdigits)
- for(iSDigit = 0; iSDigit < sdigits->GetEntriesFast(); iSDigit++)
+ TClonesArray * sdigits = gime->SDigits() ;
+ if(sdigits){
+ for(iSDigit = 0; iSDigit < sdigits->GetEntriesFast() ; iSDigit++)
{
- sdigit = (AliPHOSDigit*)sdigits->At(iSDigit) ;
+ sdigit = (AliPHOSDigit *) sdigits->At(iSDigit) ;
Int_t relid[4];
phosgeom->AbsToRelNumbering(sdigit->GetId(), relid) ;
Float_t x,z ;
phosgeom->RelPosInModule(relid,x,z) ;
- Float_t e = gime->SDigitizer(branchTitle)->Calibrate(sdigit->GetAmp()) ;
+ Float_t e = gime->SDigitizer()->Calibrate(sdigit->GetAmp()) ;
if(relid[0]==Nmod){
if(relid[1]==0) //EMC
- sdigitOccupancy->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) ;
+ emcSdigits->Fill(x,z,e) ;
+ if( relid[1]!=0 )
+ cpvSdigits->Fill(x,z,e) ;
}
}
-
+ }
+
//Plot digits
Int_t iDigit ;
AliPHOSDigit * digit ;
- TClonesArray * digits = gime->Digits(branchTitle) ;
- if(digits)
+ TClonesArray * digits = gime->Digits();
+ if(digits) {
for(iDigit = 0; iDigit < digits->GetEntriesFast(); iDigit++)
{
- digit = (AliPHOSDigit*) digits->At(iDigit) ;
+ digit = (AliPHOSDigit *) digits->At(iDigit) ;
Int_t relid[4];
phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
Float_t x,z ;
phosgeom->RelPosInModule(relid,x,z) ;
- Float_t e = gime->SDigitizer(branchTitle)->Calibrate(digit->GetAmp()) ;
+ Float_t e = gime->SDigitizer()->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) ;
+ emcDigits->Fill(x,z,e) ;
+ if( relid[1]!=0 )
+ cpvDigits->Fill(x,z,e) ;
}
}
+ }
//Plot RecPoints
Int_t irecp ;
TVector3 pos ;
- TObjArray * emcrp = gime->EmcRecPoints(branchTitle) ;
- TObjArray * cpvrp = gime->CpvRecPoints(branchTitle) ;
- if(cpvrp)
+ TObjArray * emcrp = gime->EmcRecPoints() ;
+ if(emcrp) {
for(irecp = 0; irecp < emcrp->GetEntriesFast() ; irecp ++){
- const AliPHOSEmcRecPoint * emc= (AliPHOSEmcRecPoint *) emcrp->At(irecp) ;
+ AliPHOSEmcRecPoint * emc = (AliPHOSEmcRecPoint *) emcrp->At(irecp) ;
if(emc->GetPHOSMod()==Nmod){
emc->GetLocalPosition(pos) ;
- emcOccupancy->Fill(pos.X(),pos.Z(),emc->GetEnergy());
+ emcRecPoints->Fill(pos.X(),pos.Z(),emc->GetEnergy());
}
}
+ }
- if(cpvrp)
+ TObjArray * cpvrp = gime->CpvRecPoints() ;
+ if(cpvrp) {
for(irecp = 0; irecp < cpvrp->GetEntriesFast() ; irecp ++){
- const AliPHOSRecPoint * cpv = (AliPHOSRecPoint *) cpvrp->At(irecp) ;
- if((strcmp(cpv->ClassName(),"AliPHOSPpsdRecPoint" )) == 0){ // PPSD Rec Point
- AliPHOSPpsdRecPoint * ppsd = (AliPHOSPpsdRecPoint*) cpv ;
- ppsd->GetLocalPosition(pos) ;
- 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());
- }
- }
- else{
- AliPHOSCpvRecPoint * cpv1 = (AliPHOSCpvRecPoint*) cpv ;
- if(cpv1->GetPHOSMod()==Nmod){
- cpv1->GetLocalPosition(pos) ;
- ppsdUpCl->Fill(pos.X(),pos.Z(),cpv1->GetEnergy());
- }
+ AliPHOSRecPoint * cpv = (AliPHOSCpvRecPoint *) cpvrp->At(irecp) ;
+ if(cpv->GetPHOSMod()==Nmod){
+ cpv->GetLocalPosition(pos) ;
+ cpvRecPoints->Fill(pos.X(),pos.Z(),cpv->GetEnergy());
}
}
-
-
+ }
+
//Plot RecParticles
- const AliPHOSRecParticle * recParticle ;
+ AliPHOSRecParticle * recParticle ;
Int_t iRecParticle ;
- TClonesArray * rps = gime->RecParticles(branchTitle) ;
- TClonesArray * tss = gime->TrackSegments(branchTitle) ;
- if(rps && tss && emcrp)
- for(iRecParticle = 0; iRecParticle < rps->GetEntriesFast() ;iRecParticle++ )
+ TClonesArray * rp = gime->RecParticles() ;
+ TClonesArray * ts = gime->TrackSegments() ;
+ if(rp && ts && emcrp) {
+ for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ; iRecParticle++ )
{
- recParticle = (AliPHOSRecParticle *) rps->At(iRecParticle) ;
+ recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
Int_t moduleNumberRec ;
Double_t recX, recZ ;
phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
if(moduleNumberRec == Nmod){
Double_t minDistance = 5. ;
- Int_t closestPrimary = -1 ;
-
-
+ Int_t closestPrimary = -1 ;
+
//extract list of primaries: it is stored at EMC RecPoints
- Int_t emcIndex = ((AliPHOSTrackSegment *) tss->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
+ Int_t emcIndex = ((AliPHOSTrackSegment *) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
Int_t numberofprimaries ;
- Int_t * listofprimaries = ((AliPHOSEmcRecPoint *)emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ;
+ Int_t * listofprimaries = ((AliPHOSRecPoint*) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ;
Int_t index ;
const TParticle * primary ;
Double_t distance = minDistance ;
}
}
}
-
+
+ }
//Plot made histograms
- digitOccupancy->Draw("box") ;
- sdigitOccupancy->SetLineColor(5) ;
- sdigitOccupancy->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") ;
+ emcSdigits->Draw("box") ;
+ emcDigits->SetLineColor(5) ;
+ emcDigits->Draw("box") ;
+ emcRecPoints->SetLineColor(2) ;
+ emcRecPoints->Draw("boxsame") ;
+ cpvSdigits->SetLineColor(1) ;
+ cpvSdigits->Draw("box") ;
+ charg->SetLineColor(2) ;
+ charg->Draw("boxsame") ;
nbar->SetLineColor(6) ;
nbar->Draw("boxsame") ;
Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
// for(event = 0; event < gime->MaxEvent(); event++ ){
for(event = 0; event < maxevent; event++ ){
- gime->Event(event);
+ gime->Event(event,"R"); //will read only TreeR
//copy EM RecParticles to the "total" list
const AliPHOSRecParticle * recParticle ;
Int_t iRecParticle ;
- TClonesArray * rp = gime->RecParticles(branchTitle) ;
- if(rp)
- for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ;iRecParticle++ )
- {
- recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
- if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
- (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM))
- new( (*allRecParticleList)[iRecPhot++] ) AliPHOSRecParticle(*recParticle) ;
- }
+ TClonesArray * rp = gime->RecParticles() ;
+ if(!rp){
+ cout << "AliPHOSAnalyze::InvariantMass --> Can't find RecParticles " << endl ;
+ return ;
+ }
+
+ for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ )
+ {
+ recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
+ if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)||
+ (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW))
+ new( (*allRecParticleList)[iRecPhot++] ) AliPHOSRecParticle(*recParticle) ;
+ }
Int_t mevent = event%nMixedEvents ; //event number in the "mixed" cicle
nRecParticles[mevent] = iRecPhot-1 ;
-
+
//check, if it is time to calculate invariant mass?
Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
if((mevent == 0) && (event +1 == maxevent)){
-
- // if((mevent == 0) && (event +1 == gime->MaxEvent())){
+
+ // if((mevent == 0) && (event +1 == gime->MaxEvent())){
//calculate invariant mass:
Int_t irp1,irp2 ;
if(irp2 <= nRecParticles[nCurEvent]){ //'Real' event
hRealEM->Fill(invMass,pt);
- if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&
- (rp2->GetType() == AliPHOSFastRecParticle::kGAMMA) )
+ if((rp1->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)&&
+ (rp2->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST) )
hRealPhot->Fill(invMass,pt);
}
else{
hMixedEM->Fill(invMass,pt);
- if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&
- (rp2->GetType() == AliPHOSFastRecParticle::kGAMMA) )
+ if((rp1->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)&&
+ (rp2->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST) )
hMixedPhot->Fill(invMass,pt);
} //real-mixed
Int_t ievent;
Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
- // for ( ievent=0; ievent < gime->MaxEvent() ; ievent++){
for ( ievent=0; ievent < maxevent ; ievent++){
//read the current event
const AliPHOSRecParticle * recParticle ;
Int_t iRecParticle ;
- TClonesArray * rp = gime->RecParticles(branchTitle) ;
- TClonesArray * tss = gime->TrackSegments(branchTitle) ;
- TObjArray * emcrp = gime->EmcRecPoints(branchTitle) ;
- if(rp && tss && emcrp)
- for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ; iRecParticle++ ){
- recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
-
- //find the closest primary
- Int_t moduleNumberRec ;
- Double_t recX, recZ ;
- phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
-
- Double_t minDistance = 100. ;
- Int_t closestPrimary = -1 ;
-
- //extract list of primaries: it is stored at EMC RecPoints
- Int_t emcIndex = ((AliPHOSTrackSegment *)tss->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
- Int_t numberofprimaries ;
- Int_t * listofprimaries = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ;
-
- Int_t index ;
- const TParticle * primary ;
- Double_t distance = minDistance ;
- Double_t dX, dZ;
- Double_t dXmin = 0.;
- Double_t dZmin = 0. ;
- for ( index = 0 ; index < numberofprimaries ; index++){
- primary = gime->Primary(listofprimaries[index]) ;
- Int_t moduleNumber ;
- Double_t primX, primZ ;
- phosgeom->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] ;
- }
+ TClonesArray * rp = gime->RecParticles() ;
+ if(!rp) {
+ cout << "AliPHOSAnalyze::EnergyResolution --> Event " <<ievent
+ << ", Can't find RecParticles " << endl ;
+ return ;
+ }
+ TClonesArray * ts = gime->TrackSegments() ;
+ if(!ts) {
+ cout << "AliPHOSAnalyze::EnergyResolution --> Event " <<ievent
+ << ", Can't find TrackSegments " << endl ;
+ return ;
+ }
+ TObjArray * emcrp = gime->EmcRecPoints() ;
+ if(!emcrp){
+ cout << "AliPHOSAnalyze::EnergyResolution --> Event " <<ievent
+ << ", Can't find EmcRecPoints " << endl ;
+ return ;
+ }
+
+ for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ;iRecParticle++ ){
+ recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
+
+ //find the closest primary
+ Int_t moduleNumberRec ;
+ Double_t recX, recZ ;
+ phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
+
+ Double_t minDistance = 100. ;
+ Int_t closestPrimary = -1 ;
+
+ //extract list of primaries: it is stored at EMC RecPoints
+ Int_t emcIndex = ((AliPHOSTrackSegment*) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
+ Int_t numberofprimaries ;
+ Int_t * listofprimaries = ((AliPHOSEmcRecPoint*) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ;
+
+ Int_t index ;
+ const TParticle * primary ;
+ Double_t distance = minDistance ;
+ Double_t dX, dZ;
+ Double_t dXmin = 0.;
+ Double_t dZmin = 0. ;
+ for ( index = 0 ; index < numberofprimaries ; index++){
+ primary = gime->Primary(listofprimaries[index]) ;
+ Int_t moduleNumber ;
+ Double_t primX, primZ ;
+ phosgeom->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] ;
}
}
-
- //if found primary, fill histograms
- if(closestPrimary >=0 ){
- const TParticle * primary = gime->Primary(closestPrimary) ;
- if(primary->GetPdgCode() == 22){
- hAllEnergy->Fill(primary->Energy(), recParticle->Energy()) ;
- if(recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA){
- hPhotEnergy->Fill(primary->Energy(), recParticle->Energy() ) ;
- hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ;
- }
- else
- if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)
- hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ;
+ }
+
+ //if found primary, fill histograms
+ if(closestPrimary >=0 ){
+ const TParticle * primary = gime->Primary(closestPrimary) ;
+ if(primary->GetPdgCode() == 22){
+ hAllEnergy->Fill(primary->Energy(), recParticle->Energy()) ;
+ if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){
+ hPhotEnergy->Fill(primary->Energy(), recParticle->Energy() ) ;
+ hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ;
}
+ else
+ if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)
+ hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ;
}
}
+ }
}
-
+
//write filled histograms
efile->cd() ;
hAllEnergy->Write(0,kOverwrite) ;
// efile->Write() ;
efile->Close() ;
delete efile ;
-
+
}
//____________________________________________________________________________
void AliPHOSAnalyze::PositionResolution(const char * branchTitle)
const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ;
Int_t ievent;
- Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
- // for ( ievent=0; ievent < gime->MaxEvent() ; ievent++){
+ Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
for ( ievent=0; ievent < maxevent ; ievent++){
-
+
//read the current event
gime->Event(ievent) ;
+ TClonesArray * rp = gime->RecParticles() ;
+ if(!rp) {
+ cout << "AliPHOSAnalyze::PositionResolution --> Event " <<ievent
+ << ", Can't find RecParticles " << endl ;
+ return ;
+ }
+ TClonesArray * ts = gime->TrackSegments() ;
+ if(!ts) {
+ cout << "AliPHOSAnalyze::PositionResolution --> Event " <<ievent
+ << ", Can't find TrackSegments " << endl ;
+ return ;
+ }
+ TObjArray * emcrp = gime->EmcRecPoints() ;
+ if(!emcrp){
+ cout << "AliPHOSAnalyze::PositionResolution --> Event " <<ievent
+ << ", Can't find EmcRecPoints " << endl ;
+ return ;
+ }
+
const AliPHOSRecParticle * recParticle ;
Int_t iRecParticle ;
- TClonesArray * rp = gime->RecParticles(branchTitle) ;
- TClonesArray * tss= gime->TrackSegments(branchTitle) ;
- TObjArray * emcrp = gime->EmcRecPoints(branchTitle) ;
- if( rp && tss && emcrp )
- for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ;iRecParticle++ ){
- recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
-
- //find the closest primary
- Int_t moduleNumberRec ;
- Double_t recX, recZ ;
- phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
-
- Double_t minDistance = 100. ;
- Int_t closestPrimary = -1 ;
-
- //extract list of primaries: it is stored at EMC RecPoints
- Int_t emcIndex = ((AliPHOSTrackSegment *) tss->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
- Int_t numberofprimaries ;
- Int_t * listofprimaries = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ;
-
- Int_t index ;
- const TParticle * primary ;
- Double_t distance = minDistance ;
- Double_t dX = 1000; // incredible number
- Double_t dZ = 1000; // for the case if no primary will be found
- Double_t dXmin = 0.;
- Double_t dZmin = 0. ;
- for ( index = 0 ; index < numberofprimaries ; index++){
- primary = gime->Primary(listofprimaries[index]) ;
- Int_t moduleNumber ;
- Double_t primX, primZ ;
- phosgeom->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] ;
- }
+ for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ ){
+ recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
+
+ //find the closest primary
+ Int_t moduleNumberRec ;
+ Double_t recX, recZ ;
+ phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
+
+ Double_t minDistance = 100. ;
+ Int_t closestPrimary = -1 ;
+
+ //extract list of primaries: it is stored at EMC RecPoints
+ Int_t emcIndex = ((AliPHOSTrackSegment*) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
+ Int_t numberofprimaries ;
+ Int_t * listofprimaries = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ;
+
+ Int_t index ;
+ const TParticle * primary ;
+ Double_t distance = minDistance ;
+ Double_t dX = 1000; // incredible number
+ Double_t dZ = 1000; // for the case if no primary will be found
+ Double_t dXmin = 0.;
+ Double_t dZmin = 0. ;
+ for ( index = 0 ; index < numberofprimaries ; index++){
+ primary = gime->Primary(listofprimaries[index]) ;
+ Int_t moduleNumber ;
+ Double_t primX, primZ ;
+ phosgeom->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] ;
}
}
-
- //if found primary, fill histograms
- if(closestPrimary >=0 ){
- const TParticle * primary = gime->Primary(closestPrimary) ;
- if(primary->GetPdgCode() == 22){
- hAllPosition->Fill(primary->Energy(), minDistance) ;
- hAllPositionX->Fill(primary->Energy(), dX) ;
- hAllPositionZ->Fill(primary->Energy(), dZ) ;
- if(recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA){
- hPhotPosition->Fill(primary->Energy(), minDistance ) ;
- hEMPosition->Fill(primary->Energy(), minDistance ) ;
- }
- else
- if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)
- hEMPosition->Fill(primary->Energy(), minDistance ) ;
+ }
+
+ //if found primary, fill histograms
+ if(closestPrimary >=0 ){
+ const TParticle * primary = gime->Primary(closestPrimary) ;
+ if(primary->GetPdgCode() == 22){
+ hAllPosition->Fill(primary->Energy(), minDistance) ;
+ hAllPositionX->Fill(primary->Energy(), dX) ;
+ hAllPositionZ->Fill(primary->Energy(), dZ) ;
+ if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){
+ hPhotPosition->Fill(primary->Energy(), minDistance ) ;
+ hEMPosition->Fill(primary->Energy(), minDistance ) ;
}
+ else
+ if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)
+ hEMPosition->Fill(primary->Energy(), minDistance ) ;
}
}
+ }
}
//Write output histgrams
pfile->Write() ;
pfile->Close() ;
delete pfile ;
-
-
+
+
}
//____________________________________________________________________________
void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){
TH1F * hPrimary = 0; //spectrum (P_t distribution) of primary photons
TH1F * hAllRP = 0; //spectrum of all RecParticles in PHOS
TH1F * hPhot = 0; //spectrum of kGAMMA RecParticles
- TH1F * hPPSD = 0; //spectrum of all RecParticles with PPSD signal
TH1F * hShape = 0; //spectrum of all EM RecParticles
TH1F * hVeto = 0; //spectrum of all neutral RecParticles
//primary - photon
TH1F * hPhotReg = 0; //Registeres as photon
TH1F * hPhotEM = 0; //Registered as EM
- TH1F * hPhotPPSD= 0; //Registered as RecParticle with PPSD signal
//primary - n
TH1F * hNReg = 0; //Registeres as photon
TH1F * hNEM = 0; //Registered as EM
- TH1F * hNPPSD= 0; //Registered as RecParticle with PPSD signal
//primary - nBar
TH1F * hNBarReg = 0; //Registeres as photon
TH1F * hNBarEM = 0; //Registered as EM
- TH1F * hNBarPPSD= 0; //Registered as RecParticle with PPSD signal
//primary - charged hadron (pBar excluded)
TH1F * hChargedReg = 0; //Registeres as photon
TH1F * hChargedEM = 0; //Registered as EM
- TH1F * hChargedPPSD= 0; //Registered as RecParticle with PPSD signal
//primary - pBar
TH1F * hPbarReg = 0; //Registeres as photon
TH1F * hPbarEM = 0; //Registered as EM
- TH1F * hPbarPPSD= 0; //Registered as RecParticle with PPSD signal
//Reading histograms from the file
hPhot = (TH1F*)cfile->Get("hPhot") ;
if(hPhot == 0)
hPhot = new TH1F("hPhot","All kGAMMA RecParticles",100, 0., 5.);
- hPPSD = (TH1F*)cfile->Get("hPPSD") ;
- if(hPPSD == 0)
- hPPSD = new TH1F("hPPSD", "All PPSD Recparticles", 100, 0., 5.);
hShape = (TH1F*) cfile->Get("hShape") ;
if(hShape == 0)
hShape = new TH1F("hShape","All particles with EM shower",100, 0., 5.);
hPhotEM =(TH1F*)cfile->Get("hPhotEM");
if(hPhotEM== 0)
hPhotEM = new TH1F("hPhotEM", "Photon registered as EM", 100, 0., 5.);
- hPhotPPSD= (TH1F*)cfile->Get("hPhotPPSD");
- if(hPhotPPSD== 0)
- hPhotPPSD = new TH1F("hPhotPPSD","Photon registered as PPSD", 100, 0., 5.);
//primary - n
hNReg = (TH1F*)cfile->Get("hNReg");
hNEM = (TH1F*)cfile->Get("hNEM");
if(hNEM== 0)
hNEM = new TH1F("hNEM", "N registered as EM", 100, 0., 5.);
- hNPPSD=(TH1F*)cfile->Get("hNPPSD");
- if(hNPPSD== 0)
- hNPPSD = new TH1F("hNPPSD","N registered as PPSD", 100, 0., 5.);
//primary - nBar
hNBarReg =(TH1F*)cfile->Get("hNBarReg");
hNBarEM =(TH1F*)cfile->Get("hNBarEM");
if(hNBarEM== 0)
hNBarEM = new TH1F("hNBarEM", "NBar registered as EM", 100, 0., 5.);
- hNBarPPSD=(TH1F*)cfile->Get("hNBarPPSD");
- if(hNBarPPSD== 0)
- hNBarPPSD = new TH1F("hNBarPPSD","NBar registered as PPSD", 100, 0., 5.);
//primary - charged hadron (pBar excluded)
hChargedReg = (TH1F*)cfile->Get("hChargedReg");
hChargedEM = (TH1F*)cfile->Get("hChargedEM");
if(hChargedEM== 0)
hChargedEM= new TH1F("hChargedEM","Charged registered as EM",100, 0., 5.);
- hChargedPPSD= (TH1F*)cfile->Get("hChargedPPSD");
- if(hChargedPPSD== 0)
- hChargedPPSD= new TH1F("hChargedPPSD","Charged registered as PPSD",100, 0., 5.);
-
+
//primary - pBar
hPbarReg = (TH1F*)cfile->Get("hPbarReg");
if(hPbarReg== 0)
hPbarEM = (TH1F*)cfile->Get("hPbarEM");
if(hPbarEM== 0)
hPbarEM= new TH1F("hPbarEM","Pbar registered as EM",100, 0., 5.);
- hPbarPPSD= (TH1F*)cfile->Get("hPbarPPSD");
- if(hPbarPPSD== 0)
- hPbarPPSD= new TH1F("hPbarPPSD","Pbar as PPSD",100, 0., 5.);
//Now make some initializations
Int_t ievent;
Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
- // for ( ievent=0; ievent < gime->MaxEvent() ; ievent++){
for ( ievent=0; ievent < maxevent ; ievent++){
gime->Event(ievent) ;
+ TClonesArray * rp = gime->RecParticles() ;
+ if(!rp) {
+ cout << "AliPHOSAnalyze::Contamination --> Event " <<ievent
+ << ", Can't find RecParticles " << endl ;
+ return ;
+ }
+ TClonesArray * ts = gime->TrackSegments() ;
+ if(!ts) {
+ cout << "AliPHOSAnalyze::Contamination --> Event " <<ievent
+ << ", Can't find TrackSegments " << endl ;
+ return ;
+ }
+ TObjArray * emcrp = gime->EmcRecPoints() ;
+ if(!emcrp){
+ cout << "AliPHOSAnalyze::Contamination --> Event " <<ievent
+ << ", Can't find EmcRecPoints " << endl ;
+ return ;
+ }
+
+
//=========== Make spectrum of the primary photons
const TParticle * primary ;
Int_t iPrimary ;
}
}
+
//========== Now scan over RecParticles
const AliPHOSRecParticle * recParticle ;
Int_t iRecParticle ;
- TClonesArray * rp = gime->RecParticles(RecPointsTitle) ;
- TClonesArray * tss= gime->TrackSegments(RecPointsTitle) ;
- TObjArray * emcrp = gime->EmcRecPoints(RecPointsTitle) ;
- if( rp && tss && emcrp )
- for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ ){
- recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
- //fill histo spectrum of all RecParticles
- hAllRP->Fill(CorrectedEnergy(recParticle->Energy())) ;
-
- //==========find the closest primary
- Int_t moduleNumberRec ;
- Double_t recX, recZ ;
- phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
-
- Double_t minDistance = 100. ;
- Int_t closestPrimary = -1 ;
-
- //extract list of primaries: it is stored at EMC RecPoints
- Int_t emcIndex = ((AliPHOSTrackSegment *) tss->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
- Int_t numberofprimaries ;
- Int_t * listofprimaries = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ;
- Int_t index ;
- const TParticle * primary ;
- Double_t distance = minDistance ;
- Double_t dX, dZ;
- Double_t dXmin = 0.;
- Double_t dZmin = 0. ;
- for ( index = 0 ; index < numberofprimaries ; index++){
- primary = gime->Primary(listofprimaries[index]) ;
- Int_t moduleNumber ;
- Double_t primX, primZ ;
- phosgeom->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] ;
- }
+ for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ ){
+ recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
+ //fill histo spectrum of all RecParticles
+ hAllRP->Fill(CorrectedEnergy(recParticle->Energy())) ;
+
+ //==========find the closest primary
+ Int_t moduleNumberRec ;
+ Double_t recX, recZ ;
+ phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
+
+ Double_t minDistance = 100. ;
+ Int_t closestPrimary = -1 ;
+
+ //extract list of primaries: it is stored at EMC RecPoints
+ Int_t emcIndex = ((AliPHOSTrackSegment *) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
+ Int_t numberofprimaries ;
+ Int_t * listofprimaries = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ;
+ Int_t index ;
+ const TParticle * primary ;
+ Double_t distance = minDistance ;
+ Double_t dX, dZ;
+ Double_t dXmin = 0.;
+ Double_t dZmin = 0. ;
+ for ( index = 0 ; index < numberofprimaries ; index++){
+ primary = gime->Primary(listofprimaries[index]) ;
+ Int_t moduleNumber ;
+ Double_t primX, primZ ;
+ phosgeom->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] ;
}
}
-
- //===========define the "type" of closest primary
- if(closestPrimary >=0 ){
- Int_t primaryCode = -1;
- const TParticle * primary = gime->Primary(closestPrimary) ;
- Int_t primaryType = primary->GetPdgCode() ;
- if(primaryType == 22) // photon ?
- primaryCode = 0 ;
+ }
+
+ //===========define the "type" of closest primary
+ if(closestPrimary >=0 ){
+ Int_t primaryCode = -1;
+ const TParticle * primary = gime->Primary(closestPrimary) ;
+ Int_t primaryType = primary->GetPdgCode() ;
+ if(primaryType == 22) // photon ?
+ primaryCode = 0 ;
+ else
+ if(primaryType == 2112) // neutron
+ primaryCode = 1 ;
else
- if(primaryType == 2112) // neutron
- primaryCode = 1 ;
+ if(primaryType == -2112) // Anti neutron
+ primaryCode = 2 ;
else
- if(primaryType == -2112) // Anti neutron
- primaryCode = 2 ;
- else
- if(primaryType == -2122) //Anti proton
- primaryCode = 4 ;
- else {
- TParticle tempo(*primary) ;
- if(tempo.GetPDG()->Charge())
- primaryCode = 3 ;
- }
-
- //==========Now look at the type of RecParticle
- Float_t energy = CorrectedEnergy(recParticle->Energy()) ;
- if(recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA){
- hPhot->Fill(energy ) ;
- switch(primaryCode){
- case 0:
- hPhotReg->Fill(energy ) ;
- break ;
- case 1:
- hNReg->Fill(energy ) ;
- break ;
- case 2:
- hNBarReg->Fill(energy ) ;
- break ;
- case 3:
- hChargedReg->Fill(energy ) ;
- break ;
- case 4:
- hPbarReg->Fill(energy ) ;
- break ;
- default:
- break ;
- }
- }
- if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
- (recParticle->GetType() == AliPHOSFastRecParticle::kGAMMAHA)){ //with PPSD signal
- hPPSD->Fill(energy ) ;
- switch(primaryCode){
- case 0:
- hPhotPPSD->Fill(energy ) ;
- break ;
- case 1:
- hNPPSD->Fill(energy ) ;
- break ;
- case 2:
- hNBarPPSD->Fill(energy ) ;
- break ;
- case 3:
- hChargedPPSD->Fill(energy ) ;
- break ;
- case 4:
- hPbarPPSD->Fill(energy ) ;
- break ;
- default:
- break ;
- }
+ if(primaryType == -2122) //Anti proton
+ primaryCode = 4 ;
+ else {
+ TParticle tempo(*primary) ;
+ if(tempo.GetPDG()->Charge())
+ primaryCode = 3 ;
+ }
+
+ //==========Now look at the type of RecParticle
+ Float_t energy = CorrectedEnergy(recParticle->Energy()) ;
+ if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){
+ hPhot->Fill(energy ) ;
+ switch(primaryCode){
+ case 0:
+ hPhotReg->Fill(energy ) ;
+ break ;
+ case 1:
+ hNReg->Fill(energy ) ;
+ break ;
+ case 2:
+ hNBarReg->Fill(energy ) ;
+ break ;
+ case 3:
+ hChargedReg->Fill(energy ) ;
+ break ;
+ case 4:
+ hPbarReg->Fill(energy ) ;
+ break ;
+ default:
+ break ;
}
- if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
- (recParticle->GetType() == AliPHOSFastRecParticle::kELECTRON)||
- (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)||
- (recParticle->GetType() == AliPHOSFastRecParticle::kABSURDEM) ){ //with EM shower
- hShape->Fill(energy ) ;
- switch(primaryCode){
- case 0:
- hPhotEM->Fill(energy ) ;
- break ;
- case 1:
- hNEM->Fill(energy ) ;
- break ;
- case 2:
- hNBarEM->Fill(energy ) ;
- break ;
- case 3:
- hChargedEM->Fill(energy ) ;
- break ;
- case 4:
- hPbarEM->Fill(energy ) ;
- break ;
- default:
- break ;
- }
+ }
+ if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)||
+ (recParticle->GetType() == AliPHOSFastRecParticle::kCHARGEDEMFAST)||
+ (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)||
+ (recParticle->GetType() == AliPHOSFastRecParticle::kCHARGEDEMSLOW) ){ //with EM shower
+ hShape->Fill(energy ) ;
+ switch(primaryCode){
+ case 0:
+ hPhotEM->Fill(energy ) ;
+ break ;
+ case 1:
+ hNEM->Fill(energy ) ;
+ break ;
+ case 2:
+ hNBarEM->Fill(energy ) ;
+ break ;
+ case 3:
+ hChargedEM->Fill(energy ) ;
+ break ;
+ case 4:
+ hPbarEM->Fill(energy ) ;
+ break ;
+ default:
+ break ;
}
-
- if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
- (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHA) ||
- (recParticle->GetType() == AliPHOSFastRecParticle::kGAMMAHA) ||
- (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM) ) //nuetral
- hVeto->Fill(energy ) ;
-
- //fill number of primaries identified as ...
- if(primaryCode >= 0) // Primary code defined
- counter[recParticle->GetType()][primaryCode]++ ;
-
}
- } // no closest primary found
+ if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)||
+ (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHAFAST) ||
+ (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW) ||
+ (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHASLOW) ) //nuetral
+ hVeto->Fill(energy ) ;
+
+ //fill number of primaries identified as ...
+ if(primaryCode >= 0) // Primary code defined
+ counter[recParticle->GetType()][primaryCode]++ ;
+
+ }
+
+ } // no closest primary found
}
hPrimary->Write(0,kOverwrite);
hAllRP->Write(0,kOverwrite);
hPhot->Write(0,kOverwrite);
- hPPSD->Write(0,kOverwrite);
hShape->Write(0,kOverwrite);
hVeto->Write(0,kOverwrite);
hPhotReg->Write(0,kOverwrite);
hPhotEM->Write(0,kOverwrite);
- hPhotPPSD->Write(0,kOverwrite);
hNReg ->Write(0,kOverwrite);
hNEM ->Write(0,kOverwrite);
- hNPPSD->Write(0,kOverwrite);
hNBarReg ->Write(0,kOverwrite);
hNBarEM ->Write(0,kOverwrite);
- hNBarPPSD->Write(0,kOverwrite);
hChargedReg ->Write(0,kOverwrite);
hChargedEM ->Write(0,kOverwrite);
- hChargedPPSD->Write(0,kOverwrite);
hPbarReg ->Write(0,kOverwrite);
hPbarEM ->Write(0,kOverwrite);
- hPbarPPSD->Write(0,kOverwrite);
cfile->Write(0,kOverwrite);
cfile->Close();
virtual Float_t GetEmcClusteringThreshold()const = 0 ;
virtual Float_t GetEmcLocalMaxCut()const = 0 ;
virtual Float_t GetEmcLogWeight()const = 0 ;
+ virtual Float_t GetEmcTimeGate() const = 0 ;
virtual Float_t GetCpvClusteringThreshold()const = 0 ;
virtual Float_t GetCpvLocalMaxCut()const = 0 ;
virtual Float_t GetCpvLogWeight()const = 0 ;
- virtual Float_t GetPpsdClusteringThreshold()const = 0 ;
virtual char * GetRecPointsBranch() const = 0 ;
virtual const Int_t GetRecPointsInRun() const = 0 ;
virtual char * GetDigitsBranch() const = 0 ;
virtual void SetEmcClusteringThreshold(Float_t cluth) = 0 ;
virtual void SetEmcLocalMaxCut(Float_t cut) = 0 ;
virtual void SetEmcLogWeight(Float_t w) = 0 ;
+ virtual void SetEmcTimeGate(Float_t gate) = 0 ;
virtual void SetCpvClusteringThreshold(Float_t cluth) = 0 ;
virtual void SetCpvLocalMaxCut(Float_t cut) = 0 ;
virtual void SetCpvLogWeight(Float_t w) = 0 ;
- virtual void SetPpsdClusteringThreshold(Float_t cluth) = 0 ;
virtual void SetDigitsBranch(const char * title) = 0 ;
virtual void SetRecPointsBranch(const char *title) = 0 ;
virtual void SetUnfolding(Bool_t toUnfold ) = 0 ;
#include "AliPHOSDigitizer.h"
#include "AliPHOSEmcRecPoint.h"
#include "AliPHOS.h"
-#include "AliPHOSPpsdRecPoint.h"
#include "AliPHOSGetter.h"
#include "AliRun.h"
fCpvClusteringThreshold = 0.0;
fEmcClusteringThreshold = 0.2;
- fPpsdClusteringThreshold = 0.0000002 ;
fEmcLocMaxCut = 0.03 ;
fCpvLocMaxCut = 0.03 ;
fW0 = 4.5 ;
fW0CPV = 4.0 ;
+ fEmcTimeGate = 1.e-8 ;
+
+ fToUnfold = kTRUE ;
+
fHeaderFileName = "" ;
fDigitsBranchTitle = "" ;
- fRecPointsInRun = 0 ;
+ fRecPointsInRun = 0 ;
}
//____________________________________________________________________________
fCpvClusteringThreshold = 0.0;
fEmcClusteringThreshold = 0.2;
- fPpsdClusteringThreshold = 0.0000002 ;
fEmcLocMaxCut = 0.03 ;
fCpvLocMaxCut = 0.03 ;
fW0 = 4.5 ;
fW0CPV = 4.0 ;
+
+ fEmcTimeGate = 1.e-8 ;
fToUnfold = kTRUE ;
fHeaderFileName = GetTitle() ;
fDigitsBranchTitle = GetName() ;
-
- TString tempo(GetName()) ;
- tempo.Append(":") ;
- tempo.Append(Version()) ;
- SetName(tempo.Data()) ;
+
+ TString clusterizerName( GetName()) ;
+ clusterizerName.Append(":") ;
+ clusterizerName.Append(Version()) ;
+ SetName(clusterizerName) ;
fRecPointsInRun = 0 ;
Init() ;
}
-
+//____________________________________________________________________________
+ AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
+{
+}
//____________________________________________________________________________
void AliPHOSClusterizerv1::Exec(Option_t * option)
{
if(strstr(option,"print"))
Print("") ;
- //check, if the branch with name of this" already exits?
gAlice->GetEvent(0) ;
+
+ //check, if the branch with name of this" already exits?
TObjArray * lob = (TObjArray*)gAlice->TreeR()->GetListOfBranches() ;
TIter next(lob) ;
TBranch * branch = 0 ;
Bool_t phosemcfound = kFALSE, phoscpvfound = kFALSE, clusterizerfound = kFALSE ;
-
- TString taskName(GetName()) ;
- taskName.Remove(taskName.Index(Version()) -1) ;
+ TString branchname = GetName() ;
+ branchname.Remove(branchname.Index(Version())-1) ;
+
while ( (branch = (TBranch*)next()) && (!phosemcfound || !phoscpvfound || !clusterizerfound) ) {
- if ( (strcmp(branch->GetName(), "PHOSEmcRP")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) )
+ if ( (strcmp(branch->GetName(), "PHOSEmcRP")==0) && (strcmp(branch->GetTitle(), branchname.Data())==0) )
phosemcfound = kTRUE ;
- else if ( (strcmp(branch->GetName(), "PHOSCpvRP")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) )
+ else if ( (strcmp(branch->GetName(), "PHOSCpvRP")==0) && (strcmp(branch->GetTitle(), branchname.Data())==0) )
phoscpvfound = kTRUE ;
- else if ( (strcmp(branch->GetName(), "AliPHOSClusterizer")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) )
+ else if ((strcmp(branch->GetName(), "AliPHOSClusterizer")==0) && (strcmp(branch->GetTitle(), GetName())==0) )
clusterizerfound = kTRUE ;
}
if ( phoscpvfound || phosemcfound || clusterizerfound ) {
cerr << "WARNING: AliPHOSClusterizer::Exec -> Emc(Cpv)RecPoints and/or Clusterizer branch with name "
- << taskName.Data() << " already exits" << endl ;
+ << branchname.Data() << " already exits" << endl ;
return ;
}
for(ievent = 0; ievent < nevents; ievent++){
- fPedestal = gime->Digitizer(taskName)->GetPedestal() ;
- fSlope = gime->Digitizer(taskName)->GetSlope() ;
+ fPedestal = gime->Digitizer(branchname)->GetPedestal() ;
+ fSlope = gime->Digitizer(branchname)->GetSlope() ;
fNumberOfEmcClusters = 0 ;
fNumberOfCpvClusters = 0 ;
gime->Event(ievent,"D") ;
- // if(!ReadDigits(ievent)) //reads digits for event fEvent
- // continue;
-
+
+ // if(!ReadDigits(ievent)) continue; //reads digits for event ievent
+
MakeClusters() ;
- if(fToUnfold)
- MakeUnfolding() ;
+ if(fToUnfold) MakeUnfolding() ;
WriteRecPoints(ievent) ;
- if(strstr(option,"deb"))
- PrintRecPoints(option) ;
+ if(strstr(option,"deb")) PrintRecPoints(option) ;
}
if(strstr(option,"tim")){
if ( strcmp(GetTitle(), "") == 0 )
SetTitle("galice.root") ;
- TString taskName(GetName()) ;
- taskName.Remove(taskName.Index(Version()) -1) ;
+ TString branchname = GetName() ;
+ branchname.Remove(branchname.Index(Version())-1) ;
- AliPHOSGetter * gime = AliPHOSGetter::GetInstance(GetTitle(), taskName.Data()) ;
+ AliPHOSGetter * gime = AliPHOSGetter::GetInstance(GetTitle(), branchname) ;
if ( gime == 0 ) {
cerr << "ERROR: AliPHOSClusterizerv1::Init -> Could not obtain the Getter object !" << endl ;
return ;
if(!gMinuit)
gMinuit = new TMinuit(100) ;
- //add Task to //YSAlice/tasks/Reconstructioner/PHOS
gime->PostClusterizer(this) ;
// create a folder on the white board //YSAlice/WhiteBoard/RecPoints/PHOS/recpointsName
- gime->PostRecPoints(taskName.Data() ) ;
+ gime->PostRecPoints(branchname ) ;
- gime->PostDigits(taskName.Data()) ;
- gime->PostDigitizer(taskName.Data()) ;
+ gime->PostDigits(branchname) ;
+ gime->PostDigitizer(branchname) ;
}
Int_t relid2[4] ;
geom->AbsToRelNumbering(d2->GetId(), relid2) ;
- if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module and the same PPSD Module
+ if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
if (( coldiff <= 1 ) && ( rowdiff <= 1 )){
+ if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate))
rv = 1 ;
}
else {
}
else {
- if( (relid1[0] < relid2[0]) || (relid1[1] < relid2[1]) )
+ if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )
rv=2 ;
}
- //Do NOT clusterize upper PPSD
- if( IsInPpsd(d1) && IsInPpsd(d2) &&
- relid1[1] > 0 &&
- relid1[1] < geom->GetNumberOfPadsPhi()*geom->GetNumberOfPadsPhi() ) rv = 2 ;
-
return rv ;
}
return rv ;
}
-//____________________________________________________________________________
-Bool_t AliPHOSClusterizerv1::IsInPpsd(AliPHOSDigit * digit) const
-{
- // Tells if (true) or not (false) the digit is in a PHOS-PPSD module
-
- Bool_t rv = kFALSE ;
- const AliPHOSGeometry * geom = AliPHOSGetter::GetInstance()->PHOSGeometry() ;
-
- Int_t relid[4] ;
- geom->AbsToRelNumbering(digit->GetId(), relid) ;
-
- if ( relid[1] > 0 && relid[0] > geom->GetNCPVModules() ) rv = kTRUE;
-
- return rv ;
-}
-
//____________________________________________________________________________
Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
{
Int_t relid[4] ;
geom->AbsToRelNumbering(digit->GetId(), relid) ;
- if ( relid[1] > 0 && relid[0] <= geom->GetNCPVModules() ) rv = kTRUE;
+ if ( relid[1] != 0 ) rv = kTRUE;
return rv ;
}
-//____________________________________________________________________________
-Bool_t AliPHOSClusterizerv1::ReadDigits(Int_t event)
- {
- // reads digitis with specified title from TreeD
-
-
- if ( gAlice->TreeD()==0) {
- cerr << "ERROR: AliPHOSClusterizerv1::ReadDigits There is no Digit Tree" << endl;
- return kFALSE;
- }
-
- //set address of the Digits and Digitizer
- TBranch * digitsBranch = 0;
- TBranch * digitizerBranch = 0;
- TObjArray * lob = (TObjArray*)gAlice->TreeD()->GetListOfBranches() ;
- TIter next(lob) ;
- TBranch * branch = 0 ;
- Bool_t phosfound = kFALSE, digitizerfound = kFALSE ;
-
- TString taskName(GetName()) ;
- taskName.ReplaceAll(Version(), "") ;
-
- while ( (branch = (TBranch*)next()) && (!phosfound || !digitizerfound) ) {
- if ( (strcmp(branch->GetName(), "PHOS")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) ) {
- phosfound = kTRUE ;
- digitsBranch = branch ;
- }
-
- else if ( (strcmp(branch->GetName(), "AliPHOSDigitizer")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) ) {
- digitizerfound = kTRUE ;
- digitizerBranch = branch ;
- }
- }
- if ( !phosfound || !digitizerfound ) {
- cerr << "WARNING: AliPHOSClusterizerv1::ReadDigits -> Digits and/or Digitizer branch with name " << taskName.Data()
- << " not found" << endl ;
- return kFALSE ;
- }
-
- AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
-
- TClonesArray * digits = gime->Digits() ;
- digits->Clear() ;
- digitsBranch->SetAddress(&digits) ;
-
- AliPHOSDigitizer * digitizer = gime->Digitizer() ;
- digitizerBranch->SetAddress(&digitizer) ;
-
- digitsBranch->GetEntry(0) ;
- digitizerBranch->GetEntry(0) ;
-
- fPedestal = digitizer->GetPedestal() ;
- fSlope = digitizer->GetSlope() ;
-
- return kTRUE ;
-}
//____________________________________________________________________________
void AliPHOSClusterizerv1::WriteRecPoints(Int_t event)
// Creates new branches with given title
// fills and writes into TreeR.
- TString branchName(GetName()) ;
- branchName.Remove(branchName.Index(Version())-1) ;
-
+
+ TString branchName(GetName() ) ;
+ branchName.Remove(branchName.Index(Version())-1) ;
+
AliPHOSGetter *gime = AliPHOSGetter::GetInstance() ;
TObjArray * emcRecPoints = gime->EmcRecPoints(branchName) ;
TObjArray * cpvRecPoints = gime->CpvRecPoints(branchName) ;
((AliPHOSRecPoint *)cpvRecPoints->At(index))->SetIndexInList(index) ;
cpvRecPoints->Expand(cpvRecPoints->GetEntriesFast()) ;
-
- gAlice->GetEvent(event) ;
- if(gAlice->TreeR()==0)
- gAlice->MakeTree("R") ;
//Make branches in TreeR for RecPoints and Clusterizer
char * filename = 0;
Int_t index ;
if (( IsInEmc (digit) && Calibrate(digit->GetAmp()) > fEmcClusteringThreshold ) ||
- ( IsInPpsd(digit) && Calibrate(digit->GetAmp()) > fPpsdClusteringThreshold ) ||
( IsInCpv (digit) && Calibrate(digit->GetAmp()) > fCpvClusteringThreshold ) ) {
Int_t iDigitInCluster = 0 ;
} else {
- // start a new PPSD/CPV cluster
+ // start a new CPV cluster
if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
cpvRecPoints->Expand(2*fNumberOfCpvClusters+1);
- if(IsInPpsd(digit))
- cpvRecPoints->AddAt(new AliPHOSPpsdRecPoint(),fNumberOfCpvClusters) ;
- else
- cpvRecPoints->AddAt(new AliPHOSCpvRecPoint(), fNumberOfCpvClusters) ;
+ cpvRecPoints->AddAt(new AliPHOSCpvRecPoint(), fNumberOfCpvClusters) ;
- clu = (AliPHOSPpsdRecPoint *) cpvRecPoints->At(fNumberOfCpvClusters) ;
+ clu = (AliPHOSCpvRecPoint *) cpvRecPoints->At(fNumberOfCpvClusters) ;
fNumberOfCpvClusters++ ;
clu->AddDigit(*digit, Calibrate(digit->GetAmp()) ) ;
clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
} // while digit
delete digitsC ;
-
+
}
//____________________________________________________________________________
void AliPHOSClusterizerv1::MakeUnfolding()
{
// Unfolds clusters using the shape of an ElectroMagnetic shower
- // Performs unfolding of all EMC/CPV but NOT ppsd clusters
+ // Performs unfolding of all EMC/CPV clusters
AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
// Unfold now CPV clusters
if(fNumberOfCpvClusters > 0){
- Int_t nModulesToUnfold = geom->GetNCPVModules() ;
+ Int_t nModulesToUnfold = geom->GetNModules() ;
Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
Int_t index ;
<< " CPV Clustering threshold = " << fCpvClusteringThreshold << endl
<< " CPV Local Maximum cut = " << fCpvLocMaxCut << endl
<< " CPV Logarothmic weight = " << fW0CPV << endl
- << endl
- << " PPSD Clustering threshold = " << fPpsdClusteringThreshold << endl;
+ << endl ;
if(fToUnfold)
cout << " Unfolding on " << endl ;
else
if(strstr(option,"all")) {
cout << "EMC clusters " << endl ;
- cout << " Index "
- << " Ene(MeV) "
- << " Multi "
- << " Module "
- << " X "
- << " Y "
- << " Z "
- << " Lambda 1 "
- << " Lambda 2 "
- << " MaxEnergy "
- << " # of prim "
- << " Primaries list " << endl;
+ cout << " Index Ene(MeV) Multi Module X Y Z Lambda 1 Lambda 2 # of prim Primaries list " << endl;
Int_t index ;
for (index = 0 ; index < emcRecPoints->GetEntries() ; index++) {
AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )emcRecPoints->At(index) ;
-
- cout << setw(6) << rp->GetIndexInList() << " ";
- cout << setw(6) << rp->GetEnergy() << " ";
- cout << setw(6) << rp->GetMultiplicity()<< " ";
- cout << setw(6) << rp->GetPHOSMod() << " ";
-
TVector3 locpos;
rp->GetLocalPosition(locpos);
- cout << setw(8) << locpos.X() << " ";
- cout << setw(8) << locpos.Y() << " ";
- cout << setw(8) << locpos.Z() << " ";
-
Float_t lambda[2];
rp->GetElipsAxis(lambda);
- cout << setw(10)<< lambda[0] << " ";
- cout << setw(10)<< lambda[1] << " ";
-
-
Int_t * primaries;
Int_t nprimaries;
primaries = rp->GetPrimaries(nprimaries);
- cout << setw(8) << primaries << " ";
+ cout << setw(4) << rp->GetIndexInList() << " "
+ << setw(7) << setprecision(3) << rp->GetEnergy() << " "
+ << setw(3) << rp->GetMultiplicity() << " "
+ << setw(1) << rp->GetPHOSMod() << " "
+ << setw(6) << setprecision(2) << locpos.X() << " "
+ << setw(6) << setprecision(2) << locpos.Y() << " "
+ << setw(6) << setprecision(2) << locpos.Z() << " "
+ << setw(4) << setprecision(2) << lambda[0] << " "
+ << setw(4) << setprecision(2) << lambda[1] << " "
+ << setw(2) << nprimaries << " " ;
+
for (Int_t iprimary=0; iprimary<nprimaries; iprimary++)
- cout << setw(4) << primaries[iprimary] << " ";
- cout << endl;
+ cout << setw(4) << primaries[iprimary] << " " ;
+ cout << endl ;
}
- cout << endl ;
- //Now plot CPV/PPSD recPoints
+ //Now plot CPV recPoints
cout << "EMC clusters " << endl ;
cout << " Index "
<< " Multi "
<< " Module "
- << " Layer "
<< " X "
<< " Y "
<< " Z "
for (index = 0 ; index < cpvRecPoints->GetEntries() ; index++) {
AliPHOSRecPoint * rp = (AliPHOSRecPoint * )cpvRecPoints->At(index) ;
cout << setw(6) << rp->GetIndexInList() << " ";
- cout << setw(6) << rp->GetPHOSMod() << " ";
-
- if( (strcmp(rp->ClassName() , "AliPHOSPpsdRecPoint" )) == 0){
- AliPHOSPpsdRecPoint * ppsd = (AliPHOSPpsdRecPoint*) rp ;
- if(ppsd->GetUp())
- cout <<" CPV ";
- else
- cout <<" PPSD ";
- }
- else
- cout <<" CPV ";
+ cout << setw(6) << rp->GetPHOSMod() << " CPV ";
TVector3 locpos;
rp->GetLocalPosition(locpos);
- cout << setw(8) << locpos.X() << " ";
- cout << setw(8) << locpos.Y() << " ";
- cout << setw(8) << locpos.Z() << " ";
+ cout << setw(6) << locpos.X() << " ";
+ cout << setw(6) << locpos.Y() << " ";
+ cout << setw(6) << locpos.Z() << " ";
Int_t * primaries;
- Int_t nprimaries;
+ Int_t nprimaries ;
primaries = rp->GetPrimaries(nprimaries);
- cout << setw(8) << primaries << " ";
+ cout << setw(6) << nprimaries << " ";
for (Int_t iprimary=0; iprimary<nprimaries; iprimary++)
cout << setw(4) << primaries[iprimary] << " ";
}
- cout << "-------------------------------------------------"<<endl ;
+ cout << "-----------------------------------------------------------------------"<<endl ;
}
}
public:
- AliPHOSClusterizerv1() ; // ctor
+ AliPHOSClusterizerv1() ;
AliPHOSClusterizerv1(const char * headerFile, const char * name = "Default");
- virtual ~AliPHOSClusterizerv1(){} // dtor
+ virtual ~AliPHOSClusterizerv1() ;
Int_t AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const ;
// Checks if digits are in neighbour cells
virtual Float_t GetEmcClusteringThreshold()const{ return fEmcClusteringThreshold;}
virtual Float_t GetEmcLocalMaxCut()const { return fEmcLocMaxCut;}
virtual Float_t GetEmcLogWeight()const { return fW0;}
+ virtual Float_t GetEmcTimeGate() const { return fEmcTimeGate ; }
virtual Float_t GetCpvClusteringThreshold()const{ return fCpvClusteringThreshold; }
virtual Float_t GetCpvLocalMaxCut()const { return fCpvLocMaxCut;}
virtual Float_t GetCpvLogWeight()const { return fW0CPV;}
- virtual Float_t GetPpsdClusteringThreshold() const { return fPpsdClusteringThreshold; }
- virtual char * GetRecPointsBranch() const { return (char*) fRecPointsBranchTitle.Data() ;}
- virtual const Int_t GetRecPointsInRun() const {return fRecPointsInRun ;}
- virtual char * GetDigitsBranch() const { return (char*) fDigitsBranchTitle.Data() ;}
+ virtual char * GetRecPointsBranch() const { return (char*) fRecPointsBranchTitle.Data() ;}
+ virtual const Int_t GetRecPointsInRun() const {return fRecPointsInRun ;}
+ virtual char * GetDigitsBranch() const { return (char*) fDigitsBranchTitle.Data() ;}
void Exec(Option_t *option); // Does the job
virtual void SetEmcClusteringThreshold(Float_t cluth) { fEmcClusteringThreshold = cluth ; }
virtual void SetEmcLocalMaxCut(Float_t cut) { fEmcLocMaxCut = cut ; }
virtual void SetEmcLogWeight(Float_t w) { fW0 = w ; }
+ virtual void SetEmcTimeGate(Float_t gate) {fEmcTimeGate = gate ;}
virtual void SetCpvClusteringThreshold(Float_t cluth) { fCpvClusteringThreshold = cluth ; }
virtual void SetCpvLocalMaxCut(Float_t cut) { fCpvLocMaxCut = cut ; }
virtual void SetCpvLogWeight(Float_t w) { fW0CPV = w ; }
- virtual void SetPpsdClusteringThreshold(Float_t cluth) { fPpsdClusteringThreshold = cluth ; }
virtual void SetDigitsBranch(const char * title) { fDigitsBranchTitle = title ;}
virtual void SetRecPointsBranch(const char *title){fRecPointsBranchTitle = title; }
virtual void SetUnfolding(Bool_t toUnfold = kTRUE ) {fToUnfold = toUnfold ;}
void Init() ;
virtual Bool_t IsInEmc (AliPHOSDigit * digit)const ; // Tells if id digit is in EMC
- virtual Bool_t IsInPpsd(AliPHOSDigit * digit)const ; // Tells if id digit is in PPSD
virtual Bool_t IsInCpv (AliPHOSDigit * digit)const ; // Tells if id digit is in CPV
virtual void MakeClusters( ) ;
virtual void MakeUnfolding() ;
- Bool_t ReadDigits(Int_t event) ;
void UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,Int_t Nmax,
int * maxAt,Float_t * maxAtEnergy ) ; //Unfolds cluster using TMinuit package
void WriteRecPoints(Int_t event) ;
Bool_t fToUnfold ; // To perform unfolding
-
Int_t fNumberOfEmcClusters ; // number of EMC clusters found
- Int_t fNumberOfCpvClusters ; // number of CPV+PPSD clusters found
+ Int_t fNumberOfCpvClusters ; // number of CPV clusters found
Float_t fPedestal ; // Calibration parameters
Float_t fSlope ; // read from Digitizer
Float_t fEmcClusteringThreshold ; // minimum energy to include a EMC digit in a cluster
- Float_t fPpsdClusteringThreshold ; // minimum energy to include a PPSD digit in a cluster
Float_t fCpvClusteringThreshold ; // minimum energy to include a CPV digit in a cluster
Float_t fEmcLocMaxCut ; // minimum energy difference to distinguish local maxima in a cluster
Float_t fW0 ; // logarithmic weight for the cluster center of gravity calculation
Float_t fCpvLocMaxCut ; // minimum energy difference to distinguish local maxima in a CPV cluster
Float_t fW0CPV ; // logarithmic weight for the CPV cluster center of gravity calculation
Int_t fRecPointsInRun ; //! Total number of recpoints in one run
-
+ Float_t fEmcTimeGate ; // Maximum time difference between the digits in ont EMC cluster
+
ClassDef(AliPHOSClusterizerv1,1) // Clusterizer implementation version 1
};
#include "TMath.h"
#include "TCanvas.h"
#include "TClonesArray.h"
-#include "AliPHOSGetter.h"
// --- Standard library ---
// --- AliRoot header files ---
#include "AliPHOSCpvRecPoint.h"
-#include "AliPHOSPpsdRecPoint.h"
-
+#include "AliPHOSGetter.h"
ClassImp(AliPHOSCpvRecPoint)
//____________________________________________________________________________
Int_t rv ;
- if( (strcmp(obj->ClassName() , "AliPHOSPpsdRecPoint" )) == 0) // PPSD Rec Point
- {
- AliPHOSPpsdRecPoint * clu = (AliPHOSPpsdRecPoint *)obj ;
- if(this->GetPHOSMod() < clu->GetPHOSMod() )
- rv = -1 ;
- else
- rv = 1 ;
- return rv ;
- }
- else
- {
- AliPHOSCpvRecPoint * clu = (AliPHOSCpvRecPoint *) obj ;
-
- Int_t phosmod1 = GetPHOSMod() ;
- Int_t phosmod2 = clu->GetPHOSMod() ;
-
- TVector3 locpos1;
- GetLocalPosition(locpos1) ;
- TVector3 locpos2;
- clu->GetLocalPosition(locpos2) ;
-
- if(phosmod1 == phosmod2 ) {
- Int_t rowdif = (Int_t)TMath::Ceil(locpos1.X()/delta)-(Int_t)TMath::Ceil(locpos2.X()/delta) ;
- if (rowdif> 0)
- rv = 1 ;
- else if(rowdif < 0)
- rv = -1 ;
- else if(locpos1.Z()>locpos2.Z())
- rv = -1 ;
- else
- rv = 1 ;
- }
-
- else {
- if(phosmod1 < phosmod2 )
- rv = -1 ;
- else
- rv = 1 ;
- }
-
- return rv ;
- }
+ AliPHOSCpvRecPoint * clu = (AliPHOSCpvRecPoint *) obj ;
+
+ Int_t phosmod1 = GetPHOSMod() ;
+ Int_t phosmod2 = clu->GetPHOSMod() ;
+
+ TVector3 locpos1;
+ GetLocalPosition(locpos1) ;
+ TVector3 locpos2;
+ clu->GetLocalPosition(locpos2) ;
+
+ if(phosmod1 == phosmod2 ) {
+ Int_t rowdif = (Int_t)TMath::Ceil(locpos1.X()/delta)-(Int_t)TMath::Ceil(locpos2.X()/delta) ;
+ if (rowdif> 0)
+ rv = 1 ;
+ else if(rowdif < 0)
+ rv = -1 ;
+ else if(locpos1.Z()>locpos2.Z())
+ rv = -1 ;
+ else
+ rv = 1 ;
+ }
+
+ else {
+ if(phosmod1 < phosmod2 )
+ rv = -1 ;
+ else
+ rv = 1 ;
+ }
+
+ return rv ;
+
}
//______________________________________________________________________________
virtual void ExecuteEvent(Int_t event, Int_t px, Int_t py) const ;
void GetClusterLengths(Int_t &lengX, Int_t &lengZ) const {lengX = fLengX ;lengZ = fLengZ ;}
- Bool_t IsEmc(void) const {return kFALSE ; } // tells that this is not a EMC
- Bool_t IsCPV(void) const {return (fPHOSMod <= ((AliPHOSGeometry*) fGeom)->GetNCPVModules()) ; }
- // true if the recpoint is in CPV
+ Bool_t IsEmc(void) const {return kFALSE ; } // tells that this is not a EMC
+ Bool_t IsCPV(void) const {return kTRUE ; } // true if the recpoint is in CPV
Bool_t IsSortable() const { return kTRUE ; } // tells that this is a sortable object
- void Print(Option_t * opt = "void") ;
+ void Print(Option_t * opt = "void") ;
AliPHOSCpvRecPoint & operator = (const AliPHOSCpvRecPoint & rvalue) {
// assignement operator requested by coding convention but not needed
// 3 identifiers for the primary particle(s) at the origine of the digit
// The digits are made in FinishEvent() by summing all the hits in a single PHOS crystal or PPSD gas cell
//
-//*-- Author: Laurent Aphecetche & Yves Schutz (SUBATECH)
+//*-- Author: Laurent Aphecetche & Yves Schutz (SUBATECH) & Dmitri Peressounko (RRC KI & SUBATECH)
// --- ROOT system ---
}
//____________________________________________________________________________
-AliPHOSDigit::AliPHOSDigit(Int_t primary, Int_t id, Int_t digEnergy, Int_t index)
+AliPHOSDigit::AliPHOSDigit(Int_t primary, Int_t id, Int_t digEnergy, Float_t time, Int_t index)
{
// ctor with all data
fNMaxPrimary = 5 ;
fAmp = digEnergy ;
+ fTime = time ;
fId = id ;
fIndexInList = index ;
if( primary != -1){
for ( i = 0; i < fNMaxPrimary ; i++)
fPrimary[i] = digit.fPrimary[i] ;
fAmp = digit.fAmp ;
+ fTime = digit.fTime ;
fId = digit.fId;
fIndexInList = digit.fIndexInList ;
fNprimary = digit.fNprimary ;
// if amplitude is larger than
fAmp += digit.fAmp ;
+ if(fTime > digit.fTime)
+ fTime = digit.fTime ;
Int_t max1 = fNprimary ;
{
// Prints the data of the digit
- out << "ID " << digit.fId << " Energy = " << digit.fAmp << endl ;
+ out << "ID " << digit.fId << " Energy = " << digit.fAmp << " Time = " << digit.fTime << endl ;
Int_t i ;
for(i=0;i<digit.fNprimary;i++)
out << "Primary " << i+1 << " = " << digit.fPrimary[i] << endl ;
public:
AliPHOSDigit() ;
- AliPHOSDigit(Int_t primary, Int_t id, Int_t DigEnergy, Int_t index = -1) ;
+ AliPHOSDigit(Int_t primary, Int_t id, Int_t DigEnergy, Float_t Time, Int_t index = -1) ;
AliPHOSDigit(const AliPHOSDigit & digit) ;
virtual ~AliPHOSDigit() ;
// returns the number of primaries
return fNprimary ; }
Int_t GetPrimary(Int_t index) const ;
+ Float_t GetTime(void) const {return fTime ;}
Bool_t IsSortable() const {
// says that AliPHOSDigits are sortable (needed for Sort method
return kTRUE ; }
void SetAmp(Int_t Amp) {
// sets the amplitude data member
fAmp=Amp ; }
+ void SetTime(Float_t Time) {fTime = Time ;}
void ShiftPrimary(Int_t shift); // shift to semarate different TreeK in merging
private:
Int_t fNprimary ; // Number of primaries
Int_t fNMaxPrimary ; //! Max Number of primaries
- Int_t fPrimary[5] ; // Array of primaries
+ Int_t fPrimary[5] ; // Array of primaries
+ Float_t fTime ;
ClassDef(AliPHOSDigit,1) // Digit in PHOS
fEMCDigitThreshold = 0.01 ;
fCPVNoise = 0.01;
fCPVDigitThreshold = 0.09 ;
- fPPSDNoise = 0.0000001;
- fPPSDDigitThreshold = 0.0000002 ;
+ fTimeResolution = 1.0e-9 ;
fDigitsInRun = 0 ;
fPedestal = 0.; // Calibration parameters
fSlope = 10000000. ;
fEMCDigitThreshold = 0.01 ;
fCPVNoise = 0.01;
fCPVDigitThreshold = 0.09 ;
- fPPSDNoise = 0.0000001;
- fPPSDDigitThreshold = 0.0000002 ;
fDigitsInRun = 0 ;
fPedestal = 0.; // Calibration parameters
fSlope = 10000000. ;
fEMCDigitThreshold = 0.01 ;
fCPVNoise = 0.01;
fCPVDigitThreshold = 0.09 ;
- fPPSDNoise = 0.0000001;
- fPPSDDigitThreshold = 0.0000002 ;
fDigitsInRun = 0 ;
fPedestal = 0.; // Calibration parameters
fSlope = 10000000. ;
Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
Int_t nCPV ;
- Int_t nPPSD ;
Int_t absID ;
TString name = geom->GetName() ;
- if ( name == "IHEP" || name == "MIXT" )
- nCPV =nEMC + geom->GetNumberOfCPVPadsZ()*geom->GetNumberOfCPVPadsPhi()*
- geom->GetNCPVModules()*geom->GetNumberOfCPVLayers() ;
- else
- nCPV = nEMC;
-
- if ( name == "GPS2" || name == "MIXT" )
- nPPSD =nCPV+2*geom->GetNPPSDModules()*geom->GetNumberOfModulesPhi()*geom->GetNumberOfModulesZ()*
- geom->GetNumberOfPadsPhi()*geom->GetNumberOfPadsZ() ;
- else
- nPPSD = nCPV;
+ nCPV = nEMC + geom->GetNumberOfCPVPadsZ()*geom->GetNumberOfCPVPadsPhi()*
+ geom->GetNModules() ;
-
- digits->Expand(nPPSD) ;
+ digits->Expand(nCPV) ;
// sdigitize random gaussian noise and add it to all cells (EMCA+CPV+PPSD)
cerr << "ERROR: AliPHOSDigitizer::Digitize -> SDigitizer with name " << GetName() << " not found " << endl ;
abort() ;
}
- for(absID = 1; absID <= nEMC; absID++){
- Float_t noise = gRandom->Gaus(0., fPinNoise) ;
- new((*digits)[absID-1]) AliPHOSDigit( -1,absID,sDigitizer->Digitize(noise) ) ;
- }
-
- for(absID = nEMC+1; absID <= nCPV; absID++){
- Float_t noise = gRandom->Gaus(0., fCPVNoise) ;
- new((*digits)[absID-1]) AliPHOSDigit( -1,absID,sDigitizer->Digitize(noise) ) ;
- }
-
- for(absID = nCPV+1; absID <= nPPSD; absID++){
- Float_t noise = gRandom->Gaus(0., fPPSDNoise) ;
- new((*digits)[absID-1]) AliPHOSDigit( -1,absID,sDigitizer->Digitize(noise) ) ;
- }
-
+
// loop through the sdigits posted to the White Board and add them to the noise
TCollection * folderslist = ((TFolder*)gROOT->FindObjectAny("YSAlice/WhiteBoard/SDigits/PHOS"))->GetListOfFolders() ;
TIter next(folderslist) ;
TFolder * folder = 0 ;
TClonesArray * sdigits = 0 ;
Int_t input = 0 ;
- while ( (folder = (TFolder*)next()) ) {
+ TObjArray * sdigArray = new TObjArray(2) ;
+ while ( (folder = (TFolder*)next()) )
if ( (sdigits = (TClonesArray*)folder->FindObject(GetName()) ) ) {
cout << "INFO: AliPHOSDigitizer::Digitize -> Adding SDigits "
<< GetName() << " from " << folder->GetName() << endl ;
- Int_t primaryoffset ;
- if(fARD)
- primaryoffset = fARD->GetMask(input) ;
- else
- primaryoffset = 10000000*input ;
-
- Int_t index ;
- AliPHOSDigit * curSDigit ;
- AliPHOSDigit * digit ;
- for ( index = 0 ; index < sdigits->GetEntriesFast(); index++) {
- curSDigit = (AliPHOSDigit*)sdigits->At(index) ;
- curSDigit->ShiftPrimary(primaryoffset) ;
- digit = (AliPHOSDigit*)digits->At(curSDigit->GetId() - 1 ) ;
- *digit = *digit + *curSDigit ;
+ sdigArray->AddAt(sdigits, input) ;
+ input++ ;
+ }
+
+ //Find the first crystall with sygnal
+ Int_t nextSig = 200000 ;
+ Int_t i;
+ for(i=0; i<input; i++){
+ sdigits = (TClonesArray *)sdigArray->At(i) ;
+ Int_t curNext = ((AliPHOSDigit *)sdigits->At(0))->GetId() ;
+ if(curNext < nextSig) nextSig = curNext ;
+ }
+
+ TArrayF * energies = new TArrayF(input);
+ TArrayF * times = new TArrayF(input) ;
+ TArrayI index(input) ;
+ index.Reset() ; //Set all indexes to zero
+
+ AliPHOSDigit * digit ;
+ AliPHOSDigit * curSDigit ;
+
+ //Put Noise contribution
+ for(absID = 1; absID <= nEMC; absID++){
+ Float_t noise = gRandom->Gaus(0., fPinNoise) ;
+ new((*digits)[absID-1]) AliPHOSDigit( -1,absID,sDigitizer->Digitize(noise), TimeOfNoise() ) ;
+ //look if we have to add signal?
+ if(absID==nextSig){
+ //Add SDigits from all inputs
+ digit = (AliPHOSDigit *) digits->At(absID-1) ;
+ Int_t contrib = 0 ;
+ energies[contrib] = digit->GetAmp() ;
+ times[contrib] = digit->GetTime() ;
+ //loop over inputs
+ for(i=0; i<input; i++){
+ curSDigit = (AliPHOSDigit*)((TClonesArray *)sdigArray->At(i))->At(index[i]) ;
+ //May be several digits will contribute from the same input
+ while(curSDigit && curSDigit->GetId() == absID){
+ //Shift primary to separate primaries belonging different inputs
+ Int_t primaryoffset ;
+ if(fARD)
+ primaryoffset = fARD->GetMask(i) ;
+ else
+ primaryoffset = 10000000*i ;
+ curSDigit->ShiftPrimary(primaryoffset) ;
+
+ *digit = *digit + *curSDigit ; //add energies
+ energies[contrib] = curSDigit->GetAmp() ;
+ times[contrib] = curSDigit->GetTime() ;
+ contrib++ ;
+ index[i]++ ;
+ curSDigit = (AliPHOSDigit*)((TClonesArray *)sdigArray->At(i))->At(index[i]) ;
+ }
+ }
+
+ //calculate and set time
+ Float_t time = FrontEdgeTime(energies, times) ;
+ digit->SetTime(time) ;
+ energies->Reset() ;
+ times->Reset() ;
+ //Find next signal module
+ for(i=0; i<input; i++){
+ sdigits = ((TClonesArray *)sdigArray->At(i)) ;
+ Int_t curNext = ((AliPHOSDigit *) sdigits->At(index[i]))->GetId() ;
+ if(curNext < nextSig) nextSig = curNext ;
}
}
}
+
+
+ //Now CPV digits (different noise and no timing)
+ for(absID = nEMC+1; absID <= nCPV; absID++){
+ Float_t noise = gRandom->Gaus(0., fCPVNoise) ;
+ new((*digits)[absID-1]) AliPHOSDigit( -1,absID,sDigitizer->Digitize(noise), TimeOfNoise() ) ;
+ //look if we have to add signal?
+ if(absID==nextSig){
+ digit = (AliPHOSDigit *) digits->At(absID-1) ;
+ //Add SDigits from all inputs
+ for(i=0; i<input; i++){
+ curSDigit = (AliPHOSDigit*)((TClonesArray *)sdigArray->At(i))->At(index[i]) ;
+ //May be several digits will contribute from the same input
+ while(curSDigit && curSDigit->GetId() == absID){
+ //Shift primary to separate primaries belonging different inputs
+ Int_t primaryoffset ;
+ if(fARD)
+ primaryoffset = fARD->GetMask(i) ;
+ else
+ primaryoffset = 10000000*i ;
+ curSDigit->ShiftPrimary(primaryoffset) ;
+
+ //add energies
+ *digit = *digit + *curSDigit ;
+ index[i]++ ;
+ curSDigit = (AliPHOSDigit*)((TClonesArray *)sdigArray->At(i))->At(index[i]) ;
+ }
+ }
+
+ //Find next signal module
+ for(i=0; i<input; i++){
+ sdigits = (TClonesArray *)sdigArray->At(i) ;
+ Int_t curNext = ((AliPHOSDigit *) sdigits->At(index[i]))->GetId() ;
+ if(curNext < nextSig) nextSig = curNext ;
+ }
+
+ }
+ }
+ delete energies ;
+ delete times ;
+ delete sdigArray ; //We should not delete its contents
+
+
+
//remove digits below thresholds
for(absID = 0; absID < nEMC ; absID++)
if(sDigitizer->Calibrate(((AliPHOSDigit*)digits->At(absID))->GetAmp()) < fEMCDigitThreshold)
for(absID = nEMC; absID < nCPV ; absID++)
if(sDigitizer->Calibrate(((AliPHOSDigit*)digits->At(absID))->GetAmp()) < fCPVDigitThreshold)
digits->RemoveAt(absID) ;
-
- for(absID = nCPV; absID < nPPSD ; absID++)
- if(sDigitizer->Calibrate(((AliPHOSDigit *)digits->At(absID))->GetAmp()) < fPPSDDigitThreshold)
- digits->RemoveAt(absID) ;
-
+
digits->Compress() ;
Int_t ndigits = digits->GetEntriesFast() ;
digits->Expand(ndigits) ;
//Set indexes in list of digits
- Int_t i ;
for (i = 0 ; i < ndigits ; i++) {
AliPHOSDigit * digit = (AliPHOSDigit *) digits->At(i) ;
digit->SetIndexInList(i) ;
}
+//____________________________________________________________________________
+Float_t AliPHOSDigitizer::FrontEdgeTime(TArrayF * energies, TArrayF * times)
+{ //
+ Float_t curtime = times->At(0) ;
+ Float_t time = curtime ;
+ Int_t i = 1 ;
+ while(curtime){
+ if(time > curtime) time = curtime ;
+ curtime = times->At(i++) ;
+ }
+ return time ;
+
+}
//____________________________________________________________________________
Bool_t AliPHOSDigitizer::Init()
{
cout << " Threshold in EMC (fEMCDigitThreshold) = " << fEMCDigitThreshold << endl ; ;
cout << " Noise in CPV (fCPVNoise) = " << fCPVNoise << endl ;
cout << " Threshold in CPV (fCPVDigitThreshold) = " << fCPVDigitThreshold << endl ;
- cout << " Noise in PPSD (fPPSDNoise) = " << fPPSDNoise << endl ;
- cout << " Threshold in PPSD (fPPSDDigitThreshold) = " << fPPSDDigitThreshold << endl ;
cout << "---------------------------------------------------" << endl ;
}
else
cout << "AliPHOSDigitizer not initialized " << endl ;
}
+
//__________________________________________________________________
void AliPHOSDigitizer::PrintDigits(Option_t * option){
// Print a table of digits
cout << "AliPHOSDigitiser: event " << gAlice->GetEvNumber() << endl ;
cout << " Number of entries in Digits list " << digits->GetEntriesFast() << endl ;
cout << endl ;
-
- fDigitsInRun += digits->GetEntriesFast() ;
-
- if(strstr(option,"all")){
+ if(strstr(option,"all")||strstr(option,"EMC")){
//loop over digits
AliPHOSDigit * digit;
- cout << "Digit Id " << " Amplitude " << " Index " << " Nprim " << " Primaries list " << endl;
+ cout << "Digit Id Amplitude Index " << " Nprim " << " Primaries list " << endl;
+ Int_t maxEmc = gime->PHOSGeometry()->GetNModules()*gime->PHOSGeometry()->GetNCristalsInModule() ;
Int_t index ;
- for (index = 0 ; index < digits->GetEntries() ; index++) {
+ for (index = 0 ; (index < digits->GetEntries()) &&
+ (((AliPHOSDigit * ) digits->At(index))->GetId() <= maxEmc) ; index++) {
digit = (AliPHOSDigit * ) digits->At(index) ;
- cout << setw(8) << digit->GetId() << " " << setw(3) << digit->GetAmp() << " "
- << setw(6) << digit->GetIndexInList() << " "
- << setw(5) << digit->GetNprimary() <<" ";
+ if(digit->GetNprimary() == 0) continue;
+ cout << setw(6) << digit->GetId() << " " << setw(10) << digit->GetAmp() << " "
+ << setw(6) << digit->GetIndexInList() << " "
+ << setw(5) << digit->GetNprimary() <<" ";
Int_t iprimary;
for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++)
- cout << setw(5) << digit->GetPrimary(iprimary+1) << " ";
+ cout << setw(5) << digit->GetPrimary(iprimary+1) << " ";
cout << endl;
- }
+ }
+ cout << endl;
+ }
+
+ if(strstr(option,"all")||strstr(option,"CPV")){
+ //loop over CPV digits
+ AliPHOSDigit * digit;
+ cout << "Digit Id " << " Amplitude " << " Index " << " Nprim " << " Primaries list " << endl;
+ Int_t maxEmc = gime->PHOSGeometry()->GetNModules()*gime->PHOSGeometry()->GetNCristalsInModule() ;
+ Int_t index ;
+ for (index = 0 ; index < digits->GetEntries(); index++) {
+ digit = (AliPHOSDigit * ) digits->At(index) ;
+ if(digit->GetId() > maxEmc){
+ cout << setw(6) << digit->GetId() << " " << setw(10) << digit->GetAmp() << " "
+ << setw(6) << digit->GetIndexInList() << " "
+ << setw(5) << digit->GetNprimary() <<" ";
+
+ Int_t iprimary;
+ for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++)
+ cout << setw(5) << digit->GetPrimary(iprimary+1) << " ";
+ cout << endl;
+ }
+ }
}
+
}
//__________________________________________________________________
AliPHOSGetter::GetInstance()->SDigits()->SetName(title) ;
+}
+//__________________________________________________________________
+Float_t AliPHOSDigitizer::TimeOfNoise(void)
+{ // Calculates the time signal generated by noise
+ //to be rewritten, now returns just big number
+ return 1. ;
+
}
//____________________________________________________________________________
void AliPHOSDigitizer::Reset()
const Float_t GetEMCThreshold() const { return fEMCDigitThreshold;}
const Float_t GetPedestal() const { return fPedestal; }
const Float_t GetPinNoise() const { return fPinNoise;}
- const Float_t GetPPSDNoise() const { return fPPSDNoise ;}
- const Float_t GetPPSDThreshold()const { return fPPSDDigitThreshold ;}
const Float_t GetSlope() const { return fSlope; }
- // const TArrayI * GetCurrentEvents()const { return fIevent ;}
+ const Float_t GetTimeResolution() const { return fTimeResolution ; }
void MixWith(const char* HeaderFile) ; // Add another one file to mix
void Print(Option_t* option)const ;
void SetCPVThreshold(Float_t CPVThreshold) {fCPVDigitThreshold= CPVThreshold;}
void SetEMCThreshold(Float_t EMCThreshold) {fEMCDigitThreshold = EMCThreshold;}
void SetPinNoise(Float_t PinNoise ) {fPinNoise = PinNoise;}
- void SetPPSDNoise(Float_t PPSDNoise) {fPPSDNoise = PPSDNoise;}
- void SetPPSDThreshold(Float_t PPSDThreshold){fPPSDDigitThreshold = PPSDThreshold;}
void SetSDigitsBranch(const char* file) ;
Bool_t Init() ;
void PrintDigits(Option_t * option) ;
void WriteDigits(Int_t evt) ; // Writes Digits for particular event
-
+ Float_t TimeOfNoise(void) ; // Calculate time signal generated by noise
+ Float_t FrontEdgeTime(TArrayF *energies, TArrayF * times) ;
+ Int_t DigitizeEnergy(Int_t amp, Int_t absId) ;
+ //Calculate the time of crossing of the threshold by front edge
private:
- Float_t fPedestal ; // Calibration parameters
- Float_t fSlope ; // read from SDigitizer
+ Float_t fPedestal ; // Calibration parameters read from SDigitizer
+ Float_t fSlope ; // (nothing common with real digitization parameters)
AliRunDigitizer * fARD ; //! Pointer to the Digitization Manager class
+ Int_t fEmcCrystals ; // Number of EMC crystalls in the given geometry
Float_t fPinNoise ; // Electronics noise in EMC
Float_t fEMCDigitThreshold ; // Threshold for storing digits in EMC
Float_t fCPVNoise ; // Noise in CPV
Float_t fCPVDigitThreshold ; // Threshold for storing digits in CPV
- Float_t fPPSDNoise ; // Noise in PPSD
- Float_t fPPSDDigitThreshold ; // Threshold for storing digits in PPSD
Int_t fDigitsInRun ; //! Total number of digits in one run
+ Float_t fTimeResolution ; // Time resolution of FEE electronics
+
+ Float_t fTimeThreshold ; // Threshold to start timing for given crystall
+ Float_t fTimeSignalLength ; // Length of the timing signal
+
+ Int_t fNADCbitsEmc ;
+ Float_t fADCpedestalEmc ;
+ Int_t fNADCbitsCPV ;
+ Float_t fADCperestalCPV ;
+
ClassDef(AliPHOSDigitizer,1) // description
fAmp = 0. ;
fCoreEnergy = 0 ;
fEnergyList = 0 ;
+ fTime = -1. ;
fLocPos.SetX(1000000.) ; //Local position should be evaluated
}
return rv ;
}
-
//______________________________________________________________________________
void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py) const
{
-// Commented by Dmitri Peressounko: there is no possibility to ensure,
-// that AliPHOSGetter keeps the correct information.
-
-// // Execute action corresponding to one event
-// // This member function is called when a AliPHOSRecPoint is clicked with the locator
-// //
-// // If Left button is clicked on AliPHOSRecPoint, the digits are switched on
-// // and switched off when the mouse button is released.
-// //
-
-// // static Int_t pxold, pyold;
-
-// AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
-// static TGraph * digitgraph = 0 ;
+ // Execute action corresponding to one event
+ // This member function is called when a AliPHOSRecPoint is clicked with the locator
+ //
+ // If Left button is clicked on AliPHOSRecPoint, the digits are switched on
+ // and switched off when the mouse button is released.
+
+
+ AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
+ if(!gime) return ;
+ AliPHOSGeometry * phosgeom = (AliPHOSGeometry*)gime->PHOSGeometry();
-// if (!gPad->IsEditable()) return;
+ static TGraph * digitgraph = 0 ;
-// TH2F * histo = 0 ;
-// TCanvas * histocanvas ;
+ if (!gPad->IsEditable()) return;
-// switch (event) {
+ TH2F * histo = 0 ;
+ TCanvas * histocanvas ;
+
+ TClonesArray * digits = gime->Digits() ;
+
+ switch (event) {
-// case kButton1Down: {
-// AliPHOSDigit * digit ;
-// AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
-// AliPHOSGeometry * phosgeom = (AliPHOSGeometry*)gime->PHOSGeometry();
-// Int_t iDigit;
-// Int_t relid[4] ;
+ case kButton1Down: {
+ AliPHOSDigit * digit ;
+ Int_t iDigit;
+ Int_t relid[4] ;
-// const Int_t kMulDigit = AliPHOSEmcRecPoint::GetDigitsMultiplicity() ;
-// Float_t * xi = new Float_t[kMulDigit] ;
-// Float_t * zi = new Float_t[kMulDigit] ;
+ const Int_t kMulDigit = AliPHOSEmcRecPoint::GetDigitsMultiplicity() ;
+ Float_t * xi = new Float_t[kMulDigit] ;
+ Float_t * zi = new Float_t[kMulDigit] ;
-// // create the histogram for the single cluster
-// // 1. gets histogram boundaries
-// Float_t ximax = -999. ;
-// Float_t zimax = -999. ;
-// Float_t ximin = 999. ;
-// Float_t zimin = 999. ;
+ // create the histogram for the single cluster
+ // 1. gets histogram boundaries
+ Float_t ximax = -999. ;
+ Float_t zimax = -999. ;
+ Float_t ximin = 999. ;
+ Float_t zimin = 999. ;
-// for(iDigit=0; iDigit<kMulDigit; iDigit++) {
-// digit = (AliPHOSDigit *) ( gime->Digit(fDigitsList[iDigit]) ) ;
-// phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
-// phosgeom->RelPosInModule(relid, xi[iDigit], zi[iDigit]);
-// if ( xi[iDigit] > ximax )
-// ximax = xi[iDigit] ;
-// if ( xi[iDigit] < ximin )
-// ximin = xi[iDigit] ;
-// if ( zi[iDigit] > zimax )
-// zimax = zi[iDigit] ;
-// if ( zi[iDigit] < zimin )
-// zimin = zi[iDigit] ;
-// }
-// ximax += phosgeom->GetCrystalSize(0) / 2. ;
-// zimax += phosgeom->GetCrystalSize(2) / 2. ;
-// ximin -= phosgeom->GetCrystalSize(0) / 2. ;
-// zimin -= phosgeom->GetCrystalSize(2) / 2. ;
-// Int_t xdim = (int)( (ximax - ximin ) / phosgeom->GetCrystalSize(0) + 0.5 ) ;
-// Int_t zdim = (int)( (zimax - zimin ) / phosgeom->GetCrystalSize(2) + 0.5 ) ;
+ for(iDigit=0; iDigit<kMulDigit; iDigit++) {
+ digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit]) ;
+ phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ phosgeom->RelPosInModule(relid, xi[iDigit], zi[iDigit]);
+ if ( xi[iDigit] > ximax )
+ ximax = xi[iDigit] ;
+ if ( xi[iDigit] < ximin )
+ ximin = xi[iDigit] ;
+ if ( zi[iDigit] > zimax )
+ zimax = zi[iDigit] ;
+ if ( zi[iDigit] < zimin )
+ zimin = zi[iDigit] ;
+ }
+ ximax += phosgeom->GetCrystalSize(0) / 2. ;
+ zimax += phosgeom->GetCrystalSize(2) / 2. ;
+ ximin -= phosgeom->GetCrystalSize(0) / 2. ;
+ zimin -= phosgeom->GetCrystalSize(2) / 2. ;
+ Int_t xdim = (int)( (ximax - ximin ) / phosgeom->GetCrystalSize(0) + 0.5 ) ;
+ Int_t zdim = (int)( (zimax - zimin ) / phosgeom->GetCrystalSize(2) + 0.5 ) ;
-// // 2. gets the histogram title
+ // 2. gets the histogram title
-// Text_t title[100] ;
-// sprintf(title,"Energy=%1.2f GeV ; Digits ; %d ", GetEnergy(), GetDigitsMultiplicity()) ;
+ Text_t title[100] ;
+ sprintf(title,"Energy=%1.2f GeV ; Digits ; %d ", GetEnergy(), GetDigitsMultiplicity()) ;
-// if (!histo) {
-// delete histo ;
-// histo = 0 ;
-// }
-// histo = new TH2F("cluster3D", title, xdim, ximin, ximax, zdim, zimin, zimax) ;
+ if (!histo) {
+ delete histo ;
+ histo = 0 ;
+ }
+ histo = new TH2F("cluster3D", title, xdim, ximin, ximax, zdim, zimin, zimax) ;
-// Float_t x, z ;
-// for(iDigit=0; iDigit<kMulDigit; iDigit++) {
-// digit = (AliPHOSDigit *) ( gime->Digit(fDigitsList[iDigit]) ) ;
-// phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
-// phosgeom->RelPosInModule(relid, x, z);
-// histo->Fill(x, z, fEnergyList[iDigit] ) ;
-// }
+ Float_t x, z ;
+ for(iDigit=0; iDigit<kMulDigit; iDigit++) {
+ digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit]) ;
+ phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ phosgeom->RelPosInModule(relid, x, z);
+ histo->Fill(x, z, fEnergyList[iDigit] ) ;
+ }
-// if (!digitgraph) {
-// digitgraph = new TGraph(kMulDigit,xi,zi);
-// digitgraph-> SetMarkerStyle(5) ;
-// digitgraph-> SetMarkerSize(1.) ;
-// digitgraph-> SetMarkerColor(1) ;
-// digitgraph-> Paint("P") ;
-// }
+ if (!digitgraph) {
+ digitgraph = new TGraph(kMulDigit,xi,zi);
+ digitgraph-> SetMarkerStyle(5) ;
+ digitgraph-> SetMarkerSize(1.) ;
+ digitgraph-> SetMarkerColor(1) ;
+ digitgraph-> Paint("P") ;
+ }
-// Print() ;
-// histocanvas = new TCanvas("cluster", "a single cluster", 600, 500) ;
-// histocanvas->Draw() ;
-// histo->Draw("lego1") ;
+ // Print() ;
+ histocanvas = new TCanvas("cluster", "a single cluster", 600, 500) ;
+ histocanvas->Draw() ;
+ histo->Draw("lego1") ;
-// delete[] xi ;
-// delete[] zi ;
+ delete[] xi ;
+ delete[] zi ;
-// break;
-// }
+ break;
+ }
-// case kButton1Up:
-// if (digitgraph) {
-// delete digitgraph ;
-// digitgraph = 0 ;
-// }
-// break;
+ case kButton1Up:
+ if (digitgraph) {
+ delete digitgraph ;
+ digitgraph = 0 ;
+ }
+ break;
-// }
+ }
}
//____________________________________________________________________________
}
return iDigitN ;
}
-
-
-
+//____________________________________________________________________________
+void AliPHOSEmcRecPoint::EvalTime(TClonesArray * digits){
+
+ Float_t maxE = 0;
+ Int_t maxAt = 0;
+ for(Int_t idig=0; idig < fMulDigit; idig++){
+ if(fEnergyList[idig] > maxE){
+ maxE = fEnergyList[idig] ;
+ maxAt = idig;
+ }
+ }
+ fTime = ((AliPHOSDigit*) digits->At(fDigitsList[maxAt]))->GetTime() ;
+
+}
//____________________________________________________________________________
void AliPHOSEmcRecPoint::Print(Option_t * option)
{
Int_t Compare(const TObject * obj) const; // method for sorting
virtual void EvalAll(Float_t logWeight,TClonesArray * digits) ;
- void EvalCoreEnergy(Float_t logWeight,TClonesArray * digits) ;
- virtual void EvalLocalPosition(Float_t logWeight,TClonesArray * digits) ;// computes the position in the PHOS module
- virtual void EvalDispersion(Float_t logWeight,TClonesArray * digits) ; // computes the dispersion of the shower
- virtual void EvalElipsAxis(Float_t logWeight, TClonesArray * digits ); // computes the axis of shower ellipsoide
+
virtual void ExecuteEvent(Int_t event, Int_t px, Int_t py) const;
Float_t GetCoreEnergy()const {return fCoreEnergy ;}
virtual Int_t GetNumberOfLocalMax(Int_t * maxAt, Float_t * maxAtEnergy,
Float_t locMaxCut,TClonesArray * digits ) const ;
// searches for the local maxima
+ Float_t GetTime(void) const{return fTime ; }
Bool_t IsEmc(void) const { return kTRUE ; } // true if the recpoint is in EMC
Bool_t IsSortable() const {return kTRUE ; } // says that emcrecpoints are sortable objects
void Print(Option_t * opt = "void") ;
}
protected:
-
+ void EvalCoreEnergy(Float_t logWeight,TClonesArray * digits) ;
+ virtual void EvalLocalPosition(Float_t logWeight,TClonesArray * digits) ;// computes the position in the PHOS module
+ virtual void EvalDispersion(Float_t logWeight,TClonesArray * digits) ; // computes the dispersion of the shower
+ virtual void EvalElipsAxis(Float_t logWeight, TClonesArray * digits ); // computes the axis of shower ellipsoide
+ void EvalTime( TClonesArray * digits );
virtual Bool_t AreNeighbours(AliPHOSDigit * digit1, AliPHOSDigit * digit2 ) const ;
Float_t fCoreEnergy ; // energy in a shower core
Float_t fLambda[2] ; // shower ellipse axes
Float_t fDispersion ; // shower dispersion
Float_t *fEnergyList ; //[fMulDigit] energy of digits
+ Float_t fTime ; // Time of the digit with maximal energy deposition
ClassDef(AliPHOSEmcRecPoint,1) // EMC RecPoint (cluster)
TString name ;
switch (fType) {
- case kGAMMA:
+ case kNEUTRALEMFAST:
name = "PHOTON" ;
break ;
- case kELECTRON:
+ case kCHARGEDEMFAST:
name = "ELECTRON" ;
break ;
- case kCHARGEDHA:
- name = "CHARGED_HA" ;
+ case kCHARGEDHAFAST:
+ name = "CHARGED_HA_FAST" ;
break ;
- case kNEUTRALHA:
- name = "NEUTRAL_HA" ;
+ case kNEUTRALHASLOW:
+ name = "NEUTRAL_HA_SLOW" ;
break ;
- case kNEUTRALEM:
- name = "NEUTRAL_EM" ;
+ case kNEUTRALEMSLOW:
+ name = "NEUTRAL_EM_SLOW" ;
break ;
- case kGAMMAHA:
- name = "PHOTON_HA" ;
+ case kNEUTRALHAFAST:
+ name = "NEUTRAL_HA_FAST" ;
break ;
}
fIndexInList = val ;
}
- enum EParticleType { kUNDEFINED=-1, kNEUTRALEM, kNEUTRALHA, kGAMMA , kGAMMAHA ,
- kABSURDEM, kABSURDHA , kELECTRON, kCHARGEDHA } ;
+ enum EParticleType { kUNDEFINED=-1,
+ kNEUTRALEMFAST, kNEUTRALHAFAST, kNEUTRALEMSLOW, kNEUTRALHASLOW,
+ kCHARGEDEMFAST, kCHARGEDHAFAST, kCHARGEDEMSLOW, kCHARGEDHASLOW } ;
typedef TClonesArray FastRecParticlesList ;
// Hits class for PHOS
// A hit in PHOS is the sum of all hits in a single crystal
//*--
-//*-- Author: Maxime Volkov (RRC KI) & Yves Schutz (SUBATECH)
+//*-- Author: Maxime Volkov (RRC KI) & Yves Schutz (SUBATECH) & Dmitri Peressounko (RRC KI & SUBATECH)
// --- ROOT system ---
#include "AliPHOSGeometry.h"
#include "AliPHOS.h"
-
-
ClassImp(AliPHOSHit)
-
-//____________________________________________________________________________
-AliPHOSHit::AliPHOSHit(const AliPHOSHit & hit)
+
+ //____________________________________________________________________________
+ AliPHOSHit::AliPHOSHit(const AliPHOSHit & hit)
{
- // copy ctor
-
+ // copy ctor
fX = hit.fX ;
fY = hit.fY ;
fZ = hit.fZ ;
fELOS = hit.fELOS ;
fPrimary = hit.fPrimary ;
fTrack = hit.fTrack ;
-
-
+ fTime = hit.fTime ;
+
}
//____________________________________________________________________________
//
// Create a CPV hit object
//
-
- fId = id ;
+
fX = hits[0] ;
fY = hits[1] ;
fZ = hits[2] ;
- fELOS = hits[3] ;
+ fTime = hits[3] ;
+ fId = id ;
+ fELOS = hits[4] ;
fPrimary = primary ;
}
-
//____________________________________________________________________________
Float_t AliPHOSHit::X() const
{
+ // if(fX < -1000.){
TVector3 pos ;
AliPHOS * phos = static_cast<AliPHOS*> (gAlice->GetDetector("PHOS")) ;
phos->GetGeometry() ->RelPosInAlice(GetId(), pos) ;
return pos.X() ;
+ // fX = pos.X() ;
+ // fY = pos.Y() ;
+ // fZ = pos.Z() ;
+ // }
+ // return fX;
}
-
//____________________________________________________________________________
Float_t AliPHOSHit::Y() const
{
+ // if(fY < -1000.){
TVector3 pos ;
AliPHOS * phos = static_cast<AliPHOS*> (gAlice->GetDetector("PHOS")) ;
phos->GetGeometry() ->RelPosInAlice(GetId(), pos) ;
- return pos.Y() ;
+ return pos.Y();
+ // fX = pos.X() ;
+ // fY = pos.Y() ;
+ // fZ = pos.Z() ;
+ // }
+ // return fY;
}
//____________________________________________________________________________
Float_t AliPHOSHit::Z() const
{
+ // if(fY < -1000.){
TVector3 pos ;
AliPHOS * phos = static_cast<AliPHOS*> (gAlice->GetDetector("PHOS")) ;
phos->GetGeometry() ->RelPosInAlice(GetId(), pos) ;
return pos.Z() ;
+ // fX = pos.X() ;
+ // fY = pos.Y() ;
+ // fZ = pos.Z() ;
+ // }
+ // return fZ;
}
-
//____________________________________________________________________________
Bool_t AliPHOSHit::operator==(AliPHOSHit const &rValue) const
{
// Add the energy of the hit
fELOS += rValue.GetEnergy() ;
+
+ if(rValue.GetTime() < fTime)
+ fTime = rValue.GetTime() ;
return *this;
{
// Print out Id and energy
- out << "AliPHOSHit = " << hit.GetId() << " " << hit.GetEnergy() << endl ;
+ out << "AliPHOSHit = " << hit.GetId() << " " << hit.GetEnergy() << " " << hit.GetTime() << endl ;
return out ;
}
// returns the primary particle id at the origine of this hit
return fPrimary ;
}
+
+ Float_t GetTime(void) const {
+ // returns the time of the first energy deposition
+ return fTime ;
+ }
+
virtual Float_t X() const ;
virtual Float_t Y() const ;
virtual Float_t Z() const ;
-
Bool_t operator == (AliPHOSHit const &rValue) const ;
AliPHOSHit operator + (const AliPHOSHit& rValue) ;
Int_t fId ; // Absolute Id number of PHOS Xtal or PPSD pad
Float_t fELOS ; // Energy deposited
Int_t fPrimary ; // Primary particles at the origine of the hit
+ Float_t fTime ; // Time of the energy deposition
ClassDef(AliPHOSHit,1) // Hit for PHOS
virtual void Init()= 0 ;
virtual void Print(Option_t * option) const = 0 ;
virtual void PlotDispersionCuts()const = 0;
- virtual Bool_t ReadTrackSegments(Int_t event)= 0 ;
virtual void SetIdentificationMethod(char * option) = 0 ;
virtual void SetShowerProfileCut(char * formula) = 0 ;
virtual void SetDispersionCut(Float_t cut) = 0 ;
virtual void SetCpvtoEmcDistanceCut(Float_t cut ) = 0;
+ virtual void SetTimeGate(Float_t gate) = 0 ;
virtual void SetTrackSegmentsBranch(const char* title) = 0 ;
virtual void SetRecParticlesBranch (const char* title) = 0 ;
virtual const char * Version() const = 0 ;
//
//*-- Author: Yves Schutz (SUBATECH) & Gines Martinez (SUBATECH) &
// Dmitri Peressounko (SUBATECH & Kurchatov Institute)
-// Complitely redesined by Dmitri Peressounko, March 2001
+// Completely redesined by Dmitri Peressounko, March 2001
// --- ROOT system ---
#include "TROOT.h"
fFormula = 0 ;
fDispersion = 0. ;
fCpvEmcDistance = 0 ;
+ fTimeGate = 2.e-9 ;
fHeaderFileName = "" ;
fTrackSegmentsTitle= "" ;
fRecPointsTitle = "" ;
fRecParticlesTitle = "" ;
- fIDOptions = "" ;
+ fIDOptions = "dis time" ;
fRecParticlesInRun = 0 ;
}
fFormula = new TFormula("LambdaCuts","(x>1)*(x<2.5)*(y>0)*(y<x)") ;
fDispersion = 2.0 ;
fCpvEmcDistance = 3.0 ;
+ fTimeGate = 2.e-9 ;
fHeaderFileName = GetTitle() ;
fTrackSegmentsTitle = GetName() ;
fRecPointsTitle = GetName() ;
fRecParticlesTitle = GetName() ;
- fIDOptions = "" ;
-
+ fIDOptions = "dis time" ;
+
TString tempo(GetName()) ;
tempo.Append(":") ;
tempo.Append(Version()) ;
SetName(tempo) ;
fRecParticlesInRun = 0 ;
-
+
Init() ;
}
//____________________________________________________________________________
AliPHOSPIDv1::~AliPHOSPIDv1()
{
- //dtor
+ if(fTrackSegments) fTrackSegments->Delete() ;
+ if(fEmcRecPoints) fEmcRecPoints->Delete() ;
+ if(fCpvRecPoints) fCpvRecPoints->Delete() ;
+ if(fRecParticles) fRecParticles->Delete() ;
}
-//____________________________________________________________________________
-Bool_t AliPHOSPIDv1::ReadTrackSegments(Int_t event)
-{
- // Reads TrackSegments an extracts the title of the RecPoints
- // branch from which TS were made of.
- // Then reads both TrackSegments and RecPoints.
-
- //Fist read Track Segment Branch and extract RecPointsBranch from fTSMaker
-
- // Get TreeR header from file
- if(gAlice->TreeR()==0){
- cerr << "ERROR: AliPHOSPIDv1::ReadTrackSegments -> There is no Reconstruction Tree" << endl;
- return kFALSE;
- }
-
- // Find TrackSegments
- TBranch * tsbranch = 0;
- TBranch * tsmakerbranch = 0;
- TObjArray * lob = (TObjArray*)gAlice->TreeR()->GetListOfBranches() ;
- TIter next(lob) ;
- TBranch * branch = 0 ;
- Bool_t phostsfound = kFALSE, tsmakerfound = kFALSE ;
-
- TString taskName(GetName()) ;
- taskName.ReplaceAll(Version(), "") ;
-
- while ( (branch = (TBranch*)next()) && (!phostsfound || !tsmakerfound) ) {
- if ( (strcmp(branch->GetName(), "PHOSTS")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) ) {
- phostsfound = kTRUE ;
- tsbranch = branch ;
-
- } else if ( (strcmp(branch->GetName(), "AliPHOSTrackSegmentMaker")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) ) {
- tsmakerfound = kTRUE ;
- tsmakerbranch = branch ;
- }
- }
- if ( !phostsfound || !tsmakerfound ) {
- cerr << "WARNING: AliPHOSPIDv1::ReadTrackSegments -> TrackSegments and/or TrackSegmentMaker branch with name " << taskName.Data()
- << " not found" << endl ;
- return kFALSE ;
- }
-
- AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
-
- TClonesArray * trackSegments = gime->TrackSegments() ;
- trackSegments->Clear() ;
- tsbranch->SetAddress(&trackSegments) ;
-
- AliPHOSTrackSegmentMaker * tsmaker = 0 ;
- tsmakerbranch->SetAddress(&tsmaker) ;
- tsmakerbranch->GetEntry(0) ;
- TString tsmakerName( fRecParticlesTitle ) ;
- tsmakerName.Append(tsmaker->Version()) ;
- tsmaker = gime->TrackSegmentMaker(tsmakerName) ;
-
- tsbranch->GetEntry(0) ;
- tsmakerbranch->GetEntry(0) ;
-
- fRecPointsTitle = tsmaker->GetRecPointsBranch() ;
-
- // Find RecPoints
- TBranch * emcbranch = 0;
- TBranch * cpvbranch = 0;
- TBranch * clusterizerbranch = 0;
- lob = (TObjArray*)gAlice->TreeR()->GetListOfBranches() ;
- TIter next2(lob) ;
- branch = 0 ;
- Bool_t phosemcfound = kFALSE, phoscpvfound = kFALSE, clusterizerfound = kFALSE ;
-
- while ( (branch = (TBranch*)next2()) && (!phosemcfound || !phoscpvfound || !clusterizerfound) ) {
- if ( (strcmp(branch->GetName(), "PHOSEmcRP")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) ) {
- phosemcfound = kTRUE ;
- emcbranch = branch ;
- }
-
- else if ( (strcmp(branch->GetName(), "PHOSCpvRP")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) ) {
- phoscpvfound = kTRUE ;
- cpvbranch = branch ;
-
- } else if ( (strcmp(branch->GetName(), "AliPHOSClusterizer")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) ) {
- clusterizerfound = kTRUE ;
- clusterizerbranch = branch ;
- }
- }
- if ( !phosemcfound || !phoscpvfound || !clusterizerfound ) {
- cerr << "WARNING: AliPHOSTrackPIDv1::ReadTrackSegments -> emc(cpv)RecPoints and/or Clusterizer branch with name " << taskName.Data()
- << " not found" << endl ;
- return kFALSE ;
- }
-
- TObjArray * emcRecPoints = gime->EmcRecPoints() ;
- emcRecPoints->Clear() ;
- emcbranch->SetAddress(&emcRecPoints) ;
-
- TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
- cpvRecPoints->Clear() ;
- cpvbranch->SetAddress(&cpvRecPoints) ;
-
-
- AliPHOSClusterizer * clusterizer = 0 ;
- clusterizerbranch->SetAddress(&clusterizer) ;
- clusterizerbranch->GetEntry(0) ;
- TString clusterizerName( fTrackSegmentsTitle ) ;
- clusterizerName.Append(clusterizer->Version()) ;
- clusterizer = gime->Clusterizer(clusterizerName) ;
-
- emcbranch->GetEntry(0) ;
- cpvbranch->GetEntry(0) ;
- clusterizerbranch->GetEntry(0) ;
-
- return kTRUE ;
-
-}
//____________________________________________________________________________
Float_t AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cpv, Option_t * Axis)const
void AliPHOSPIDv1::Exec(Option_t * option)
{
//Steering method
-
+
if( strcmp(GetName(), "")== 0 )
Init() ;
Print("") ;
return ;
}
-
+
gAlice->GetEvent(0) ;
//check, if the branch with name of this" already exits?
TObjArray * lob = (TObjArray*)gAlice->TreeR()->GetListOfBranches() ;
TString taskName(GetName()) ;
taskName.Remove(taskName.Index(Version())-1) ;
-
+
while ( (branch = (TBranch*)next()) && (!phospidfound || !pidfound) ) {
if ( (strcmp(branch->GetName(), "PHOSPID")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) )
phospidfound = kTRUE ;
else if ( (strcmp(branch->GetName(), "AliPHOSPID")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) )
pidfound = kTRUE ;
}
-
+
if ( phospidfound || pidfound ) {
cerr << "WARNING: AliPHOSPIDv1::Exec -> RecParticles and/or PIDtMaker branch with name "
<< taskName.Data() << " already exits" << endl ;
AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
for(ievent = 0; ievent < nevents; ievent++){
-
gime->Event(ievent,"R") ;
-
- // if(!ReadTrackSegments(ievent)) //reads TrackSegments for event ievent
- // continue ;
-
+
MakeRecParticles() ;
-
+
WriteRecParticles(ievent);
-
+
if(strstr(option,"deb"))
PrintRecParticles(option) ;
}
-
+
if(strstr(option,"tim")){
gBenchmark->Stop("PHOSPID");
cout << "AliPHOSPID:" << endl ;
<< gBenchmark->GetCpuTime("PHOSPID")/nevents << " seconds per event " << endl ;
cout << endl ;
}
-
+
}
//____________________________________________________________________________
void AliPHOSPIDv1::Init()
return ;
}
- //add Task to //YSAlice/tasks/Reconstructioner/PHOS
- // TTask * aliceRe = (TTask*)gROOT->FindObjectAny("YSAlice/tasks/Reconstructioner") ;
- // TTask * phosRe = (TTask*)aliceRe->GetListOfTasks()->FindObject("PHOS") ;
- // phosRe->Add(this) ;
gime->PostPID(this) ;
// create a folder on the white board //YSAlice/WhiteBoard/RecParticles/PHOS/recparticlesName
gime->PostRecParticles(taskName.Data() ) ;
Bool_t ellips = fIDOptions.Contains("ell",TString::kIgnoreCase ) ;
Bool_t disp = fIDOptions.Contains("dis",TString::kIgnoreCase ) ;
+ Bool_t time = fIDOptions.Contains("tim",TString::kIgnoreCase ) ;
while ( (ts = (AliPHOSTrackSegment *)next()) ) {
new( (*recParticles)[index] ) AliPHOSRecParticle() ;
rp = (AliPHOSRecParticle *)recParticles->At(index) ;
rp->SetTraskSegment(index) ;
+ rp->SetIndexInList(index) ;
AliPHOSEmcRecPoint * emc = 0 ;
if(ts->GetEmcIndex()>=0)
if(ts->GetCpvIndex()>=0)
cpv = (AliPHOSRecPoint *) cpvRecPoints->At(ts->GetCpvIndex()) ;
- AliPHOSRecPoint * ppsd = 0 ;
- if(ts->GetPpsdIndex()>=0)
- ppsd= (AliPHOSRecPoint *) cpvRecPoints->At(ts->GetPpsdIndex()) ;
-
//set momentum and energy first
Float_t e = emc->GetEnergy() ;
- TVector3 dir = GetMomentumDirection(emc,cpv,ppsd) ;
+ TVector3 dir = GetMomentumDirection(emc,cpv) ;
dir.SetMag(e) ;
rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),e) ;
if(emc->GetDispersion() > fDispersion )
showerprofile = 1 ; // not narrow
-
- // Looking at the photon conversion detector
- Int_t pcdetector= 0 ; //1 hit and 0 no hit
- if(ppsd)
- if(GetDistance(emc, ppsd, "R") < fCpvEmcDistance)
- pcdetector = 1 ;
-
+ Int_t slow = 0 ;
+ if(time)
+ if(emc->GetTime() > fTimeGate )
+ slow = 0 ;
+
// Looking at the CPV detector
Int_t cpvdetector= 0 ; //1 hit and 0 no hit
if(cpv)
if(GetDistance(emc, cpv, "R") < fCpvEmcDistance)
cpvdetector = 1 ;
- Int_t type = showerprofile + 2 * pcdetector + 4 * cpvdetector ;
- rp->SetType(type) ;
- rp->SetIndexInList(index) ;
+ Int_t type = showerprofile + 2 * slow + 4 * cpvdetector ;
+ rp->SetType(type) ;
index++ ;
}
cout << " TrackSegments Branch title: " << fTrackSegmentsTitle.Data() << endl ;
cout << " RecParticles Branch title " << fRecParticlesTitle.Data() << endl;
cout << "with parameters: " << endl ;
- cout << " Maximal EMC - CPV (PPSD) distance (cm) " << fCpvEmcDistance << endl ;
+ cout << " Maximal EMC - CPV distance (cm) " << fCpvEmcDistance << endl ;
if(fIDOptions.Contains("dis",TString::kIgnoreCase ))
- cout << " dispersion cut " << fDispersion << endl ;
+ cout << " dispersion cut " << fDispersion << endl ;
if(fIDOptions.Contains("ell",TString::kIgnoreCase )){
- cout << "Eliptic cuts function: " << endl ;
+ cout << " Eliptic cuts function: " << endl ;
cout << fFormula->GetTitle() << endl ;
}
+ if(fIDOptions.Contains("tim",TString::kIgnoreCase ))
+ cout << " Time Gate uzed: " << fTimeGate << endl ;
cout << "============================================" << endl ;
}
taskName.Remove(taskName.Index(Version())-1) ;
TClonesArray * recParticles = gime->RecParticles(taskName) ;
recParticles->Expand(recParticles->GetEntriesFast() ) ;
-
+
//Make branch in TreeR for RecParticles
char * filename = 0;
if(gSystem->Getenv("CONFIG_SPLIT_FILE")!=0){ //generating file name
}
cwd->cd();
}
-
+
//second, pid
Int_t splitlevel = 0 ;
AliPHOSPIDv1 * pid = this ;
TBranch * pidBranch = gAlice->TreeR()->Branch("AliPHOSPID","AliPHOSPIDv1",&pid,bufferSize,splitlevel);
- pidBranch->SetTitle(fRecParticlesTitle);
+ pidBranch->SetTitle(fRecParticlesTitle.Data());
if (filename) {
pidBranch->SetFile(filename);
TIter next( pidBranch->GetListOfBranches());
rpBranch->Fill() ;
pidBranch->Fill() ;
-
+
gAlice->TreeR()->Write(0,kOverwrite) ;
}
}
//____________________________________________________________________________
-TVector3 AliPHOSPIDv1::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * cpv,AliPHOSRecPoint * ppsd)const
+TVector3 AliPHOSPIDv1::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * cpv)const
{
// Calculates the momentum direction:
// 1. if only a EMC RecPoint, direction is given by IP and this RecPoint
- // 2. if a EMC RecPoint and one PPSD RecPoint, direction is given by the line through the 2 recpoints
- // 3. if a EMC RecPoint and two PPSD RecPoints, dirrection is given by the average line through
- // the 2 pairs of recpoints
- // However because of the poor position resolution of PPSD the direction is always taken as if we were
+ // 2. if a EMC RecPoint and CPV RecPoint, direction is given by the line through the 2 recpoints
+ // However because of the poor position resolution of PPSD the direction is always taken as if we were
// in case 1.
TVector3 dir(0,0,0) ;
TString taskName(GetName()) ;
taskName.Remove(taskName.Index(Version())-1) ;
TClonesArray * recParticles = gime->RecParticles(taskName) ;
-
-
+
cout << "AliPHOSPIDv1: event "<<gAlice->GetEvNumber() << endl ;
cout << " found " << recParticles->GetEntriesFast() << " RecParticles " << endl ;
fRecParticlesInRun += recParticles->GetEntriesFast() ;
-
+
if(strstr(option,"all")) { // printing found TS
cout << " PARTICLE "
Text_t particle[11];
switch(rp->GetType()) {
- case AliPHOSFastRecParticle::kNEUTRALEM:
- strcpy( particle, "NEUTRAL_EM");
+ case AliPHOSFastRecParticle::kNEUTRALEMFAST:
+ strcpy( particle, "NEUTRAL EM FAST");
break;
- case AliPHOSFastRecParticle::kNEUTRALHA:
- strcpy(particle, "NEUTRAL_HA");
+ case AliPHOSFastRecParticle::kNEUTRALHAFAST:
+ strcpy(particle, "NEUTRAL HA FAST");
break;
- case AliPHOSFastRecParticle::kGAMMA:
- strcpy(particle, "GAMMA");
+ case AliPHOSFastRecParticle::kNEUTRALEMSLOW:
+ strcpy(particle, "NEUTRAL EM SLOW");
break ;
- case AliPHOSFastRecParticle::kGAMMAHA:
- strcpy(particle, "GAMMA_H");
+ case AliPHOSFastRecParticle::kNEUTRALHASLOW:
+ strcpy(particle, "NEUTRAL HA SLOW");
break ;
- case AliPHOSFastRecParticle::kABSURDEM:
- strcpy(particle, "ABSURD_EM") ;
+ case AliPHOSFastRecParticle::kCHARGEDEMFAST:
+ strcpy(particle, "CHARGED EM FAST") ;
break ;
- case AliPHOSFastRecParticle::kABSURDHA:
- strcpy(particle, "ABSURD_HA") ;
+ case AliPHOSFastRecParticle::kCHARGEDHAFAST:
+ strcpy(particle, "CHARGED HA FAST") ;
break ;
- case AliPHOSFastRecParticle::kELECTRON:
- strcpy(particle, "ELECTRON") ;
+ case AliPHOSFastRecParticle::kCHARGEDEMSLOW:
+ strcpy(particle, "CHARGEDEMSLOW") ;
break ;
- case AliPHOSFastRecParticle::kCHARGEDHA:
- strcpy(particle, "CHARGED_HA") ;
+ case AliPHOSFastRecParticle::kCHARGEDHASLOW:
+ strcpy(particle, "CHARGED HA SLOW") ;
break ;
}
// Int_t nprimaries;
// primaries = rp->GetPrimaries(nprimaries);
- cout << setw(15) << particle << " "
- << setw(3) << rp->GetIndexInList() << " " ;
+ cout << setw(10) << particle << " "
+ << setw(5) << rp->GetIndexInList() << " " ;
// << setw(4) << nprimaries << " ";
// for (Int_t iprimary=0; iprimary<nprimaries; iprimary++)
// cout << setw(4) << primaries[iprimary] << " ";
virtual char * GetRecParticlesBranch()const {return (char*) fRecParticlesTitle.Data() ;}
virtual char * GetTrackSegmentsBranch()const{return (char*) fTrackSegmentsTitle.Data(); }
virtual const Int_t GetRecParticlesInRun() const {return fRecParticlesInRun ;}
-
+
virtual void Init() ;
virtual void PlotDispersionCuts()const ;
virtual void Print(Option_t * option)const ;
virtual void SetIdentificationMethod(char * option = "CPV DISP" ){fIDOptions = option ;}
virtual void SetShowerProfileCut(char * formula = "0.35*0.35 - (x-1.386)*(x-1.386) - 1.707*1.707*(y-1.008)*(y-1.008)") ;
virtual void SetDispersionCut(Float_t cut){fDispersion = cut ; }
- virtual void SetCpvtoEmcDistanceCut(Float_t cut ) {fCpvEmcDistance = cut ;}
+ virtual void SetCpvtoEmcDistanceCut(Float_t cut ) {fCpvEmcDistance = cut ;}
+ virtual void SetTimeGate(Float_t gate) {fTimeGate = gate ;}
virtual void SetTrackSegmentsBranch(const char* title) { fTrackSegmentsTitle = title;}
virtual void SetRecParticlesBranch (const char* title) { fRecParticlesTitle = title;}
virtual const char * Version() const { return "pid-v1" ; }
private:
void MakeRecParticles(void ) ;
- Float_t GetDistance(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * cpv, Option_t * Axis)const ; // Relative Distance PPSD-EMC
- TVector3 GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * cpv, AliPHOSRecPoint * ppsd)const ;
+ Float_t GetDistance(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * cpv, Option_t * Axis)const ; // Relative Distance CPV-EMC
+ TVector3 GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * cpv)const ;
void PrintRecParticles(Option_t * option) ;
- virtual Bool_t ReadTrackSegments(Int_t event) ;
virtual void WriteRecParticles(Int_t event) ;
private:
TFormula * fFormula ; // formula to define cut on the shouer elips axis
Float_t fDispersion ; // dispersion cut
Float_t fCpvEmcDistance ; // Max EMC-CPV distance
+ Float_t fTimeGate ; // Time of the latest EmcRecPoint accepted as EM
Int_t fRecParticlesInRun ; //! Total number of recparticles in one run
+
ClassDef( AliPHOSPIDv1,1) // Particle identifier implementation version 1
};
TTask * aliceQA = (TTask*)gROOT->FindObjectAny("YSAlice/tasks/QA") ;
TTask * phosQA = (TTask*)aliceQA->GetListOfTasks()->FindObject("PHOS") ;
if (phosQA) // PHOS QA Tasks container exists
- phosQA->Add(this) ;
- else // create //YSAlice/tasks/QA/PHOS
- aliceQA->Add(this) ;
+ phosQA->Add(this) ;
+ else // create //YSAlice/tasks/QA/PHOS
+ aliceQA->Add(this) ;
fCheckablesList = new TList() ;
}
public:
- AliPHOSQAChecker(){
- fCheckable = 0 ;
- fCheckablesList = 0 ;
- } ; // default ctor (not to be used)
+ AliPHOSQAChecker(){} ; // default ctor (not to be used)
AliPHOSQAChecker(const char * name, const char * title) ; // ctor
AliPHOSQAChecker(AliPHOSQAChecker& obj) {assert(0==1);}
virtual ~AliPHOSQAChecker() ; // dtor
void StatusAll() { ExecuteTask("S") ; }
friend void AliPHOSQAVirtualCheckable::AddChecker(AliPHOSQAChecker * ch) ;
- // friend AliPHOSQAVirtualCheckable::AliPHOSQAVirtualCheckable(const char * name) ;
+ friend AliPHOSQAVirtualCheckable::AliPHOSQAVirtualCheckable(const char * name) ;
private:
{
// print the name
- TString status = "inactive" ;
- if ( IsActive() )
- status = "active" ;
- cout << "Checker : " << GetName() << " : " << GetTitle() << " : Mean = " << fMean << " Rms = " << fRms
- << " : Status = " << status << endl ;
+ cout << "Checker : " << GetName() << " : " << GetTitle() << " : Mean = " << fMean << " Rms = " << fRms << endl ;
}
public:
- AliPHOSQAVirtualCheckable(){
- fAlarms = 0 ;
- fChecker = 0 ;
- } // default ctor not to be used
+ AliPHOSQAVirtualCheckable(){} // default ctor not to be used
AliPHOSQAVirtualCheckable(const char * name) ; // ctor
AliPHOSQAVirtualCheckable(AliPHOSQAVirtualCheckable& obj) {assert(0==1);}
virtual ~AliPHOSQAVirtualCheckable() ; // dtor
for ( i = 0 ; i < hits->GetEntries() ; i++ ) {
AliPHOSHit * hit = (AliPHOSHit *) hits->At(i) ;
// Assign primary number only if contribution is significant
-
+
if( hit->GetEnergy() > fPrimThreshold)
- new((*sdigits)[nSdigits]) AliPHOSDigit( hit->GetPrimary(), hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
+ new((*sdigits)[nSdigits]) AliPHOSDigit(hit->GetPrimary(),hit->GetId(),
+ Digitize(hit->GetEnergy()), hit->GetTime()) ;
else
- new((*sdigits)[nSdigits]) AliPHOSDigit( -1 , hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
-
- nSdigits++ ;
+ new((*sdigits)[nSdigits]) AliPHOSDigit( -1 , hit->GetId(),
+ Digitize(hit->GetEnergy()), hit->GetTime()) ;
+ nSdigits++ ;
- }
-
+ }
} // loop over tracks
sdigits->Sort() ;
nSdigits = sdigits->GetEntriesFast() ;
-
+
sdigits->Expand(nSdigits) ;
Int_t i ;
for (i = 0 ; i < nSdigits ; i++) {
//____________________________________________________________________________
AliPHOSTrackSegment::AliPHOSTrackSegment( AliPHOSEmcRecPoint * emc ,
- AliPHOSRecPoint * ppsdrp1,
- AliPHOSRecPoint * ppsdrp2 )
+ AliPHOSRecPoint * ppsdrp1)
{
// ctor
else
fPpsdUpRecPoint = -1 ;
- if( ppsdrp2 )
- fPpsdLowRecPoint = ppsdrp2->GetIndexInList() ;
- else
- fPpsdLowRecPoint = -1 ;
fIndexInList = -1 ;
}
TObject::Copy(obj) ;
( (AliPHOSTrackSegment &)obj ).fEmcRecPoint = fEmcRecPoint ;
- ( (AliPHOSTrackSegment &)obj ).fPpsdLowRecPoint = fPpsdLowRecPoint ;
( (AliPHOSTrackSegment &)obj ).fPpsdUpRecPoint = fPpsdUpRecPoint ;
( (AliPHOSTrackSegment &)obj ).fIndexInList = fIndexInList ;
}
else
cout << "No CPV RecPoint " << endl ;
- if(fPpsdLowRecPoint >= 0)
- cout << "PPSD RecPoint # " << fPpsdLowRecPoint << endl ;
- else
- cout << "No PPSD RecPoint " << endl ;
cout << "------------------------------------ " << endl ;
AliPHOSTrackSegment() {}
AliPHOSTrackSegment(AliPHOSEmcRecPoint * EmcRecPoint ,
- AliPHOSRecPoint * PpsdUp,
- AliPHOSRecPoint * PpsdLow ) ; // ctor
+ AliPHOSRecPoint * PpsdUp) ;
AliPHOSTrackSegment(const AliPHOSTrackSegment & ts) ; // ctor
virtual ~AliPHOSTrackSegment() { }
Int_t GetIndexInList() const { return fIndexInList ; }
Int_t GetEmcIndex() const { return fEmcRecPoint ; }
- Int_t GetPpsdIndex() const { return fPpsdLowRecPoint;}
Int_t GetCpvIndex() const { return fPpsdUpRecPoint; }
virtual void Print(Option_t * option) const;
Int_t fEmcRecPoint ; // The EMC reconstructed point index in array stored in TreeR/PHOSEmcRP
Int_t fIndexInList ; // the index of this TrackSegment in the list stored in TreeR (to be set by analysis)
- Int_t fPpsdLowRecPoint ; // The PPSD reconstructed point from the lower layer index in array stored in TreeR/PHOSPpsdRP
- Int_t fPpsdUpRecPoint ; // The PPSD reconstructed point from the upper layer index in array stored in TreeR/PHOSPpsdRP
+ Int_t fPpsdUpRecPoint ; // The CPV reconstructed point from the upper layer index in array stored in TreeR/PHOSPpsdRP
ClassDef(AliPHOSTrackSegment,1) // Track segment in PHOS
// virtual void Set... // method to choose recPoints: along z only, along x ...???
// virtual void SetChoosingAlgirithm() = 0 ;
// virtual void SetMaxEmcCpvDistance(Float_t r) = 0 ;
- virtual Bool_t ReadRecPoints(Int_t ievent) = 0 ;
virtual void SetRecPointsBranch(const char * title) = 0 ;
virtual void SetTrackSegmentsBranch(const char * title) = 0 ;
virtual const char * Version() const = 0 ;
#include "AliPHOSClusterizerv1.h"
#include "AliPHOSTrackSegment.h"
#include "AliPHOSCpvRecPoint.h"
-#include "AliPHOSPpsdRecPoint.h"
#include "AliPHOSLink.h"
#include "AliPHOSGetter.h"
#include "AliPHOS.h"
fEmcLast = 0 ;
fCpvFirst = 0 ;
fCpvLast = 0 ;
- fPpsdFirst = 0 ;
- fPpsdLast = 0 ;
- fLinkLowArray = 0 ;
fLinkUpArray = 0 ;
fHeaderFileName = "" ;
fRecPointsBranchTitle = "" ;
fEmcLast = 0 ;
fCpvFirst = 0 ;
fCpvLast = 0 ;
- fPpsdFirst = 0 ;
- fPpsdLast = 0 ;
fHeaderFileName = GetTitle() ;
fRecPointsBranchTitle = GetName() ;
fTrackSegmentsBranchTitle = GetName() ;
fTrackSegmentsInRun = 0 ;
- TString tempo(GetName()) ;
- tempo.Append(":") ;
- tempo.Append(Version()) ;
- SetName(tempo.Data()) ;
-
+ TString tsmName( GetName()) ;
+ tsmName.Append(":") ;
+ tsmName.Append(Version()) ;
+ SetName(tsmName) ;
+
Init() ;
}
AliPHOSTrackSegmentMakerv1::~AliPHOSTrackSegmentMakerv1()
{
// dtor
- if(fLinkLowArray) delete fLinkLowArray ;
if(fLinkUpArray) delete fLinkUpArray ;
}
{
// Finds first and last indexes between which
// clusters from one PHOS module are
-
TString taskName(GetName()) ;
taskName.Remove(taskName.Index(Version())-1) ;
AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
- const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
TObjArray * emcRecPoints = gime->EmcRecPoints(taskName) ;
TObjArray * cpvRecPoints = gime->CpvRecPoints(taskName) ;
-
+
//First EMC clusters
Int_t totalEmc = emcRecPoints->GetEntriesFast() ;
for(fEmcFirst = fEmcLast; (fEmcLast < totalEmc) &&
//Now CPV clusters
Int_t totalCpv = cpvRecPoints->GetEntriesFast() ;
- if(fModule <= geom->GetNCPVModules()){ // in CPV geometry
-
for(fCpvFirst = fCpvLast; (fCpvLast < totalCpv) &&
(((AliPHOSRecPoint *) cpvRecPoints->At(fCpvLast))->GetPHOSMod() == fModule );
fCpvLast ++) ;
-
- fPpsdFirst = fCpvLast ; //To avoid scanning RecPoints between fPpsdFirst and fPpsdLast
- fPpsdLast = fCpvLast ; //and to be ready to switch to mixed geometry
- }
- else{ //in PPSD geometry
- fCpvLast = fPpsdLast ;
- //Upper layer first
- for(fCpvFirst = fCpvLast; (fCpvLast < totalCpv) &&
- (((AliPHOSPpsdRecPoint *) cpvRecPoints->At(fCpvLast))->GetPHOSMod() == fModule ) &&
- (((AliPHOSPpsdRecPoint *) cpvRecPoints->At(fCpvLast))->GetUp()) ;
- fCpvLast ++) ;
-
- fPpsdLast= fCpvLast ;
- for(fPpsdFirst = fPpsdLast; (fPpsdLast < totalCpv) &&
- (((AliPHOSPpsdRecPoint *) cpvRecPoints->At(fPpsdLast))->GetPHOSMod() == fModule ) &&
- (!((AliPHOSPpsdRecPoint *) cpvRecPoints->At(fPpsdLast))->GetUp()) ;
- fPpsdLast ++) ;
- }
-
+
}
//____________________________________________________________________________
if ( strcmp(GetTitle(), "") == 0 )
SetTitle("galice.root") ;
-
- TString taskName(GetName()) ;
- taskName.Remove(taskName.Index(Version())-1) ;
-
- AliPHOSGetter * gime = AliPHOSGetter::GetInstance(GetTitle(), taskName) ;
+
+ TString branchname = GetName() ;
+ branchname.Remove(branchname.Index(Version())-1) ;
+
+ AliPHOSGetter * gime = AliPHOSGetter::GetInstance(GetTitle(), branchname) ;
if ( gime == 0 ) {
cerr << "ERROR: AliPHOSTrackSegmentMakerv1::Init -> Could not obtain the Getter object !" << endl ;
return ;
}
- fLinkLowArray = new TClonesArray("AliPHOSLink", 1000);
fLinkUpArray = new TClonesArray("AliPHOSLink", 1000);
//add Task to //YSAlice/tasks/Reconstructioner/PHOS
gime->PostTrackSegmentMaker(this) ;
// create a folder on the white board //YSAlice/WhiteBoard/RecPoints/PHOS/trackSegmentsName
- gime->PostTrackSegments(taskName) ;
+ gime->PostTrackSegments(branchname) ;
}
// Finds distances (links) between all EMC and PPSD clusters,
// which are not further apart from each other than fR0
// and sort them in accordance with this distance
-
+
TString taskName(GetName()) ;
taskName.Remove(taskName.Index(Version())-1) ;
-
AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
TObjArray * emcRecPoints = gime->EmcRecPoints(taskName) ;
TObjArray * cpvRecPoints = gime->CpvRecPoints(taskName) ;
fLinkUpArray->Clear() ;
- fLinkLowArray->Clear() ;
- AliPHOSRecPoint * ppsd ;
AliPHOSRecPoint * cpv ;
AliPHOSEmcRecPoint * emcclu ;
- Int_t iLinkLow = 0 ;
Int_t iLinkUp = 0 ;
Int_t iEmcRP;
for(iEmcRP = fEmcFirst; iEmcRP < fEmcLast; iEmcRP++ ) {
emcclu = (AliPHOSEmcRecPoint *) emcRecPoints->At(iEmcRP) ;
- Bool_t toofar ;
- Int_t iPpsd ;
- for(iPpsd = fPpsdFirst; iPpsd < fPpsdLast;iPpsd++ ) {
-
- ppsd = (AliPHOSRecPoint *) cpvRecPoints->At(iPpsd) ;
- Float_t r = GetDistanceInPHOSPlane(emcclu, ppsd, toofar) ;
-
- if(toofar)
- break ;
- if(r < fR0)
- new ((*fLinkLowArray)[iLinkLow++]) AliPHOSLink(r, iEmcRP, iPpsd) ;
- }
-
+ Bool_t toofar ;
Int_t iCpv = 0 ;
for(iCpv = fCpvFirst; iCpv < fCpvLast;iCpv++ ) {
}
}
- fLinkLowArray->Sort() ; //first links with smallest distances
- fLinkUpArray->Sort() ;
+ fLinkUpArray->Sort() ; //first links with smallest distances
}
//____________________________________________________________________________
TObjArray * emcRecPoints = gime->EmcRecPoints(taskName) ;
TObjArray * cpvRecPoints = gime->CpvRecPoints(taskName) ;
TClonesArray * trackSegments = gime->TrackSegments(taskName) ;
-
+
//Make arrays to mark clusters already chosen
Int_t * emcExist = 0;
if(fEmcLast > fEmcFirst)
for(index = 0; index <fCpvLast-fCpvFirst; index ++)
cpvExist[index] = kTRUE ;
- Bool_t * ppsdExist = 0;
- if(fPpsdLast > fPpsdFirst)
- ppsdExist = new Bool_t[fPpsdLast-fPpsdFirst] ;
- for(index = 0; index <fPpsdLast-fPpsdFirst; index ++)
- ppsdExist[index] = kTRUE ;
// Finds the smallest links and makes pairs of CPV and EMC clusters with smallest distance
- TIter nextLow(fLinkLowArray) ;
TIter nextUp(fLinkUpArray) ;
- AliPHOSLink * linkLow ;
AliPHOSLink * linkUp ;
-
-
+
AliPHOSRecPoint * nullpointer = 0 ;
-
- while ( (linkLow = (AliPHOSLink *)nextLow() ) ){
-
- if( (emcExist[linkLow->GetEmc()-fEmcFirst]> 0) &&
- ppsdExist[linkLow->GetPpsd()-fPpsdFirst] ){ // RecPoints not removed yet
- new ((*trackSegments)[fNTrackSegments]) AliPHOSTrackSegment((AliPHOSEmcRecPoint *) emcRecPoints->At(linkLow->GetEmc()),
- nullpointer,
- (AliPHOSRecPoint *)cpvRecPoints->At(linkLow->GetPpsd()) ) ;
-
- ((AliPHOSTrackSegment* )trackSegments->At(fNTrackSegments))->SetIndexInList(fNTrackSegments);
- //replace index of emc to negative and shifted index of TS
- emcExist[linkLow->GetEmc()-fEmcFirst] = -2 - fNTrackSegments ;
- //mark ppsd as used
- ppsdExist[linkLow->GetPpsd()-fPpsdFirst] = kFALSE ;
- fNTrackSegments++ ;
- }
- }
-
-
+
while ( (linkUp = (AliPHOSLink *)nextUp() ) ){
+
if(emcExist[linkUp->GetEmc()-fEmcFirst] != -1){ //without ppsd Up yet
if(cpvExist[linkUp->GetPpsd()-fCpvFirst]){ //CPV still exist
- if(emcExist[linkUp->GetEmc()-fEmcFirst] > 0){ //without ppsd Low => create new TS
-
- new ((* trackSegments)[fNTrackSegments]) AliPHOSTrackSegment((AliPHOSEmcRecPoint *) emcRecPoints->At(linkUp->GetEmc()) ,
- (AliPHOSRecPoint *)cpvRecPoints->At(linkUp->GetPpsd()),
- nullpointer) ;
- ((AliPHOSTrackSegment *) trackSegments->At(fNTrackSegments))->SetIndexInList(fNTrackSegments);
- fNTrackSegments++ ;
- }
- else{ // append ppsd Up to existing TS
- ((AliPHOSTrackSegment *)trackSegments->At(-2-emcExist[linkUp->GetEmc()-fEmcFirst]))->SetCpvRecPoint((AliPHOSCpvRecPoint *)cpvRecPoints->At(linkUp->GetPpsd()));
- }
-
+ new ((* trackSegments)[fNTrackSegments])
+ AliPHOSTrackSegment((AliPHOSEmcRecPoint *) emcRecPoints->At(linkUp->GetEmc()) ,
+ (AliPHOSRecPoint *)cpvRecPoints->At(linkUp->GetPpsd())) ;
+ ((AliPHOSTrackSegment *) trackSegments->At(fNTrackSegments))->SetIndexInList(fNTrackSegments);
+ fNTrackSegments++ ;
+
emcExist[linkUp->GetEmc()-fEmcFirst] = -1 ; //Mark emc that Cpv was found
//mark CPV recpoint as already used
- cpvExist[linkUp->GetPpsd()-fCpvFirst] = kFALSE ;
+ cpvExist[linkUp->GetPpsd()-fCpvFirst] = kFALSE ;
} //if ppsdUp still exist
}
}
Int_t iEmcRP ;
for(iEmcRP = 0; iEmcRP < fEmcLast-fEmcFirst ; iEmcRP++ ){
if(emcExist[iEmcRP] > 0 ){
- new ((*trackSegments)[fNTrackSegments]) AliPHOSTrackSegment((AliPHOSEmcRecPoint *)emcRecPoints->At(iEmcRP+fEmcFirst),
- nullpointer,
- nullpointer ) ;
+ new ((*trackSegments)[fNTrackSegments])
+ AliPHOSTrackSegment((AliPHOSEmcRecPoint *)emcRecPoints->At(iEmcRP+fEmcFirst),
+ nullpointer) ;
((AliPHOSTrackSegment *) trackSegments->At(fNTrackSegments))->SetIndexInList(fNTrackSegments);
fNTrackSegments++;
}
return ;
}
- //check, if the branch with name of this" already exits?
gAlice->GetEvent(0) ;
+ //check, if the branch with name of this" already exits?
TObjArray * lob = (TObjArray*)gAlice->TreeR()->GetListOfBranches() ;
TIter next(lob) ;
TBranch * branch = 0 ;
Bool_t phostsfound = kFALSE, tracksegmentmakerfound = kFALSE ;
- TString taskName(GetName()) ;
- taskName.Remove(taskName.Index(Version())-1) ;
+ TString branchname = GetName() ;
+ branchname.Remove(branchname.Index(Version())-1) ;
while ( (branch = (TBranch*)next()) && (!phostsfound || !tracksegmentmakerfound) ) {
- if ( (strcmp(branch->GetName(), "PHOSTS")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) )
+ if ( (strcmp(branch->GetName(), "PHOSTS")==0) && (strcmp(branch->GetTitle(), branchname.Data())==0) )
phostsfound = kTRUE ;
- else if ( (strcmp(branch->GetName(), "AliPHOSTrackSegmentMaker")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) )
+ else if ( (strcmp(branch->GetName(), "AliPHOSTrackSegmentMaker")==0) && (strcmp(branch->GetTitle(), GetName())==0) )
tracksegmentmakerfound = kTRUE ;
}
if ( phostsfound || tracksegmentmakerfound ) {
cerr << "WARNING: AliPHOSTrackSegmentMakerv1::Exec -> TrackSegments and/or TrackSegmentMaker branch with name "
- << taskName.Data() << " already exits" << endl ;
+ << branchname.Data() << " already exits" << endl ;
return ;
}
- AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
+ AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
Int_t nevents = (Int_t) gAlice->TreeE()->GetEntries() ;
Int_t ievent ;
for(ievent = 0; ievent < nevents; ievent++){
- // if(!ReadRecPoints(ievent)) //reads RecPoints for event ievent
- // continue;
gime->Event(ievent,"R") ;
-
//Make some initializations
fNTrackSegments = 0 ;
fEmcFirst = 0 ;
fEmcLast = 0 ;
fCpvFirst = 0 ;
fCpvLast = 0 ;
- fPpsdFirst= 0 ;
- fPpsdLast = 0 ;
- gime->TrackSegments(taskName)->Clear() ;
+ gime->TrackSegments(branchname)->Clear() ;
+ // if(!ReadRecPoints(ievent)) continue; //reads RecPoints for event ievent
+
for(fModule = 1; fModule <= geom->GetNModules() ; fModule++ ){
FillOneModule() ;
else
cout << "AliPHOSTrackSegmentMakerv1 not initialized " << endl ;
}
-//____________________________________________________________________________
-Bool_t AliPHOSTrackSegmentMakerv1::ReadRecPoints(Int_t event)
-{
- // Reads Emc and CPV recPoints
- // made previously with Clusterizer.
-
-
- //Make some initializations
-
- fNTrackSegments = 0 ;
- fEmcFirst = 0 ;
- fEmcLast = 0 ;
- fCpvFirst = 0 ;
- fCpvLast = 0 ;
- fPpsdFirst= 0 ;
- fPpsdLast = 0 ;
-
- // Get TreeR header from file
- if(gAlice->TreeR()==0){
- cerr << "ERROR: AliPHOSTrackSegmentMakerv1::ReadRecPoints -> There is no Reconstruction Tree" << endl;
- return kFALSE;
- }
-
-
- // Find RecPoints
- TBranch * emcbranch = 0;
- TBranch * cpvbranch = 0;
- TBranch * clusterizerbranch = 0;
- TObjArray * lob = (TObjArray*)gAlice->TreeR()->GetListOfBranches() ;
- TIter next(lob) ;
- TBranch * branch = 0 ;
- Bool_t phosemcfound = kFALSE, phoscpvfound = kFALSE, clusterizerfound = kFALSE ;
-
- TString taskName(GetName()) ;
- taskName.ReplaceAll(Version(), "") ;
-
- while ( (branch = (TBranch*)next()) && (!phosemcfound || !phoscpvfound || !clusterizerfound) ) {
- if ( (strcmp(branch->GetName(), "PHOSEmcRP")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) ) {
- phosemcfound = kTRUE ;
- emcbranch = branch ;
- }
-
- else if ( (strcmp(branch->GetName(), "PHOSCpvRP")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) ) {
- phoscpvfound = kTRUE ;
- cpvbranch = branch ;
-
- } else if ( (strcmp(branch->GetName(), "AliPHOSClusterizer")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) ) {
- clusterizerfound = kTRUE ;
- clusterizerbranch = branch ;
- }
- }
- if ( !phosemcfound || !phoscpvfound || !clusterizerfound ) {
- cerr << "WARNING: AliPHOSTrackSegmentMakerv1::ReadRecPoints -> emc(cpv)RecPoints and/or Clusterizer branch with name " << taskName.Data()
- << " not found" << endl ;
- return kFALSE ;
- }
-
- AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
-
- TObjArray * emcRecPoints = gime->EmcRecPoints() ;
- emcRecPoints->Clear() ;
- emcbranch->SetAddress(&emcRecPoints) ;
-
- TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
- cpvRecPoints->Clear() ;
- cpvbranch->SetAddress(&cpvRecPoints) ;
-
- TClonesArray * trackSegments = gime->TrackSegments() ;
- trackSegments->Clear() ;
-
- AliPHOSClusterizer * clusterizer = 0 ;
- clusterizerbranch->SetAddress(&clusterizer) ;
- clusterizerbranch->GetEntry(0) ;
- TString clusterizerName( fTrackSegmentsBranchTitle ) ;
- clusterizerName.Append(clusterizer->Version()) ;
- clusterizer = gime->Clusterizer(clusterizerName) ;
-
- emcbranch->GetEntry(0) ;
- cpvbranch->GetEntry(0) ;
-
- clusterizerbranch->GetEntry(0) ;
-
- printf("ReadRecPoint: EMC=%d, CPV=%d\n",emcRecPoints->GetEntries(),cpvRecPoints->GetEntries());
-
- return kTRUE ;
-
-}
//____________________________________________________________________________
void AliPHOSTrackSegmentMakerv1::WriteTrackSegments(Int_t event)
// If yes - exits without writing.
AliPHOSGetter *gime = AliPHOSGetter::GetInstance() ;
- TString taskName(GetName()) ;
- taskName.Remove(taskName.Index(Version())-1) ;
- TClonesArray * trackSegments = gime->TrackSegments(taskName) ;
+ TString branchName(GetName() ) ;
+ branchName.Remove(branchName.Index(Version())-1) ;
+
+ TClonesArray * trackSegments = gime->TrackSegments(branchName) ;
trackSegments->Expand(trackSegments->GetEntriesFast()) ;
//Make branch in TreeR for TrackSegments
//First TS
Int_t bufferSize = 32000 ;
TBranch * tsBranch = gAlice->TreeR()->Branch("PHOSTS",&trackSegments,bufferSize);
- tsBranch->SetTitle(fTrackSegmentsBranchTitle);
+ tsBranch->SetTitle(branchName);
if (filename) {
tsBranch->SetFile(filename);
TIter next( tsBranch->GetListOfBranches());
AliPHOSTrackSegmentMakerv1 * ts = this ;
TBranch * tsMakerBranch = gAlice->TreeR()->Branch("AliPHOSTrackSegmentMaker","AliPHOSTrackSegmentMakerv1",
&ts,bufferSize,splitlevel);
- tsMakerBranch->SetTitle(fTrackSegmentsBranchTitle);
+ tsMakerBranch->SetTitle(branchName);
if (filename) {
tsMakerBranch->SetFile(filename);
TIter next( tsMakerBranch->GetListOfBranches());
taskName.Remove(taskName.Index(Version())-1) ;
TClonesArray * trackSegments = AliPHOSGetter::GetInstance()->TrackSegments(taskName) ;
+
cout << "AliPHOSTrackSegmentMakerv1: event "<<gAlice->GetEvNumber() << endl ;
cout << " Found " << trackSegments->GetEntriesFast() << " trackSegments " << endl ;
fTrackSegmentsInRun += trackSegments->GetEntriesFast() ;
if(strstr(option,"all")) { // printing found TS
- cout << "TrackSegment # " << " EMC RP# " << " CPV RP# " << " PPSD RP#" << endl ;
+ cout << "TrackSegment # " << " EMC RP# " << " CPV RP# " << endl ;
Int_t index;
for (index = 0 ; index <trackSegments->GetEntriesFast() ; index++) {
AliPHOSTrackSegment * ts = (AliPHOSTrackSegment * )trackSegments->At(index) ;
cout<<" "<< setw(4) << ts->GetIndexInList() << " "
<<setw(4) << ts->GetEmcIndex()<< " "
- <<setw(4) << ts->GetCpvIndex()<< " "
- <<setw(4) << ts->GetPpsdIndex()<< endl ;
+ <<setw(4) << ts->GetCpvIndex()<< " " << endl ;
}
cout << "-------------------------------------------------------"<< endl ;
virtual ~ AliPHOSTrackSegmentMakerv1() ; // dtor
- virtual char* GetRecPointsBranch (void)const {return (char*)fRecPointsBranchTitle.Data() ;}
- virtual char* GetTrackSegmentsBranch(void)const {return (char*)fTrackSegmentsBranchTitle.Data() ;}
+ virtual char* GetRecPointsBranch (void)const{return (char*)fRecPointsBranchTitle.Data() ;}
+ virtual char* GetTrackSegmentsBranch(void)const{return (char*)fTrackSegmentsBranchTitle.Data() ;}
virtual const Int_t GetTrackSegmentsInRun()const {return fTrackSegmentsInRun ;}
virtual void Exec(Option_t * option) ;
Float_t GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * EmcClu , AliPHOSRecPoint * Ppsd , Bool_t & TooFar )const ; // see R0
void Init() ;
void PrintTrackSegments(Option_t *option) ;
- virtual Bool_t ReadRecPoints(Int_t event) ;
virtual void WriteTrackSegments(Int_t event) ;
private:
Float_t fR0 ; // Maximum distance between a EMC RecPoint and a PPSD RecPoint
- TClonesArray * fLinkLowArray ; //!
TClonesArray * fLinkUpArray ; //!
-
Int_t fEmcFirst; //! Index of first EMC RecPoint belonging to currect PHOS module
Int_t fEmcLast ; //!
Int_t fCpvFirst; //! Cpv upper layer
Int_t fCpvLast; //!
- Int_t fPpsdFirst; //! Cpv low layer
- Int_t fPpsdLast; //!
Int_t fModule ; //! number of module being processed
Int_t fTrackSegmentsInRun ; //! Total number of track segments in one run
fQAHitsMulB = 0 ;
fQATotEner = 0 ;
fQATotEnerB = 0 ;
+
+ fLightYieldMean = 0. ;
+ fIntrinsicPINEfficiency = 0. ;
+ fLightYieldAttenuation = 0. ;
+ fRecalibrationFactor = 0. ;
+ fElectronsPerGeV = 0. ;
+ fAPDGain = 0. ;
+
}
//____________________________________________________________________________
AliPHOSv1::AliPHOSv1(const char *name, const char *title):
AliPHOSv0(name,title)
{
- // ctor : title is used to identify the layout
- // GPS2 = 5 modules (EMC + PPSD)
- // IHEP = 5 modules (EMC + CPV )
- // MIXT = 4 modules (EMC + CPV ) and 1 module (EMC + PPSD)
//
// We store hits :
// - fHits (the "normal" one), which retains the hits associated with
fIshunt = 1 ; // All hits are associated with primary particles
+ //Photoelectron statistics:
+ // The light yield is a poissonian distribution of the number of
+ // photons created in the PbWo4 crystal, calculated using following formula
+ // NumberOfPhotons = EnergyLost * LightYieldMean* APDEfficiency *
+ // exp (-LightYieldAttenuation * DistanceToPINdiodeFromTheHit);
+ // LightYieldMean is parameter calculated to be over 47000 photons per GeV
+ // APDEfficiency is 0.02655
+ // k_0 is 0.0045 from Valery Antonenko
+ // The number of electrons created in the APD is
+ // NumberOfElectrons = APDGain * LightYield
+ // The APD Gain is 300
+ fLightYieldMean = 47000;
+ fIntrinsicPINEfficiency = 0.02655 ; //APD= 0.1875/0.1271 * 0.018 (PIN)
+ fLightYieldAttenuation = 0.0045 ;
+ fRecalibrationFactor = 13.418/ fLightYieldMean ;
+ fElectronsPerGeV = 2.77e+8 ;
+ fAPDGain= 300. ;
+
Int_t nb = GetGeometry()->GetNModules() ;
// create checkables
void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id, Float_t * hits)
{
// 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
+ // A PHOS hit is the sum of all hits in a single crystal from one primary and within soem taime gate
Int_t hitCounter ;
AliPHOSHit *newHit ;
for ( hitCounter = fNhits-1 ; hitCounter >= 0 && !deja ; hitCounter-- ) {
curHit = (AliPHOSHit*) (*fHits)[hitCounter] ;
- if(curHit->GetPrimary() != primary) break ; // We add hits with the same primary, while GEANT treats primaries succesively
+ if(curHit->GetPrimary() != primary) break ;
+ // We add hits with the same primary, while GEANT treats primaries succesively
if( *curHit == *newHit ) {
*curHit + *newHit ;
deja = kTRUE ;
if ( !deja ) {
new((*fHits)[fNhits]) AliPHOSHit(*newHit) ;
// get the block Id number
- Int_t * relid = new Int_t[geom->GetNModules()] ;
+ Int_t relid[4] ;
geom->AbsToRelNumbering(Id, relid) ;
// and fill the relevant QA checkable (only if in PbW04)
if ( relid[1] == 0 ) {
fQAHitsMul->Update(1) ;
((AliPHOSQAIntCheckable*)(*fQAHitsMulB)[relid[0]-1])->Update(1) ;
}
- delete relid ;
fNhits++ ;
}
}
}
}
-
//____________________________________________________________________________
void AliPHOSv1::StepManager(void)
{
-
- // Accumulates hits as long as the track stays in a single crystal or PPSD gas Cell
+ // Accumulates hits as long as the track stays in a single crystal or CPV gas Cell
Int_t relid[4] ; // (box, layer, row, column) indices
Int_t absid ; // absolute cell ID number
- Float_t xyze[4]={0,0,0,0} ; // position wrt MRS and energy deposited
+ Float_t xyze[5]={-1000,-1000,-1000,0,0} ; // position wrt MRS, time and energy deposited
TLorentzVector pos ; // Lorentz vector of the track current position
Int_t copy ;
Int_t primary = gAlice->GetPrimary( gAlice->CurrentTrack() );
TString name = GetGeometry()->GetName() ;
-
- if ( name == "GPS2" || name == "MIXT" ) { // ======> CPV is a GPS' PPSD
-
- if( gMC->CurrentVolID(copy) == gMC->VolId("PPCE") ) // We are inside a gas cell
- {
- gMC->TrackPosition(pos) ;
- xyze[0] = pos[0] ;
- xyze[1] = pos[1] ;
- xyze[2] = pos[2] ;
- xyze[3] = gMC->Edep() ;
-
- 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] += GetGeometry()->GetNModules() - GetGeometry()->GetNPPSDModules();
- }
- gMC->CurrentVolOffID(3, relid[1]) ; // get the Micromegas Module number
- // 1-> GetGeometry()->GetNumberOfModulesPhi() * GetGeometry()->GetNumberOfModulesZ() upper
- // > GetGeometry()->GetNumberOfModulesPhi() * GetGeometry()->GetNumberOfModulesZ() lower
- gMC->CurrentVolOffID(1, relid[2]) ; // get the row number of the cell
- gMC->CurrentVolID(relid[3]) ; // get the column number
-
- // get the absolute Id number
-
- GetGeometry()->RelToAbsNumbering(relid, absid) ;
-
- // add current hit to the hit list
- AddHit(fIshunt, primary, tracknumber, absid, xyze);
-
-
- } // there is deposited energy
- } // We are inside the gas of the CPV
- } // GPS2 configuration
-
- if ( name == "IHEP" || name == "MIXT" ) { // ======> CPV is a IHEP's one
-
- // Yuri Kharlov, 28 September 2000
-
- if( gMC->CurrentVolID(copy) == gMC->VolId("PCPQ") &&
- (gMC->IsTrackEntering() ) &&
- gMC->TrackCharge() != 0) {
-
- gMC -> TrackPosition(pos);
- Float_t xyzm[3], xyzd[3] ;
- Int_t i;
- for (i=0; i<3; i++) xyzm[i] = pos[i];
- gMC -> Gmtod (xyzm, xyzd, 1); // transform coordinate from master to daughter system
-
- Float_t xyd[3]={0,0,0} ; //local posiiton of the entering
- xyd[0] = xyzd[0];
- xyd[1] =-xyzd[1];
- xyd[2] =-xyzd[2];
-
+ Int_t moduleNumber ;
+
+ if( gMC->CurrentVolID(copy) == gMC->VolId("PCPQ") &&
+ (gMC->IsTrackEntering() ) &&
+ gMC->TrackCharge() != 0) {
- // Current momentum of the hit's track in the local ref. system
- TLorentzVector pmom ; //momentum of the particle initiated hit
- 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];
-
- // Digitize the current CPV hit:
-
- // 1. find pad response and
+ gMC -> TrackPosition(pos);
+
+ Float_t xyzm[3], xyzd[3] ;
+ Int_t i;
+ for (i=0; i<3; i++) xyzm[i] = pos[i];
+ gMC -> Gmtod (xyzm, xyzd, 1); // transform coordinate from master to daughter system
+
+ Float_t xyd[3]={0,0,0} ; //local posiiton of the entering
+ xyd[0] = xyzd[0];
+ xyd[1] =-xyzd[2];
+ xyd[2] =-xyzd[1];
+
+ // Current momentum of the hit's track in the local ref. system
+ TLorentzVector pmom ; //momentum of the particle initiated hit
+ 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[2];
+ pmom[2] =-pd[1];
- Int_t moduleNumber;
- gMC->CurrentVolOffID(3,moduleNumber);
- moduleNumber--;
-
-
- TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0); // array of digits for current hit
- CPVDigitize(pmom,xyd,moduleNumber,cpvDigits);
+ // Digitize the current CPV hit:
+
+ // 1. find pad response and
+ gMC->CurrentVolOffID(3,moduleNumber);
+ moduleNumber--;
+
+ TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0); // array of digits for current hit
+ CPVDigitize(pmom,xyd,moduleNumber,cpvDigits);
- Float_t xmean = 0;
- Float_t zmean = 0;
- Float_t qsum = 0;
- Int_t idigit,ndigits;
-
- // 2. go through the current digit list and sum digits in pads
-
- ndigits = cpvDigits->GetEntriesFast();
- for (idigit=0; idigit<ndigits-1; idigit++) {
- AliPHOSCPVDigit *cpvDigit1 = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(idigit);
- Float_t x1 = cpvDigit1->GetXpad() ;
- Float_t z1 = cpvDigit1->GetYpad() ;
- for (Int_t jdigit=idigit+1; jdigit<ndigits; jdigit++) {
- AliPHOSCPVDigit *cpvDigit2 = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(jdigit);
- Float_t x2 = cpvDigit2->GetXpad() ;
- Float_t z2 = cpvDigit2->GetYpad() ;
- if (x1==x2 && z1==z2) {
- Float_t qsum = cpvDigit1->GetQpad() + cpvDigit2->GetQpad() ;
- cpvDigit2->SetQpad(qsum) ;
- cpvDigits->RemoveAt(idigit) ;
- }
+ Float_t xmean = 0;
+ Float_t zmean = 0;
+ Float_t qsum = 0;
+ Int_t idigit,ndigits;
+
+ // 2. go through the current digit list and sum digits in pads
+
+ ndigits = cpvDigits->GetEntriesFast();
+ for (idigit=0; idigit<ndigits-1; idigit++) {
+ AliPHOSCPVDigit *cpvDigit1 = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(idigit);
+ Float_t x1 = cpvDigit1->GetXpad() ;
+ Float_t z1 = cpvDigit1->GetYpad() ;
+ for (Int_t jdigit=idigit+1; jdigit<ndigits; jdigit++) {
+ AliPHOSCPVDigit *cpvDigit2 = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(jdigit);
+ Float_t x2 = cpvDigit2->GetXpad() ;
+ Float_t z2 = cpvDigit2->GetYpad() ;
+ if (x1==x2 && z1==z2) {
+ Float_t qsum = cpvDigit1->GetQpad() + cpvDigit2->GetQpad() ;
+ cpvDigit2->SetQpad(qsum) ;
+ cpvDigits->RemoveAt(idigit) ;
}
}
- cpvDigits->Compress() ;
-
- // 3. add digits to temporary hit list fTmpHits
+ }
+ cpvDigits->Compress() ;
+
+ // 3. add digits to temporary hit list fTmpHits
+
+ ndigits = cpvDigits->GetEntriesFast();
+ for (idigit=0; idigit<ndigits; idigit++) {
+ AliPHOSCPVDigit *cpvDigit = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(idigit);
+ relid[0] = moduleNumber + 1 ; // CPV (or PHOS) module number
+ relid[1] =-1 ; // means CPV
+ relid[2] = cpvDigit->GetXpad() ; // column number of a pad
+ relid[3] = cpvDigit->GetYpad() ; // row number of a pad
+
+ // get the absolute Id number
+ GetGeometry()->RelToAbsNumbering(relid, absid) ;
+
+ // add current digit to the temporary hit list
- ndigits = cpvDigits->GetEntriesFast();
- for (idigit=0; idigit<ndigits; idigit++) {
- AliPHOSCPVDigit *cpvDigit = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(idigit);
- relid[0] = moduleNumber + 1 ; // CPV (or PHOS) module number
- relid[1] =-1 ; // means CPV
- relid[2] = cpvDigit->GetXpad() ; // column number of a pad
- relid[3] = cpvDigit->GetYpad() ; // row number of a pad
-
- // get the absolute Id number
- GetGeometry()->RelToAbsNumbering(relid, absid) ;
-
- // add current digit to the temporary hit list
- xyze[0] = 0. ;
- xyze[1] = 0. ;
- xyze[2] = 0. ;
- xyze[3] = cpvDigit->GetQpad() ; // amplitude in a pad
- primary = -1; // No need in primary for CPV
- AddHit(fIshunt, primary, tracknumber, absid, xyze);
-
- if (cpvDigit->GetQpad() > 0.02) {
- xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5);
- zmean += cpvDigit->GetQpad() * (cpvDigit->GetYpad() + 0.5);
- qsum += cpvDigit->GetQpad();
- }
+ xyze[3] = gMC->TrackTime() ;
+ xyze[4] = cpvDigit->GetQpad() ; // amplitude in a pad
+ primary = -1; // No need in primary for CPV
+ AddHit(fIshunt, primary, tracknumber, absid, xyze);
+
+ if (cpvDigit->GetQpad() > 0.02) {
+ xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5);
+ zmean += cpvDigit->GetQpad() * (cpvDigit->GetYpad() + 0.5);
+ qsum += cpvDigit->GetQpad();
}
- delete cpvDigits;
}
- } // end of IHEP configuration
-
+ delete cpvDigits;
+ }
+
+
if(gMC->CurrentVolID(copy) == gMC->VolId("PXTL") ) { // We are inside a PBWO crystal
+
gMC->TrackPosition(pos) ;
xyze[0] = pos[0] ;
xyze[1] = pos[1] ;
xyze[2] = pos[2] ;
- xyze[3] = gMC->Edep() ;
-
+ Float_t global[3], local[3] ;
+ global[0] = pos[0] ;
+ global[1] = pos[1] ;
+ global[2] = pos[2] ;
+ Float_t lostenergy = gMC->Edep();
- if ( xyze[3] != 0 ) { // Track is inside the crystal and deposits some energy
+ if ( lostenergy != 0 ) { // Track is inside the crystal and deposits some energy
+
+ xyze[3] = gMC->TrackTime() ;
- gMC->CurrentVolOffID(10, relid[0]) ; // get the PHOS module number ;
+ gMC->CurrentVolOffID(10, moduleNumber) ; // get the PHOS module number ;
// fill the relevant QA Checkables
- fQATotEner->Update( xyze[3] ) ; // total energy in PHOS
- ((AliPHOSQAFloatCheckable*)(*fQATotEnerB)[relid[0]-1])->Update( xyze[3] ) ; // energy in this block
-
- if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(10),"PHO1") == 0 )
- relid[0] += GetGeometry()->GetNModules() - GetGeometry()->GetNPPSDModules();
+ fQATotEner->Update( xyze[4] ) ; // total energy in PHOS
+ ((AliPHOSQAFloatCheckable*)(*fQATotEnerB)[moduleNumber-1])->Update( xyze[4] ) ; // energy in this block
- 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
- GetGeometry()->RelToAbsNumbering(relid, absid) ;
+ Int_t strip ;
+ gMC->CurrentVolOffID(3, strip);
+ Int_t cell ;
+ gMC->CurrentVolOffID(2, cell);
- // add current hit to the hit list
- AddHit(fIshunt, primary,tracknumber, absid, xyze);
+ Int_t row = 1 + GetGeometry()->GetNZ() - strip % GetGeometry()->GetNZ() ;
+ Int_t col = (Int_t) TMath::Ceil((Double_t) strip/GetGeometry()->GetNZ()) -1 ;
+ absid = (moduleNumber-1)*GetGeometry()->GetNCristalsInModule() +
+ row + (col*GetGeometry()->GetEMCAGeometry()->GetNCellsInStrip() + cell-1)*GetGeometry()->GetNZ() ;
+ gMC->Gmtod(global, local, 1) ;
+
+ //Calculates the light yield, the number of photns produced in the
+ //crystal
+ Float_t lightYield = gRandom->Poisson(fLightYieldMean * lostenergy *
+ fIntrinsicPINEfficiency *
+ exp(-fLightYieldAttenuation *
+ (local[1]+GetGeometry()->GetCrystalSize(1)/2.0 ))
+ ) ;
+ //Calculates de energy deposited in the crystal
+ xyze[4] = (fRecalibrationFactor/100.) * fAPDGain * lightYield ;
+
+ // add current hit to the hit list
+ AddHit(fIshunt, primary,tracknumber, absid, xyze);
+
+
} // there is deposited energy
} // we are inside a PHOS Xtal
+
}
//____________________________________________________________________________
}
virtual void StepManager(void) ;
- virtual TString Version(void){
- // returns the version number
- return TString("v1") ;
- }
+ virtual TString Version(void){ return TString("v1") ; }
AliPHOSv1 & operator = (const AliPHOSv1 & rvalue) {
// assignement operator requested by coding convention but not needed
Float_t CPVPadResponseFunction(Float_t qhit, Float_t zg, Float_t xg) ;
Double_t CPVCumulPadResponse(Double_t x, Double_t y) ;
+ //Variables conserning light yeild and APD efficiency
+ Float_t GetLightYieldMean() const { return fLightYieldMean ;}
+ Float_t GetLightYieldAttenuation() const { return fLightYieldAttenuation ;}
+ Float_t GetRecalibrationFactor() const { return fRecalibrationFactor ;}
+ Float_t GetAPDGain() const { return fAPDGain ;}
+ Float_t GetIntrinsicPINEfficiency() const { return fIntrinsicPINEfficiency ;}
+ Float_t GetElectronsPerGeV() const { return fElectronsPerGeV ;}
+
+ void SetLightYieldMean(Float_t LightYieldMean)
+ {fLightYieldMean = LightYieldMean;}
+ void SetLightYieldAttenuation(Float_t LightYieldAttenuation)
+ {fLightYieldAttenuation = LightYieldAttenuation;}
+ void SetIntrinsicPINEfficiency(Float_t IntrinsicPINEfficiency)
+ {fIntrinsicPINEfficiency = IntrinsicPINEfficiency;}
+ void SetRecalibrationFactor(Float_t RecalibrationFactor)
+ {fRecalibrationFactor = RecalibrationFactor;}
+ void SetElectronsPerGeV(Float_t ElectronsPerGeV)
+ {fElectronsPerGeV = ElectronsPerGeV;}
+ void SetAPDGain(Float_t APDGain) {fAPDGain = APDGain;}
+
protected:
-
AliPHOSQAIntCheckable * fQAHitsMul ; // QA Hits Multiplicity checkable
- TClonesArray * fQAHitsMulB ; // QA Hits Multiplicity / Block checkable
+ TClonesArray * fQAHitsMulB ; // QA Hits Multiplicity / Block checkable
AliPHOSQAFloatCheckable * fQATotEner ; // QA Total Energy checkable
- TClonesArray * fQATotEnerB ; // QA Total Energy / Block checkable
+ TClonesArray * fQATotEnerB ; // QA Total Energy / Block checkable
+
+ Float_t fLightYieldMean ; // Mean lightyield in the PbOW4 xtal per GeV (Poisson distribution)
+ Float_t fIntrinsicPINEfficiency ; // Photo efficiency of the PIN diode
+ Float_t fLightYieldAttenuation ; // Attenuation of the light through the crystal
+ Float_t fRecalibrationFactor ; // Recalibration factor
+ Float_t fElectronsPerGeV ; // Number of electrons per GeV created in the PIN by a ionizing particle
+ Float_t fAPDGain ; // APD Gain
ClassDef(AliPHOSv1,1) // Implementation of PHOS manager class for layout EMC+PPSD
ClassImp(AliPHOSvImpacts)
//____________________________________________________________________________
-AliPHOSvImpacts::AliPHOSvImpacts()
+AliPHOSvImpacts::AliPHOSvImpacts():AliPHOSv1()
{
// ctor
}
AliPHOSv1(name,title)
{
// ctor : title is used to identify the layout
- // GPS2 = 5 modules (EMC + PPSD)
- // IHEP = 5 modules (EMC + CPV )
- // MIXT = 4 modules (EMC + CPV ) and 1 module (EMC + PPSD)
//
// We store hits :
// - fHits (the "normal" one), which retains the hits associated with
// This part inherits from AliPHOSv1
//
// We store impacts :
- // - fEMCImpacts, fCPVImpacts, fPPSDImpacts which are
- // TList of EMC, CPV and PPSD modules respectively, each
+ // - fEMCImpacts, fCPVImpacts which are
+ // TList of EMC and CPV modules respectively, each
// modules contains TClonesArray of AliPHOSImpacts
fEMCImpacts = new TList();
fCPVImpacts = new TList();
- fPPSDImpacts = new TList();
Int_t nPHOSModules = GetGeometry()->GetNModules();
- Int_t nCPVModules = GetGeometry()->GetNCPVModules();
- Int_t nPPSDModules = GetGeometry()->GetNPPSDModules();
+ Int_t nCPVModules = GetGeometry()->GetNModules();
Int_t iPHOSModule;
TClonesArray * impacts;
fNCPVImpacts[iPHOSModule] = 0;
impacts = (TClonesArray *)fCPVImpacts->At(iPHOSModule);
}
- for (iPHOSModule=0; iPHOSModule<nPPSDModules; iPHOSModule++) {
- fPPSDImpacts->Add(new TClonesArray("AliPHOSImpact",200)) ;
- fNPPSDImpacts[iPHOSModule] = 0;
- impacts = (TClonesArray *)fPPSDImpacts->At(iPHOSModule);
- }
}
fHits = 0 ;
}
- // Delete impacts in EMC, CPV and PPSD
+ // Delete impacts in EMC, CPV
if ( fEMCImpacts ) {
fEMCImpacts->Delete() ;
delete fEMCImpacts ;
delete fCPVImpacts ;
fCPVImpacts = 0 ;
}
- if ( fPPSDImpacts ) {
- fPPSDImpacts->Delete() ;
- delete fPPSDImpacts ;
- fPPSDImpacts = 0 ;
- }
}
//____________________________________________________________________________
nImpacts= fNCPVImpacts[module];
fNCPVImpacts[module]++ ;
}
- else if (strcmp(det,"PPSD")==0) {
- impacts = (TClonesArray *)fPPSDImpacts->At(module);
- nImpacts= fNPPSDImpacts[module];
- fNPPSDImpacts[module]++ ;
- }
new((*impacts)[nImpacts]) AliPHOSImpact(shunt,primary,track,pid,p,xyz) ;
Int_t splitlevel = 0 ;
gAlice->TreeH()->Branch("PHOSEmcImpacts" , "TList", &fEMCImpacts , bufferSize, splitlevel);
gAlice->TreeH()->Branch("PHOSCpvImpacts" , "TList", &fCPVImpacts , bufferSize, splitlevel);
- gAlice->TreeH()->Branch("PHOSPpsdImpacts", "TList", &fPPSDImpacts, bufferSize, splitlevel);
}
fNEMCImpacts[i] = 0 ;
}
- if ( strcmp(GetGeometry()->GetName(),"IHEP") == 0 || strcmp(GetGeometry()->GetName(),"MIXT") == 0 ) {
- for (i=0; i<GetGeometry()->GetNCPVModules(); i++) {
- ((TClonesArray*)fCPVImpacts->At(i)) -> Clear();
- fNCPVImpacts[i] = 0 ;
- }
- }
-
- if ( strcmp(GetGeometry()->GetName(),"GPS2") == 0 || strcmp(GetGeometry()->GetName(),"MIXT") == 0 ) {
- for (i=0; i<GetGeometry()->GetNPPSDModules(); i++) {
- ((TClonesArray*)fPPSDImpacts->At(i)) -> Clear();
- fNPPSDImpacts[i] = 0 ;
- }
+ for (i=0; i<GetGeometry()->GetNModules(); i++) {
+ ((TClonesArray*)fCPVImpacts->At(i)) -> Clear();
+ fNCPVImpacts[i] = 0 ;
}
}
//_____________________________________________________________________________
void AliPHOSvImpacts::StepManager(void)
{
- // Find impacts (tracks which enter the EMC, CPV or PPSD)
+ // Find impacts (tracks which enter the EMC, CPV)
// and add them to to the impact lists
AliPHOSv1::StepManager();
Int_t pid = gMC->TrackPid();
Int_t module;
gMC->CurrentVolOffID(10,module);
- if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(10),"PHO1") == 0 )
- module += GetGeometry()->GetNModules() - GetGeometry()->GetNPPSDModules();
module--;
AddImpact("EMC ",fIshunt, primary,tracknumber, module, pid, pmom, xyzm);
}
// Add impact to CPV
- if( (name == "IHEP" || name == "MIXT") &&
- gMC->CurrentVolID(copy) == gMC->VolId("PCPQ") &&
+ if( gMC->CurrentVolID(copy) == gMC->VolId("PCPQ") &&
gMC->IsTrackEntering() ) {
gMC->TrackMomentum(pmom);
gMC->TrackPosition(pos) ;
module--;
AddImpact("CPV ",fIshunt, primary,tracknumber, module, pid, pmom, xyzm);
}
-
- // Add impact to PPSD
-
- if( (name == "GPS2" || name == "MIXT") &&
- gMC->CurrentVolID(copy) == gMC->VolId("PPCE") &&
- gMC->IsTrackEntering() ) {
- gMC->TrackMomentum(pmom);
- gMC->TrackPosition(pos) ;
-
- Int_t i;
- for (i=0; i<3; i++) xyzm[i] = pos[i];
-
- for (i=0; i<3; i++) {
- xyzm[i] = pos[i] ;
- pm[i] = pmom[i];
- }
- gMC -> Gmtod (xyzm, xyzd, 1); // transform coordinate from master to daughter system
- gMC -> Gmtod (pm, pd, 2); // transform 3-momentum from master to daughter system
-
- // Select tracks coming to the crystal from up or down sides
- if (pd[1]<0 && xyzd[1] > (GetGeometry()->GetConversionGap() + GetGeometry()->GetAvalancheGap())/2-0.001 ||
- pd[1]>0 && xyzd[1] < -(GetGeometry()->GetConversionGap() + GetGeometry()->GetAvalancheGap())/2+0.001) {
- Int_t pid = gMC->TrackPid();
- Int_t module;
- gMC->CurrentVolOffID(5,module);
- module--;
- AddImpact("PPSD",fIshunt, primary,tracknumber, module, pid, pmom, xyzm);
- }
- }
+
}
# C++ sources
-SRCS = AliPHOS.cxx AliPHOSv0.cxx AliPHOSv1.cxx AliPHOSv2.cxx \
- AliPHOSv3.cxx AliPHOSv4.cxx AliPHOSvImpacts.cxx \
+SRCS = AliPHOS.cxx AliPHOSv0.cxx AliPHOSv1.cxx \
+ AliPHOSvFast.cxx AliPHOSvImpacts.cxx \
AliPHOSImpact.cxx \
AliPHOSHit.cxx \
AliPHOSGeometry.cxx \
AliPHOSEMCAGeometry.cxx \
AliPHOSCPVGeometry.cxx AliPHOSCPVBaseGeometry.cxx \
- AliPHOSPPSDGeometry.cxx AliPHOSSupportGeometry.cxx \
+ AliPHOSSupportGeometry.cxx \
AliPHOSCPVDigit.cxx AliPHOSDigit.cxx \
AliPHOSRecPoint.cxx AliPHOSEmcRecPoint.cxx \
- AliPHOSPpsdRecPoint.cxx AliPHOSCpvRecPoint.cxx \
+ AliPHOSCpvRecPoint.cxx \
AliPHOSClusterizer.cxx AliPHOSClusterizerv1.cxx \
AliPHOSLink.cxx AliPHOSSDigitizer.cxx AliPHOSDigitizer.cxx\
AliPHOSReconstructioner.cxx AliPHOSTrackSegment.cxx \
-include tgt_$(ALICE_TARGET)/Make-depend
-test:
- @echo " ____________________________________________________________ "
- @echo " "
- @echo " Starting the test of the simulation/reconstruction software. Please don't take the warning messages into account. "
- @echo " ____________________________________________________________ "
- @aliroot -b -q testsim.C > out
- @aliroot -b -q testsimexam.C > out
- @aliroot -b -q testreconSDigits.C > out
- @aliroot -b -q testreconDigits.C > out
- @aliroot -b -q testreconRecPoints.C > out
- @aliroot -b -q testreconTrackSegments.C > out
- @aliroot -b -q testreconRecParticles.C > out
- @rm out
-# @rm galice.root
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff