From 9688c1ddeb34093dbf8265987926402e715ba192 Mon Sep 17 00:00:00 2001 From: schutz Date: Wed, 26 Sep 2001 12:25:54 +0000 Subject: [PATCH] Timing added to reconstruction --- PHOS/AliPHOSAnalyze.cxx | 784 ++++++++++++++-------------- PHOS/AliPHOSClusterizer.h | 4 +- PHOS/AliPHOSClusterizerv1.cxx | 257 +++------ PHOS/AliPHOSClusterizerv1.h | 23 +- PHOS/AliPHOSCpvRecPoint.cxx | 77 ++- PHOS/AliPHOSCpvRecPoint.h | 7 +- PHOS/AliPHOSDigit.cxx | 10 +- PHOS/AliPHOSDigit.h | 7 +- PHOS/AliPHOSDigitizer.cxx | 240 ++++++--- PHOS/AliPHOSDigitizer.h | 28 +- PHOS/AliPHOSEmcRecPoint.cxx | 195 +++---- PHOS/AliPHOSEmcRecPoint.h | 13 +- PHOS/AliPHOSFastRecParticle.cxx | 20 +- PHOS/AliPHOSFastRecParticle.h | 5 +- PHOS/AliPHOSHit.cxx | 52 +- PHOS/AliPHOSHit.h | 8 +- PHOS/AliPHOSPID.h | 2 +- PHOS/AliPHOSPIDv1.cxx | 246 +++------ PHOS/AliPHOSPIDv1.h | 12 +- PHOS/AliPHOSQAChecker.cxx | 6 +- PHOS/AliPHOSQAChecker.h | 7 +- PHOS/AliPHOSQAMeanChecker.cxx | 6 +- PHOS/AliPHOSQAVirtualCheckable.h | 5 +- PHOS/AliPHOSSDigitizer.cxx | 16 +- PHOS/AliPHOSTrackSegment.cxx | 12 +- PHOS/AliPHOSTrackSegment.h | 7 +- PHOS/AliPHOSTrackSegmentMaker.h | 1 - PHOS/AliPHOSTrackSegmentMakerv1.cxx | 265 ++-------- PHOS/AliPHOSTrackSegmentMakerv1.h | 9 +- PHOS/AliPHOSv1.cxx | 309 +++++------ PHOS/AliPHOSv1.h | 37 +- PHOS/AliPHOSvImpacts.cxx | 83 +-- PHOS/Makefile | 23 +- 33 files changed, 1212 insertions(+), 1564 deletions(-) diff --git a/PHOS/AliPHOSAnalyze.cxx b/PHOS/AliPHOSAnalyze.cxx index 9c74233fd22..a34e9c3f3c8 100644 --- a/PHOS/AliPHOSAnalyze.cxx +++ b/PHOS/AliPHOSAnalyze.cxx @@ -89,7 +89,6 @@ #include "AliPHOSTrackSegment.h" #include "AliPHOSRecParticle.h" #include "AliPHOSCpvRecPoint.h" -#include "AliPHOSPpsdRecPoint.h" #include "AliPHOSGetter.h" @@ -128,24 +127,32 @@ void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod,const char * branchName,c //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 ; @@ -154,7 +161,8 @@ void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod,const char * branchName,c { 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) ; @@ -178,114 +186,98 @@ void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod,const char * branchName,c } } } + 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 ; @@ -316,20 +308,19 @@ void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod,const char * branchName,c } } } - + + } //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") ; @@ -422,29 +413,33 @@ void AliPHOSAnalyze::Ls(){ 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 ; @@ -474,14 +469,14 @@ void AliPHOSAnalyze::Ls(){ 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 @@ -543,7 +538,6 @@ void AliPHOSAnalyze::Ls(){ 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 @@ -551,67 +545,82 @@ void AliPHOSAnalyze::Ls(){ 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 " <TrackSegments() ; + if(!ts) { + cout << "AliPHOSAnalyze::EnergyResolution --> Event " <EmcRecPoints() ; + if(!emcrp){ + cout << "AliPHOSAnalyze::EnergyResolution --> Event " <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) ; @@ -620,7 +629,7 @@ void AliPHOSAnalyze::Ls(){ // efile->Write() ; efile->Close() ; delete efile ; - + } //____________________________________________________________________________ void AliPHOSAnalyze::PositionResolution(const char * branchTitle) @@ -666,77 +675,91 @@ 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 " <TrackSegments() ; + if(!ts) { + cout << "AliPHOSAnalyze::PositionResolution --> Event " <EmcRecPoints() ; + if(!emcrp){ + cout << "AliPHOSAnalyze::PositionResolution --> Event " <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 @@ -749,8 +772,8 @@ void AliPHOSAnalyze::PositionResolution(const char * branchTitle) pfile->Write() ; pfile->Close() ; delete pfile ; - - + + } //____________________________________________________________________________ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){ @@ -762,7 +785,6 @@ 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 @@ -770,27 +792,22 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){ //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 @@ -806,9 +823,6 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){ 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.); @@ -824,9 +838,6 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){ 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"); @@ -835,9 +846,6 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){ 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"); @@ -846,9 +854,6 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){ 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"); @@ -857,10 +862,7 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){ 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) @@ -868,9 +870,6 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){ 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 @@ -888,11 +887,30 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){ 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 " <TrackSegments() ; + if(!ts) { + cout << "AliPHOSAnalyze::Contamination --> Event " <EmcRecPoints() ; + if(!emcrp){ + cout << "AliPHOSAnalyze::Contamination --> Event " <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 } @@ -1074,24 +1066,18 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){ 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(); diff --git a/PHOS/AliPHOSClusterizer.h b/PHOS/AliPHOSClusterizer.h index fd28665031b..7ecfaf91f98 100644 --- a/PHOS/AliPHOSClusterizer.h +++ b/PHOS/AliPHOSClusterizer.h @@ -29,10 +29,10 @@ public: 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 ; @@ -43,10 +43,10 @@ public: 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 ; diff --git a/PHOS/AliPHOSClusterizerv1.cxx b/PHOS/AliPHOSClusterizerv1.cxx index ea7367a2e59..aeaa9cb9c13 100644 --- a/PHOS/AliPHOSClusterizerv1.cxx +++ b/PHOS/AliPHOSClusterizerv1.cxx @@ -78,7 +78,6 @@ #include "AliPHOSDigitizer.h" #include "AliPHOSEmcRecPoint.h" #include "AliPHOS.h" -#include "AliPHOSPpsdRecPoint.h" #include "AliPHOSGetter.h" #include "AliRun.h" @@ -94,7 +93,6 @@ ClassImp(AliPHOSClusterizerv1) fCpvClusteringThreshold = 0.0; fEmcClusteringThreshold = 0.2; - fPpsdClusteringThreshold = 0.0000002 ; fEmcLocMaxCut = 0.03 ; fCpvLocMaxCut = 0.03 ; @@ -102,9 +100,13 @@ ClassImp(AliPHOSClusterizerv1) fW0 = 4.5 ; fW0CPV = 4.0 ; + fEmcTimeGate = 1.e-8 ; + + fToUnfold = kTRUE ; + fHeaderFileName = "" ; fDigitsBranchTitle = "" ; - fRecPointsInRun = 0 ; + fRecPointsInRun = 0 ; } //____________________________________________________________________________ @@ -119,29 +121,33 @@ AliPHOSClusterizerv1::AliPHOSClusterizerv1(const char* headerFile,const char* na 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) { @@ -156,30 +162,31 @@ 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 ; } @@ -190,26 +197,24 @@ void AliPHOSClusterizerv1::Exec(Option_t * option) 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")){ @@ -324,10 +329,10 @@ void AliPHOSClusterizerv1::Init() 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 ; @@ -336,13 +341,12 @@ void AliPHOSClusterizerv1::Init() 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) ; } @@ -366,11 +370,12 @@ Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)c 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 { @@ -381,16 +386,11 @@ Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)c } 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 ; } @@ -411,22 +411,6 @@ Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const 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 { @@ -438,66 +422,10 @@ 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) @@ -505,9 +433,10 @@ 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) ; @@ -535,10 +464,6 @@ void AliPHOSClusterizerv1::WriteRecPoints(Int_t event) ((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; @@ -635,7 +560,6 @@ void AliPHOSClusterizerv1::MakeClusters() 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 ; @@ -655,16 +579,13 @@ void AliPHOSClusterizerv1::MakeClusters() } 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() ; @@ -721,14 +642,14 @@ void AliPHOSClusterizerv1::MakeClusters() } // 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() ; @@ -774,7 +695,7 @@ void AliPHOSClusterizerv1::MakeUnfolding() // Unfold now CPV clusters if(fNumberOfCpvClusters > 0){ - Int_t nModulesToUnfold = geom->GetNCPVModules() ; + Int_t nModulesToUnfold = geom->GetNModules() ; Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ; Int_t index ; @@ -1052,8 +973,7 @@ void AliPHOSClusterizerv1::Print(Option_t * option)const << " 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 @@ -1081,57 +1001,40 @@ void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option) 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; iprimaryGetEntries() ; 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; iprimaryClassName() , "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 ; + } //______________________________________________________________________________ diff --git a/PHOS/AliPHOSCpvRecPoint.h b/PHOS/AliPHOSCpvRecPoint.h index 8902a62e182..0e14b3bada9 100644 --- a/PHOS/AliPHOSCpvRecPoint.h +++ b/PHOS/AliPHOSCpvRecPoint.h @@ -45,11 +45,10 @@ public: 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 diff --git a/PHOS/AliPHOSDigit.cxx b/PHOS/AliPHOSDigit.cxx index 2bf6f24a688..9ffb54cfe2f 100644 --- a/PHOS/AliPHOSDigit.cxx +++ b/PHOS/AliPHOSDigit.cxx @@ -21,7 +21,7 @@ // 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 --- @@ -48,12 +48,13 @@ ClassImp(AliPHOSDigit) } //____________________________________________________________________________ -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){ @@ -80,6 +81,7 @@ AliPHOSDigit::AliPHOSDigit(const AliPHOSDigit & digit) 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 ; @@ -159,6 +161,8 @@ AliPHOSDigit& AliPHOSDigit::operator+(AliPHOSDigit const & digit) // if amplitude is larger than fAmp += digit.fAmp ; + if(fTime > digit.fTime) + fTime = digit.fTime ; Int_t max1 = fNprimary ; @@ -188,7 +192,7 @@ ostream& operator << ( ostream& out , const AliPHOSDigit & digit) { // 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;iGetNumberOfCPVPadsZ()*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) @@ -200,49 +184,130 @@ void AliPHOSDigitizer::Digitize(const Int_t event) 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; iAt(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; iAt(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; iAt(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; iAt(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; iAt(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) @@ -251,18 +316,13 @@ void AliPHOSDigitizer::Digitize(const Int_t event) 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) ; @@ -369,6 +429,19 @@ void AliPHOSDigitizer::Exec(Option_t *option) } +//____________________________________________________________________________ +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() { @@ -478,14 +551,13 @@ void AliPHOSDigitizer::Print(Option_t* option)const { 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 @@ -496,28 +568,51 @@ void AliPHOSDigitizer::PrintDigits(Option_t * option){ 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; iprimaryGetNprimary(); 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; iprimaryGetNprimary(); iprimary++) + cout << setw(5) << digit->GetPrimary(iprimary+1) << " "; + cout << endl; + } + } } + } //__________________________________________________________________ @@ -529,6 +624,13 @@ void AliPHOSDigitizer::SetSDigitsBranch(const char* title) 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() diff --git a/PHOS/AliPHOSDigitizer.h b/PHOS/AliPHOSDigitizer.h index 9dafc4452af..6cbdacb729b 100644 --- a/PHOS/AliPHOSDigitizer.h +++ b/PHOS/AliPHOSDigitizer.h @@ -44,10 +44,8 @@ public: 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 ; @@ -57,8 +55,6 @@ public: 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) ; @@ -72,21 +68,33 @@ private: 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 diff --git a/PHOS/AliPHOSEmcRecPoint.cxx b/PHOS/AliPHOSEmcRecPoint.cxx index 2ab02e4d710..0cd2c50275c 100644 --- a/PHOS/AliPHOSEmcRecPoint.cxx +++ b/PHOS/AliPHOSEmcRecPoint.cxx @@ -51,6 +51,7 @@ AliPHOSEmcRecPoint::AliPHOSEmcRecPoint() : AliPHOSRecPoint() fAmp = 0. ; fCoreEnergy = 0 ; fEnergyList = 0 ; + fTime = -1. ; fLocPos.SetX(1000000.) ; //Local position should be evaluated } @@ -174,117 +175,114 @@ Int_t AliPHOSEmcRecPoint::Compare(const TObject * obj) const 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; iDigitDigit(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; iDigitAt(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; iDigitDigit(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; iDigitAt(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; -// } + } } //____________________________________________________________________________ @@ -631,9 +629,20 @@ Int_t AliPHOSEmcRecPoint::GetNumberOfLocalMax(Int_t * maxAt, Float_t * maxAtEn } 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) { diff --git a/PHOS/AliPHOSEmcRecPoint.h b/PHOS/AliPHOSEmcRecPoint.h index 0e19624ca97..32b9dc8e3bf 100644 --- a/PHOS/AliPHOSEmcRecPoint.h +++ b/PHOS/AliPHOSEmcRecPoint.h @@ -40,10 +40,7 @@ public: 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 ;} @@ -59,6 +56,7 @@ public: 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") ; @@ -70,13 +68,18 @@ public: } 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) diff --git a/PHOS/AliPHOSFastRecParticle.cxx b/PHOS/AliPHOSFastRecParticle.cxx index 177d018993b..501136f07a7 100644 --- a/PHOS/AliPHOSFastRecParticle.cxx +++ b/PHOS/AliPHOSFastRecParticle.cxx @@ -164,23 +164,23 @@ TString AliPHOSFastRecParticle::Name() 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 ; } diff --git a/PHOS/AliPHOSFastRecParticle.h b/PHOS/AliPHOSFastRecParticle.h index d754271cdc9..8be17966ee9 100644 --- a/PHOS/AliPHOSFastRecParticle.h +++ b/PHOS/AliPHOSFastRecParticle.h @@ -60,8 +60,9 @@ class AliPHOSFastRecParticle : public TParticle { 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 ; diff --git a/PHOS/AliPHOSHit.cxx b/PHOS/AliPHOSHit.cxx index c43c7a31c22..246d2c1be32 100644 --- a/PHOS/AliPHOSHit.cxx +++ b/PHOS/AliPHOSHit.cxx @@ -19,7 +19,7 @@ // 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 --- @@ -37,15 +37,12 @@ #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 ; @@ -53,8 +50,8 @@ AliPHOSHit::AliPHOSHit(const AliPHOSHit & hit) fELOS = hit.fELOS ; fPrimary = hit.fPrimary ; fTrack = hit.fTrack ; - - + fTime = hit.fTime ; + } //____________________________________________________________________________ @@ -63,41 +60,57 @@ AliPHOSHit::AliPHOSHit(Int_t shunt, Int_t primary, Int_t track, Int_t id, Float_ // // 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 (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 (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 (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 { @@ -117,6 +130,9 @@ AliPHOSHit AliPHOSHit::operator+(const AliPHOSHit &rValue) // Add the energy of the hit fELOS += rValue.GetEnergy() ; + + if(rValue.GetTime() < fTime) + fTime = rValue.GetTime() ; return *this; @@ -127,7 +143,7 @@ ostream& operator << (ostream& out, const AliPHOSHit& hit) { // Print out Id and energy - out << "AliPHOSHit = " << hit.GetId() << " " << hit.GetEnergy() << endl ; + out << "AliPHOSHit = " << hit.GetId() << " " << hit.GetEnergy() << " " << hit.GetTime() << endl ; return out ; } diff --git a/PHOS/AliPHOSHit.h b/PHOS/AliPHOSHit.h index c23dd225e15..e278613cfe6 100644 --- a/PHOS/AliPHOSHit.h +++ b/PHOS/AliPHOSHit.h @@ -49,11 +49,16 @@ class AliPHOSHit : public AliHit { // 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) ; @@ -63,6 +68,7 @@ class AliPHOSHit : public AliHit { 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 diff --git a/PHOS/AliPHOSPID.h b/PHOS/AliPHOSPID.h index 84a8711c7e6..f7ec42b4619 100644 --- a/PHOS/AliPHOSPID.h +++ b/PHOS/AliPHOSPID.h @@ -40,11 +40,11 @@ public: 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 ; diff --git a/PHOS/AliPHOSPIDv1.cxx b/PHOS/AliPHOSPIDv1.cxx index cb383bdefbe..7fa241589a0 100644 --- a/PHOS/AliPHOSPIDv1.cxx +++ b/PHOS/AliPHOSPIDv1.cxx @@ -54,7 +54,7 @@ // //*-- 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" @@ -93,11 +93,12 @@ AliPHOSPIDv1::AliPHOSPIDv1():AliPHOSPID() fFormula = 0 ; fDispersion = 0. ; fCpvEmcDistance = 0 ; + fTimeGate = 2.e-9 ; fHeaderFileName = "" ; fTrackSegmentsTitle= "" ; fRecPointsTitle = "" ; fRecParticlesTitle = "" ; - fIDOptions = "" ; + fIDOptions = "dis time" ; fRecParticlesInRun = 0 ; } @@ -109,19 +110,20 @@ AliPHOSPIDv1::AliPHOSPIDv1(const char * headerFile,const char * name) : AliPHOSP fFormula = new TFormula("LambdaCuts","(x>1)*(x<2.5)*(y>0)*(yDelete() ; + 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 @@ -278,7 +170,7 @@ Float_t AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cp void AliPHOSPIDv1::Exec(Option_t * option) { //Steering method - + if( strcmp(GetName(), "")== 0 ) Init() ; @@ -289,7 +181,7 @@ void AliPHOSPIDv1::Exec(Option_t * option) Print("") ; return ; } - + gAlice->GetEvent(0) ; //check, if the branch with name of this" already exits? TObjArray * lob = (TObjArray*)gAlice->TreeR()->GetListOfBranches() ; @@ -299,7 +191,7 @@ void AliPHOSPIDv1::Exec(Option_t * option) 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 ; @@ -307,7 +199,7 @@ void AliPHOSPIDv1::Exec(Option_t * option) 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 ; @@ -319,20 +211,16 @@ void AliPHOSPIDv1::Exec(Option_t * option) 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 ; @@ -340,7 +228,7 @@ void AliPHOSPIDv1::Exec(Option_t * option) << gBenchmark->GetCpuTime("PHOSPID")/nevents << " seconds per event " << endl ; cout << endl ; } - + } //____________________________________________________________________________ void AliPHOSPIDv1::Init() @@ -360,10 +248,6 @@ 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() ) ; @@ -391,12 +275,14 @@ void AliPHOSPIDv1::MakeRecParticles(){ 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) @@ -406,13 +292,9 @@ void AliPHOSPIDv1::MakeRecParticles(){ 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) ; @@ -432,22 +314,19 @@ void AliPHOSPIDv1::MakeRecParticles(){ 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++ ; } @@ -464,13 +343,15 @@ void AliPHOSPIDv1:: Print(Option_t * option) const 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 ; } @@ -492,7 +373,7 @@ void AliPHOSPIDv1::WriteRecParticles(Int_t event) 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 @@ -515,12 +396,12 @@ void AliPHOSPIDv1::WriteRecParticles(Int_t event) } 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()); @@ -533,7 +414,7 @@ void AliPHOSPIDv1::WriteRecParticles(Int_t event) rpBranch->Fill() ; pidBranch->Fill() ; - + gAlice->TreeR()->Write(0,kOverwrite) ; } @@ -572,14 +453,12 @@ void AliPHOSPIDv1::PlotDispersionCuts()const } //____________________________________________________________________________ -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) ; @@ -628,13 +507,12 @@ void AliPHOSPIDv1::PrintRecParticles(Option_t * option) TString taskName(GetName()) ; taskName.Remove(taskName.Index(Version())-1) ; TClonesArray * recParticles = gime->RecParticles(taskName) ; - - + cout << "AliPHOSPIDv1: event "<GetEvNumber() << endl ; cout << " found " << recParticles->GetEntriesFast() << " RecParticles " << endl ; fRecParticlesInRun += recParticles->GetEntriesFast() ; - + if(strstr(option,"all")) { // printing found TS cout << " PARTICLE " @@ -651,29 +529,29 @@ void AliPHOSPIDv1::PrintRecParticles(Option_t * option) 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 ; } @@ -681,8 +559,8 @@ void AliPHOSPIDv1::PrintRecParticles(Option_t * option) // 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; iprimaryFindObjectAny("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() ; } diff --git a/PHOS/AliPHOSQAChecker.h b/PHOS/AliPHOSQAChecker.h index 619d05b9e8b..c6d9671fc4e 100644 --- a/PHOS/AliPHOSQAChecker.h +++ b/PHOS/AliPHOSQAChecker.h @@ -28,10 +28,7 @@ class AliPHOSQAChecker : public TTask { 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 @@ -51,7 +48,7 @@ public: void StatusAll() { ExecuteTask("S") ; } friend void AliPHOSQAVirtualCheckable::AddChecker(AliPHOSQAChecker * ch) ; - // friend AliPHOSQAVirtualCheckable::AliPHOSQAVirtualCheckable(const char * name) ; + friend AliPHOSQAVirtualCheckable::AliPHOSQAVirtualCheckable(const char * name) ; private: diff --git a/PHOS/AliPHOSQAMeanChecker.cxx b/PHOS/AliPHOSQAMeanChecker.cxx index ad2abf4c385..c3236813ecb 100644 --- a/PHOS/AliPHOSQAMeanChecker.cxx +++ b/PHOS/AliPHOSQAMeanChecker.cxx @@ -89,9 +89,5 @@ TString AliPHOSQAMeanChecker::CheckingOperation() { // 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 ; } diff --git a/PHOS/AliPHOSQAVirtualCheckable.h b/PHOS/AliPHOSQAVirtualCheckable.h index 7a0513e6ece..ac77fdce371 100644 --- a/PHOS/AliPHOSQAVirtualCheckable.h +++ b/PHOS/AliPHOSQAVirtualCheckable.h @@ -30,10 +30,7 @@ class AliPHOSQAVirtualCheckable : public TNamed { 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 diff --git a/PHOS/AliPHOSSDigitizer.cxx b/PHOS/AliPHOSSDigitizer.cxx index 98dff99341e..1a2bbe1b9c9 100644 --- a/PHOS/AliPHOSSDigitizer.cxx +++ b/PHOS/AliPHOSSDigitizer.cxx @@ -195,22 +195,22 @@ void AliPHOSSDigitizer::Exec(Option_t *option) 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++) { diff --git a/PHOS/AliPHOSTrackSegment.cxx b/PHOS/AliPHOSTrackSegment.cxx index 979656b15a7..deeba20449c 100644 --- a/PHOS/AliPHOSTrackSegment.cxx +++ b/PHOS/AliPHOSTrackSegment.cxx @@ -36,8 +36,7 @@ ClassImp(AliPHOSTrackSegment) //____________________________________________________________________________ AliPHOSTrackSegment::AliPHOSTrackSegment( AliPHOSEmcRecPoint * emc , - AliPHOSRecPoint * ppsdrp1, - AliPHOSRecPoint * ppsdrp2 ) + AliPHOSRecPoint * ppsdrp1) { // ctor @@ -51,10 +50,6 @@ AliPHOSTrackSegment::AliPHOSTrackSegment( AliPHOSEmcRecPoint * emc , else fPpsdUpRecPoint = -1 ; - if( ppsdrp2 ) - fPpsdLowRecPoint = ppsdrp2->GetIndexInList() ; - else - fPpsdLowRecPoint = -1 ; fIndexInList = -1 ; } @@ -75,7 +70,6 @@ void AliPHOSTrackSegment::Copy(TObject & obj) TObject::Copy(obj) ; ( (AliPHOSTrackSegment &)obj ).fEmcRecPoint = fEmcRecPoint ; - ( (AliPHOSTrackSegment &)obj ).fPpsdLowRecPoint = fPpsdLowRecPoint ; ( (AliPHOSTrackSegment &)obj ).fPpsdUpRecPoint = fPpsdUpRecPoint ; ( (AliPHOSTrackSegment &)obj ).fIndexInList = fIndexInList ; } @@ -94,10 +88,6 @@ void AliPHOSTrackSegment::Print(Option_t * opt) const else cout << "No CPV RecPoint " << endl ; - if(fPpsdLowRecPoint >= 0) - cout << "PPSD RecPoint # " << fPpsdLowRecPoint << endl ; - else - cout << "No PPSD RecPoint " << endl ; cout << "------------------------------------ " << endl ; diff --git a/PHOS/AliPHOSTrackSegment.h b/PHOS/AliPHOSTrackSegment.h index 10f25b53ebb..2f6e55c5666 100644 --- a/PHOS/AliPHOSTrackSegment.h +++ b/PHOS/AliPHOSTrackSegment.h @@ -30,8 +30,7 @@ public: AliPHOSTrackSegment() {} AliPHOSTrackSegment(AliPHOSEmcRecPoint * EmcRecPoint , - AliPHOSRecPoint * PpsdUp, - AliPHOSRecPoint * PpsdLow ) ; // ctor + AliPHOSRecPoint * PpsdUp) ; AliPHOSTrackSegment(const AliPHOSTrackSegment & ts) ; // ctor virtual ~AliPHOSTrackSegment() { } @@ -39,7 +38,6 @@ public: 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; @@ -52,8 +50,7 @@ public: 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 diff --git a/PHOS/AliPHOSTrackSegmentMaker.h b/PHOS/AliPHOSTrackSegmentMaker.h index d1abd03ef69..1a8efd570c4 100644 --- a/PHOS/AliPHOSTrackSegmentMaker.h +++ b/PHOS/AliPHOSTrackSegmentMaker.h @@ -44,7 +44,6 @@ public: // 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 ; diff --git a/PHOS/AliPHOSTrackSegmentMakerv1.cxx b/PHOS/AliPHOSTrackSegmentMakerv1.cxx index d2c3d911fa4..74f92e95d8e 100644 --- a/PHOS/AliPHOSTrackSegmentMakerv1.cxx +++ b/PHOS/AliPHOSTrackSegmentMakerv1.cxx @@ -60,7 +60,6 @@ #include "AliPHOSClusterizerv1.h" #include "AliPHOSTrackSegment.h" #include "AliPHOSCpvRecPoint.h" -#include "AliPHOSPpsdRecPoint.h" #include "AliPHOSLink.h" #include "AliPHOSGetter.h" #include "AliPHOS.h" @@ -79,9 +78,6 @@ ClassImp( AliPHOSTrackSegmentMakerv1) fEmcLast = 0 ; fCpvFirst = 0 ; fCpvLast = 0 ; - fPpsdFirst = 0 ; - fPpsdLast = 0 ; - fLinkLowArray = 0 ; fLinkUpArray = 0 ; fHeaderFileName = "" ; fRecPointsBranchTitle = "" ; @@ -99,19 +95,17 @@ ClassImp( AliPHOSTrackSegmentMakerv1) 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() ; } @@ -120,7 +114,6 @@ ClassImp( AliPHOSTrackSegmentMakerv1) AliPHOSTrackSegmentMakerv1::~AliPHOSTrackSegmentMakerv1() { // dtor - if(fLinkLowArray) delete fLinkLowArray ; if(fLinkUpArray) delete fLinkUpArray ; } @@ -129,14 +122,12 @@ void AliPHOSTrackSegmentMakerv1::FillOneModule() { // 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) && @@ -146,30 +137,10 @@ void AliPHOSTrackSegmentMakerv1::FillOneModule() //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 ++) ; - } - + } //____________________________________________________________________________ @@ -215,24 +186,23 @@ void AliPHOSTrackSegmentMakerv1::Init() 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) ; } @@ -242,41 +212,25 @@ void AliPHOSTrackSegmentMakerv1::MakeLinks()const // 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++ ) { @@ -291,8 +245,7 @@ void AliPHOSTrackSegmentMakerv1::MakeLinks()const } } - fLinkLowArray->Sort() ; //first links with smallest distances - fLinkUpArray->Sort() ; + fLinkUpArray->Sort() ; //first links with smallest distances } //____________________________________________________________________________ @@ -310,7 +263,7 @@ void AliPHOSTrackSegmentMakerv1::MakePairs() 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) @@ -326,60 +279,29 @@ void AliPHOSTrackSegmentMakerv1::MakePairs() for(index = 0; index fPpsdFirst) - ppsdExist = new Bool_t[fPpsdLast-fPpsdFirst] ; - for(index = 0; index 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 } } @@ -389,9 +311,9 @@ void AliPHOSTrackSegmentMakerv1::MakePairs() 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++; } @@ -416,51 +338,48 @@ void AliPHOSTrackSegmentMakerv1::Exec(Option_t * option) 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() ; @@ -506,93 +425,6 @@ void AliPHOSTrackSegmentMakerv1::Print(Option_t * option)const 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) @@ -606,10 +438,11 @@ 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 @@ -624,7 +457,7 @@ void AliPHOSTrackSegmentMakerv1::WriteTrackSegments(Int_t event) //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()); @@ -640,7 +473,7 @@ void AliPHOSTrackSegmentMakerv1::WriteTrackSegments(Int_t event) 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()); @@ -668,6 +501,7 @@ void AliPHOSTrackSegmentMakerv1::PrintTrackSegments(Option_t * option) taskName.Remove(taskName.Index(Version())-1) ; TClonesArray * trackSegments = AliPHOSGetter::GetInstance()->TrackSegments(taskName) ; + cout << "AliPHOSTrackSegmentMakerv1: event "<GetEvNumber() << endl ; cout << " Found " << trackSegments->GetEntriesFast() << " trackSegments " << endl ; @@ -675,15 +509,14 @@ void AliPHOSTrackSegmentMakerv1::PrintTrackSegments(Option_t * option) 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 GetEntriesFast() ; index++) { AliPHOSTrackSegment * ts = (AliPHOSTrackSegment * )trackSegments->At(index) ; cout<<" "<< setw(4) << ts->GetIndexInList() << " " <GetEmcIndex()<< " " - <GetCpvIndex()<< " " - <GetPpsdIndex()<< endl ; + <GetCpvIndex()<< " " << endl ; } cout << "-------------------------------------------------------"<< endl ; diff --git a/PHOS/AliPHOSTrackSegmentMakerv1.h b/PHOS/AliPHOSTrackSegmentMakerv1.h index cbc3391cfae..5a880f96baa 100644 --- a/PHOS/AliPHOSTrackSegmentMakerv1.h +++ b/PHOS/AliPHOSTrackSegmentMakerv1.h @@ -37,8 +37,8 @@ public: 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) ; @@ -63,7 +63,6 @@ private: 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: @@ -75,16 +74,12 @@ 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 diff --git a/PHOS/AliPHOSv1.cxx b/PHOS/AliPHOSv1.cxx index 8ed90e9b2c4..442da937848 100644 --- a/PHOS/AliPHOSv1.cxx +++ b/PHOS/AliPHOSv1.cxx @@ -69,16 +69,20 @@ AliPHOSv0() 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 @@ -99,6 +103,24 @@ AliPHOSv1::AliPHOSv1(const char *name, const char *title): 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 @@ -148,8 +170,7 @@ AliPHOSv1::~AliPHOSv1() 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 ; @@ -161,7 +182,8 @@ void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id, 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 ; @@ -171,14 +193,13 @@ void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id, 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++ ; } @@ -235,16 +256,14 @@ void AliPHOSv1::FinishEvent() } } } - //____________________________________________________________________________ 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 ; @@ -252,171 +271,153 @@ void AliPHOSv1::StepManager(void) 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; idigitUncheckedAt(idigit); - Float_t x1 = cpvDigit1->GetXpad() ; - Float_t z1 = cpvDigit1->GetYpad() ; - for (Int_t jdigit=idigit+1; jdigitUncheckedAt(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; idigitUncheckedAt(idigit); + Float_t x1 = cpvDigit1->GetXpad() ; + Float_t z1 = cpvDigit1->GetYpad() ; + for (Int_t jdigit=idigit+1; jdigitUncheckedAt(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; idigitUncheckedAt(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; idigitUncheckedAt(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 + } //____________________________________________________________________________ diff --git a/PHOS/AliPHOSv1.h b/PHOS/AliPHOSv1.h index 3fd8a78653d..4b0a76e4aeb 100644 --- a/PHOS/AliPHOSv1.h +++ b/PHOS/AliPHOSv1.h @@ -42,10 +42,7 @@ public: } 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 @@ -57,13 +54,39 @@ public: 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 diff --git a/PHOS/AliPHOSvImpacts.cxx b/PHOS/AliPHOSvImpacts.cxx index 8d888eda690..1414388cb2d 100644 --- a/PHOS/AliPHOSvImpacts.cxx +++ b/PHOS/AliPHOSvImpacts.cxx @@ -43,7 +43,7 @@ ClassImp(AliPHOSvImpacts) //____________________________________________________________________________ -AliPHOSvImpacts::AliPHOSvImpacts() +AliPHOSvImpacts::AliPHOSvImpacts():AliPHOSv1() { // ctor } @@ -53,9 +53,6 @@ AliPHOSvImpacts::AliPHOSvImpacts(const char *name, const char *title): 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 @@ -64,17 +61,15 @@ AliPHOSv1(name,title) // 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; @@ -88,11 +83,6 @@ AliPHOSv1(name,title) fNCPVImpacts[iPHOSModule] = 0; impacts = (TClonesArray *)fCPVImpacts->At(iPHOSModule); } - for (iPHOSModule=0; iPHOSModuleAdd(new TClonesArray("AliPHOSImpact",200)) ; - fNPPSDImpacts[iPHOSModule] = 0; - impacts = (TClonesArray *)fPPSDImpacts->At(iPHOSModule); - } } @@ -108,7 +98,7 @@ AliPHOSvImpacts::~AliPHOSvImpacts() fHits = 0 ; } - // Delete impacts in EMC, CPV and PPSD + // Delete impacts in EMC, CPV if ( fEMCImpacts ) { fEMCImpacts->Delete() ; delete fEMCImpacts ; @@ -119,11 +109,6 @@ AliPHOSvImpacts::~AliPHOSvImpacts() delete fCPVImpacts ; fCPVImpacts = 0 ; } - if ( fPPSDImpacts ) { - fPPSDImpacts->Delete() ; - delete fPPSDImpacts ; - fPPSDImpacts = 0 ; - } } //____________________________________________________________________________ @@ -145,11 +130,6 @@ void AliPHOSvImpacts::AddImpact( char* det, Int_t shunt, Int_t primary, Int_t tr 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) ; @@ -171,7 +151,6 @@ void AliPHOSvImpacts::MakeBranch(Option_t *opt, const char *file) 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); } @@ -188,18 +167,9 @@ void AliPHOSvImpacts::ResetHits() fNEMCImpacts[i] = 0 ; } - if ( strcmp(GetGeometry()->GetName(),"IHEP") == 0 || strcmp(GetGeometry()->GetName(),"MIXT") == 0 ) { - for (i=0; iGetNCPVModules(); i++) { - ((TClonesArray*)fCPVImpacts->At(i)) -> Clear(); - fNCPVImpacts[i] = 0 ; - } - } - - if ( strcmp(GetGeometry()->GetName(),"GPS2") == 0 || strcmp(GetGeometry()->GetName(),"MIXT") == 0 ) { - for (i=0; iGetNPPSDModules(); i++) { - ((TClonesArray*)fPPSDImpacts->At(i)) -> Clear(); - fNPPSDImpacts[i] = 0 ; - } + for (i=0; iGetNModules(); i++) { + ((TClonesArray*)fCPVImpacts->At(i)) -> Clear(); + fNCPVImpacts[i] = 0 ; } } @@ -207,7 +177,7 @@ void AliPHOSvImpacts::ResetHits() //_____________________________________________________________________________ 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(); @@ -244,8 +214,6 @@ void AliPHOSvImpacts::StepManager(void) 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); } @@ -253,8 +221,7 @@ void AliPHOSvImpacts::StepManager(void) // 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) ; @@ -272,33 +239,5 @@ void AliPHOSvImpacts::StepManager(void) 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); - } - } + } diff --git a/PHOS/Makefile b/PHOS/Makefile index 8a5b3e3fe2e..52816396e65 100644 --- a/PHOS/Makefile +++ b/PHOS/Makefile @@ -9,17 +9,17 @@ PACKAGE = PHOS # 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 \ @@ -76,20 +76,7 @@ include $(ALICE_ROOT)/conf/GeneralMacros -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 + -- 2.43.0