From 037cc66d49267fe093d36d6fb90842312fc6f407 Mon Sep 17 00:00:00 2001 From: schutz Date: Mon, 5 Feb 2001 09:03:17 +0000 Subject: [PATCH] parameters have been redistributed; Hits2SDigits etc ... introduce --- PHOS/AliPHOS.cxx | 12 + PHOS/AliPHOS.h | 2 +- PHOS/AliPHOSAnalyze.cxx | 797 ++++++++++++++++++++-------------------- PHOS/AliPHOSDigit.cxx | 2 +- PHOS/AliPHOSDigit.h | 2 +- PHOS/AliPHOSHit.cxx | 29 +- PHOS/AliPHOSHit.h | 16 +- PHOS/AliPHOSv1.cxx | 518 ++++++++++++++------------ PHOS/AliPHOSv1.h | 36 +- PHOS/AliPHOSv2.cxx | 4 +- PHOS/AliPHOSv2.h | 2 +- PHOS/AliPHOSv3.cxx | 264 +++++++++---- 12 files changed, 934 insertions(+), 750 deletions(-) diff --git a/PHOS/AliPHOS.cxx b/PHOS/AliPHOS.cxx index 0c180e494a2..af959462ba5 100644 --- a/PHOS/AliPHOS.cxx +++ b/PHOS/AliPHOS.cxx @@ -32,6 +32,7 @@ class TFile; // --- Standard library --- +#include // --- AliRoot header files --- #include "AliPHOS.h" @@ -348,11 +349,22 @@ void AliPHOS::CreateMaterials() //____________________________________________________________________________ void AliPHOS::SetTreeAddress() { + // TBranch *branch; AliDetector::SetTreeAddress(); TBranch * branch ; + if(fSDigits) + fSDigits->Clear(); + else + fSDigits = new TClonesArray("AliPHOSDigit",1) ; + + if (gAlice->TreeS() && fSDigits ) { + branch = gAlice->TreeS()->GetBranch("PHOS"); + if (branch) branch->SetAddress(&fSDigits) ; + } + if(fDigits) fDigits->Clear(); else diff --git a/PHOS/AliPHOS.h b/PHOS/AliPHOS.h index 60782a82540..1761e48f985 100644 --- a/PHOS/AliPHOS.h +++ b/PHOS/AliPHOS.h @@ -70,7 +70,7 @@ class AliPHOS : public AliDetector { } protected: - + TClonesArray *fSDigits ; // List of summable digits AliPHOSRecPoint::RecPointsList *fEmcRecPoints ; // The RecPoints (clusters) list in EMC AliPHOSRecPoint::RecPointsList *fPpsdRecPoints ;// The RecPoints (clusters) list in PPSD (veto) AliPHOSTrackSegment::TrackSegmentsList *fTrackSegments ;// The TrackSegment list in PHOS diff --git a/PHOS/AliPHOSAnalyze.cxx b/PHOS/AliPHOSAnalyze.cxx index ab1c91b15fc..413a5c4b04c 100644 --- a/PHOS/AliPHOSAnalyze.cxx +++ b/PHOS/AliPHOSAnalyze.cxx @@ -325,10 +325,9 @@ void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod){ (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0) cout << "======= Analyze ======> Event " << ievent+1 << endl ; - //=========== Connects the various Tree's for evt - Int_t tracks = gAlice->GetEvent(ievent); + fPHOS->Enable() ; - fPHOS->Hit2Digit(tracks) ; + gAlice->Hits2Digits() ; //=========== Do the reconstruction fPHOS->Reconstruction(fRec); @@ -428,183 +427,183 @@ void AliPHOSAnalyze::ReadAndPrintCPV(Int_t EvFirst, Int_t EvLast) //____________________________________________________________________________ void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents) { - // - // Analyzes CPV characteristics - // Author: Yuri Kharlov - // 9 October 2000 - // - - // Book histograms - - TH1F *hDx = new TH1F("hDx" ,"CPV x-resolution@reconstruction",100,-5. , 5.); - TH1F *hDz = new TH1F("hDz" ,"CPV z-resolution@reconstruction",100,-5. , 5.); - TH1F *hDr = new TH1F("hDr" ,"CPV r-resolution@reconstruction",100, 0. , 5.); - TH1S *hNrp = new TH1S("hNrp" ,"CPV rec.point multiplicity", 21,-0.5,20.5); - TH1S *hNrpX = new TH1S("hNrpX","CPV rec.point Phi-length" , 21,-0.5,20.5); - TH1S *hNrpZ = new TH1S("hNrpZ","CPV rec.point Z-length" , 21,-0.5,20.5); - - cout << "Start CPV Analysis"<< endl ; - for ( Int_t ievent=0; ievent -// if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0) - cout << endl << "==== AnalyzeCPV ====> Event is " << ievent+1 << endl ; +// //========== Event Number> +// // if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0) +// cout << endl << "==== AnalyzeCPV ====> Event is " << ievent+1 << endl ; - //=========== Connects the various Tree's for evt - Int_t ntracks = gAlice->GetEvent(ievent); +// //=========== Connects the various Tree's for evt +// Int_t ntracks = gAlice->GetEvent(ievent); - //========== Creating branches =================================== - AliPHOSRecPoint::RecPointsList ** emcRecPoints = fPHOS->EmcRecPoints() ; - gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , emcRecPoints ) ; +// //========== Creating branches =================================== +// AliPHOSRecPoint::RecPointsList ** emcRecPoints = fPHOS->EmcRecPoints() ; +// gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , emcRecPoints ) ; - AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->PpsdRecPoints() ; - gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", cpvRecPoints ) ; +// AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->PpsdRecPoints() ; +// gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", cpvRecPoints ) ; - // Create and fill arrays of hits for each CPV module +// // Create and fill arrays of hits for each CPV module - Int_t nOfModules = fGeom->GetNModules(); - TClonesArray **hitsPerModule = new TClonesArray *[nOfModules]; - Int_t iModule = 0; - for (iModule=0; iModule < nOfModules; iModule++) - hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100); - - AliPHOSCPVModule cpvModule; - TClonesArray *cpvHits; - Int_t nCPVhits; - AliPHOSCPVHit *cpvHit; - TLorentzVector p; - Float_t xzgen[2]; - Int_t ipart; - - // First go through all primary tracks and fill the arrays - // of hits per each CPV module - - for (Int_t itrack=0; itrackResetHits(); - gAlice->TreeH()->GetEvent(itrack); - for (Int_t iModule=0; iModule < nOfModules; iModule++) { - cpvModule = fPHOS->GetCPVModule(iModule); - cpvHits = cpvModule.Hits(); - nCPVhits = cpvHits->GetEntriesFast(); - for (Int_t ihit=0; ihitUncheckedAt(ihit); - p = cpvHit->GetMomentum(); - xzgen[0] = cpvHit->X(); - xzgen[1] = cpvHit->Y(); - ipart = cpvHit->GetIpart(); - TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule]; - new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*cpvHit); - } - cpvModule.Clear(); - } - } - for (iModule=0; iModule < nOfModules; iModule++) { - Int_t nsum = hitsPerModule[iModule]->GetEntriesFast(); - printf("Module %d has %d hits\n",iModule,nsum); - } +// Int_t nOfModules = fGeom->GetNModules(); +// TClonesArray **hitsPerModule = new TClonesArray *[nOfModules]; +// Int_t iModule = 0; +// for (iModule=0; iModule < nOfModules; iModule++) +// hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100); - // Then go through reconstructed points and for each find - // the closeset hit - // The distance from the rec.point to the closest hit - // gives the coordinate resolution of the CPV - - // Get the Reconstruction Tree - gAlice->TreeR()->GetEvent(0) ; - TIter nextRP(*fPHOS->PpsdRecPoints() ) ; - AliPHOSCpvRecPoint *cpvRecPoint ; - Float_t xgen, zgen; - while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) { - TVector3 locpos; - cpvRecPoint->GetLocalPosition(locpos); - Int_t phosModule = cpvRecPoint->GetPHOSMod(); - Int_t rpMult = cpvRecPoint->GetDigitsMultiplicity(); - Int_t rpMultX, rpMultZ; - cpvRecPoint->GetClusterLengths(rpMultX,rpMultZ); - Float_t xrec = locpos.X(); - Float_t zrec = locpos.Z(); - Float_t dxmin = 1.e+10; - Float_t dzmin = 1.e+10; - Float_t r2min = 1.e+10; - Float_t r2; - - cpvHits = hitsPerModule[phosModule-1]; - Int_t nCPVhits = cpvHits->GetEntriesFast(); - for (Int_t ihit=0; ihitUncheckedAt(ihit); - xgen = cpvHit->X(); - zgen = cpvHit->Y(); - r2 = TMath::Power((xgen-xrec),2) + TMath::Power((zgen-zrec),2); - if ( r2 < r2min ) { - r2min = r2; - dxmin = xgen - xrec; - dzmin = zgen - zrec; - } - } - hDx ->Fill(dxmin); - hDz ->Fill(dzmin); - hDr ->Fill(TMath::Sqrt(r2min)); - hNrp ->Fill(rpMult); - hNrpX->Fill(rpMultX); - hNrpZ->Fill(rpMultZ); - } - delete [] hitsPerModule; - } - // Save histograms +// AliPHOSCPVModule cpvModule; +// TClonesArray *cpvHits; +// Int_t nCPVhits; +// AliPHOSCPVHit *cpvHit; +// TLorentzVector p; +// Float_t xzgen[2]; +// Int_t ipart; - Text_t outputname[80] ; - sprintf(outputname,"%s.analyzed",fRootFile->GetName()); - TFile output(outputname,"RECREATE"); - output.cd(); +// // First go through all primary tracks and fill the arrays +// // of hits per each CPV module - hDx ->Write() ; - hDz ->Write() ; - hDr ->Write() ; - hNrp ->Write() ; - hNrpX->Write() ; - hNrpZ->Write() ; - - // Plot histograms - - TCanvas *cpvCanvas = new TCanvas("CPV","CPV analysis",20,20,800,400); - gStyle->SetOptStat(111111); - gStyle->SetOptFit(1); - gStyle->SetOptDate(1); - cpvCanvas->Divide(3,2); - - cpvCanvas->cd(1); - gPad->SetFillColor(10); - hNrp->SetFillColor(16); - hNrp->Draw(); - - cpvCanvas->cd(2); - gPad->SetFillColor(10); - hNrpX->SetFillColor(16); - hNrpX->Draw(); - - cpvCanvas->cd(3); - gPad->SetFillColor(10); - hNrpZ->SetFillColor(16); - hNrpZ->Draw(); - - cpvCanvas->cd(4); - gPad->SetFillColor(10); - hDx->SetFillColor(16); - hDx->Fit("gaus"); - hDx->Draw(); - - cpvCanvas->cd(5); - gPad->SetFillColor(10); - hDz->SetFillColor(16); - hDz->Fit("gaus"); - hDz->Draw(); - - cpvCanvas->cd(6); - gPad->SetFillColor(10); - hDr->SetFillColor(16); - hDr->Draw(); - - cpvCanvas->Print("CPV.ps"); +// for (Int_t itrack=0; itrackResetHits(); +// gAlice->TreeH()->GetEvent(itrack); +// for (Int_t iModule=0; iModule < nOfModules; iModule++) { +// cpvModule = fPHOS->GetCPVModule(iModule); +// cpvHits = cpvModule.Hits(); +// nCPVhits = cpvHits->GetEntriesFast(); +// for (Int_t ihit=0; ihitUncheckedAt(ihit); +// p = cpvHit->GetMomentum(); +// xzgen[0] = cpvHit->X(); +// xzgen[1] = cpvHit->Y(); +// ipart = cpvHit->GetIpart(); +// TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule]; +// new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*cpvHit); +// } +// cpvModule.Clear(); +// } +// } +// for (iModule=0; iModule < nOfModules; iModule++) { +// Int_t nsum = hitsPerModule[iModule]->GetEntriesFast(); +// printf("Module %d has %d hits\n",iModule,nsum); +// } + +// // Then go through reconstructed points and for each find +// // the closeset hit +// // The distance from the rec.point to the closest hit +// // gives the coordinate resolution of the CPV + +// // Get the Reconstruction Tree +// gAlice->TreeR()->GetEvent(0) ; +// TIter nextRP(*fPHOS->PpsdRecPoints() ) ; +// AliPHOSCpvRecPoint *cpvRecPoint ; +// Float_t xgen, zgen; +// while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) { +// TVector3 locpos; +// cpvRecPoint->GetLocalPosition(locpos); +// Int_t phosModule = cpvRecPoint->GetPHOSMod(); +// Int_t rpMult = cpvRecPoint->GetDigitsMultiplicity(); +// Int_t rpMultX, rpMultZ; +// cpvRecPoint->GetClusterLengths(rpMultX,rpMultZ); +// Float_t xrec = locpos.X(); +// Float_t zrec = locpos.Z(); +// Float_t dxmin = 1.e+10; +// Float_t dzmin = 1.e+10; +// Float_t r2min = 1.e+10; +// Float_t r2; + +// cpvHits = hitsPerModule[phosModule-1]; +// Int_t nCPVhits = cpvHits->GetEntriesFast(); +// for (Int_t ihit=0; ihitUncheckedAt(ihit); +// xgen = cpvHit->X(); +// zgen = cpvHit->Y(); +// r2 = TMath::Power((xgen-xrec),2) + TMath::Power((zgen-zrec),2); +// if ( r2 < r2min ) { +// r2min = r2; +// dxmin = xgen - xrec; +// dzmin = zgen - zrec; +// } +// } +// hDx ->Fill(dxmin); +// hDz ->Fill(dzmin); +// hDr ->Fill(TMath::Sqrt(r2min)); +// hNrp ->Fill(rpMult); +// hNrpX->Fill(rpMultX); +// hNrpZ->Fill(rpMultZ); +// } +// delete [] hitsPerModule; +// } +// // Save histograms + +// Text_t outputname[80] ; +// sprintf(outputname,"%s.analyzed",fRootFile->GetName()); +// TFile output(outputname,"RECREATE"); +// output.cd(); + +// hDx ->Write() ; +// hDz ->Write() ; +// hDr ->Write() ; +// hNrp ->Write() ; +// hNrpX->Write() ; +// hNrpZ->Write() ; + +// // Plot histograms + +// TCanvas *cpvCanvas = new TCanvas("CPV","CPV analysis",20,20,800,400); +// gStyle->SetOptStat(111111); +// gStyle->SetOptFit(1); +// gStyle->SetOptDate(1); +// cpvCanvas->Divide(3,2); + +// cpvCanvas->cd(1); +// gPad->SetFillColor(10); +// hNrp->SetFillColor(16); +// hNrp->Draw(); + +// cpvCanvas->cd(2); +// gPad->SetFillColor(10); +// hNrpX->SetFillColor(16); +// hNrpX->Draw(); + +// cpvCanvas->cd(3); +// gPad->SetFillColor(10); +// hNrpZ->SetFillColor(16); +// hNrpZ->Draw(); + +// cpvCanvas->cd(4); +// gPad->SetFillColor(10); +// hDx->SetFillColor(16); +// hDx->Fit("gaus"); +// hDx->Draw(); + +// cpvCanvas->cd(5); +// gPad->SetFillColor(10); +// hDz->SetFillColor(16); +// hDz->Fit("gaus"); +// hDz->Draw(); + +// cpvCanvas->cd(6); +// gPad->SetFillColor(10); +// hDr->SetFillColor(16); +// hDr->Draw(); + +// cpvCanvas->Print("CPV.ps"); } @@ -723,257 +722,257 @@ void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents) //____________________________________________________________________________ void AliPHOSAnalyze::ReadAndPrintEMC(Int_t EvFirst, Int_t EvLast) { - // - // Read and print generated and reconstructed hits in EMC - // for events from EvFirst to Nevent. - // If only EvFirst is defined, print only this one event. - // Author: Yuri Kharlov - // 24 November 2000 - // - - if (EvFirst!=0 && EvLast==0) EvLast=EvFirst; - Int_t ievent; - for (ievent=EvFirst; ievent<=EvLast; ievent++) { +// // +// // Read and print generated and reconstructed hits in EMC +// // for events from EvFirst to Nevent. +// // If only EvFirst is defined, print only this one event. +// // Author: Yuri Kharlov +// // 24 November 2000 +// // + +// if (EvFirst!=0 && EvLast==0) EvLast=EvFirst; +// Int_t ievent; +// for (ievent=EvFirst; ievent<=EvLast; ievent++) { - //========== Event Number> - cout << endl << "==== ReadAndPrintEMC ====> Event is " << ievent+1 << endl ; +// //========== Event Number> +// cout << endl << "==== ReadAndPrintEMC ====> Event is " << ievent+1 << endl ; - //=========== Connects the various Tree's for evt - Int_t ntracks = gAlice->GetEvent(ievent); - fPHOS->SetTreeAddress() ; +// //=========== Connects the various Tree's for evt +// Int_t ntracks = gAlice->GetEvent(ievent); +// fPHOS->SetTreeAddress() ; - gAlice->TreeD()->GetEvent(0) ; - gAlice->TreeR()->GetEvent(0) ; +// gAlice->TreeD()->GetEvent(0) ; +// gAlice->TreeR()->GetEvent(0) ; - // Loop over reconstructed particles +// // Loop over reconstructed particles - TClonesArray ** recParticleList = fPHOS->RecParticles() ; - AliPHOSRecParticle * recParticle ; - Int_t iRecParticle ; - Int_t *primList; - Int_t nPrimary; - for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ ) { - recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ; - Float_t recE = recParticle->Energy(); - primList = recParticle->GetPrimaries(nPrimary); - Int_t moduleNumberRec ; - Double_t recX, recZ ; - fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ; - printf("Rec point: module %d, (X,Z) = (%8.4f,%8.4f) cm, E = %.3f GeV, primary = %d\n", - moduleNumberRec,recX,recZ,recE,*primList); - } +// TClonesArray ** recParticleList = fPHOS->RecParticles() ; +// AliPHOSRecParticle * recParticle ; +// Int_t iRecParticle ; +// Int_t *primList; +// Int_t nPrimary; +// for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ ) { +// recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ; +// Float_t recE = recParticle->Energy(); +// primList = recParticle->GetPrimaries(nPrimary); +// Int_t moduleNumberRec ; +// Double_t recX, recZ ; +// fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ; +// printf("Rec point: module %d, (X,Z) = (%8.4f,%8.4f) cm, E = %.3f GeV, primary = %d\n", +// moduleNumberRec,recX,recZ,recE,*primList); +// } - // Read and print EMC hits from EMCn branches +// // Read and print EMC hits from EMCn branches - AliPHOSCPVModule emcModule; - TClonesArray *emcHits; - Int_t nEMChits; - AliPHOSCPVHit *emcHit; - TLorentzVector p; - Float_t xgen, zgen; - Int_t ipart, primary; - Int_t nGenHits = 0; - for (Int_t itrack=0; itrackResetHits(); - gAlice->TreeH()->GetEvent(itrack); - Int_t iModule = 0 ; - for (iModule=0; iModule < fGeom->GetNModules(); iModule++) { - emcModule = fPHOS->GetEMCModule(iModule); - emcHits = emcModule.Hits(); - nEMChits = emcHits->GetEntriesFast(); - for (Int_t ihit=0; ihitUncheckedAt(ihit); - p = emcHit->GetMomentum(); - xgen = emcHit->X(); - zgen = emcHit->Y(); - ipart = emcHit->GetIpart(); - primary= emcHit->GetTrack(); - printf("EMC hit A: module %d, ",iModule+1); - printf(" p = (%f .4, %f .4, %f .4, %f .4) GeV,\n", - p.Px(),p.Py(),p.Pz(),p.Energy()); - printf(" (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n", - xgen,zgen,ipart,primary); - } - } - } +// AliPHOSCPVModule emcModule; +// TClonesArray *emcHits; +// Int_t nEMChits; +// AliPHOSCPVHit *emcHit; +// TLorentzVector p; +// Float_t xgen, zgen; +// Int_t ipart, primary; +// Int_t nGenHits = 0; +// for (Int_t itrack=0; itrackResetHits(); +// gAlice->TreeH()->GetEvent(itrack); +// Int_t iModule = 0 ; +// for (iModule=0; iModule < fGeom->GetNModules(); iModule++) { +// emcModule = fPHOS->GetEMCModule(iModule); +// emcHits = emcModule.Hits(); +// nEMChits = emcHits->GetEntriesFast(); +// for (Int_t ihit=0; ihitUncheckedAt(ihit); +// p = emcHit->GetMomentum(); +// xgen = emcHit->X(); +// zgen = emcHit->Y(); +// ipart = emcHit->GetIpart(); +// primary= emcHit->GetTrack(); +// printf("EMC hit A: module %d, ",iModule+1); +// printf(" p = (%f .4, %f .4, %f .4, %f .4) GeV,\n", +// p.Px(),p.Py(),p.Pz(),p.Energy()); +// printf(" (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n", +// xgen,zgen,ipart,primary); +// } +// } +// } -// // Read and print EMC hits from PHOS branch - -// for (Int_t itrack=0; itrackResetHits(); -// gAlice->TreeH()->GetEvent(itrack); -// TClonesArray *hits = fPHOS->Hits(); -// AliPHOSHit *hit ; -// Int_t ihit; -// for ( ihit = 0 ; ihit < hits->GetEntries() ; ihit++ ) { -// hit = (AliPHOSHit*)hits->At(ihit) ; -// Float_t hitXYZ[3]; -// hitXYZ[0] = hit->X(); -// hitXYZ[1] = hit->Y(); -// hitXYZ[2] = hit->Z(); -// ipart = hit->GetPid(); -// primary = hit->GetPrimary(); -// Int_t absId = hit->GetId(); -// Int_t relId[4]; -// fGeom->AbsToRelNumbering(absId, relId) ; -// Int_t module = relId[0]; -// if (relId[1]==0 && !(hitXYZ[0]==0 && hitXYZ[2]==0)) -// printf("EMC hit B: module %d, (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n", -// module,hitXYZ[0],hitXYZ[2],ipart,primary); -// } -// } +// // // Read and print EMC hits from PHOS branch + +// // for (Int_t itrack=0; itrackResetHits(); +// // gAlice->TreeH()->GetEvent(itrack); +// // TClonesArray *hits = fPHOS->Hits(); +// // AliPHOSHit *hit ; +// // Int_t ihit; +// // for ( ihit = 0 ; ihit < hits->GetEntries() ; ihit++ ) { +// // hit = (AliPHOSHit*)hits->At(ihit) ; +// // Float_t hitXYZ[3]; +// // hitXYZ[0] = hit->X(); +// // hitXYZ[1] = hit->Y(); +// // hitXYZ[2] = hit->Z(); +// // ipart = hit->GetPid(); +// // primary = hit->GetPrimary(); +// // Int_t absId = hit->GetId(); +// // Int_t relId[4]; +// // fGeom->AbsToRelNumbering(absId, relId) ; +// // Int_t module = relId[0]; +// // if (relId[1]==0 && !(hitXYZ[0]==0 && hitXYZ[2]==0)) +// // printf("EMC hit B: module %d, (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n", +// // module,hitXYZ[0],hitXYZ[2],ipart,primary); +// // } +// // } - } +// } } //____________________________________________________________________________ void AliPHOSAnalyze::AnalyzeEMC(Int_t Nevents) { - // - // Read generated and reconstructed hits in EMC for Nevents events. - // Plots the coordinate and energy resolution histograms. - // Coordinate resolution is a difference between the reconstructed - // coordinate and the exact coordinate on the face of the PHOS - // Author: Yuri Kharlov - // 27 November 2000 - // - - // Book histograms - - TH1F *hDx1 = new TH1F("hDx1" ,"EMC x-resolution", 100,-5. , 5.); - TH1F *hDz1 = new TH1F("hDz1" ,"EMC z-resolution", 100,-5. , 5.); - TH1F *hDE1 = new TH1F("hDE1" ,"EMC E-resolution", 100,-2. , 2.); - - TH2F *hDx2 = new TH2F("hDx2" ,"EMC x-resolution", 100, 0., 10., 100,-5. , 5.); - TH2F *hDz2 = new TH2F("hDz2" ,"EMC z-resolution", 100, 0., 10., 100,-5. , 5.); - TH2F *hDE2 = new TH2F("hDE2" ,"EMC E-resolution", 100, 0., 10., 100, 0. , 5.); - - cout << "Start EMC Analysis"<< endl ; - for (Int_t ievent=0; ievent - if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0) - cout << "==== AnalyzeEMC ====> Event is " << ievent+1 << endl ; +// //========== Event Number> +// if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0) +// cout << "==== AnalyzeEMC ====> Event is " << ievent+1 << endl ; - //=========== Connects the various Tree's for evt - Int_t ntracks = gAlice->GetEvent(ievent); +// //=========== Connects the various Tree's for evt +// Int_t ntracks = gAlice->GetEvent(ievent); - fPHOS->SetTreeAddress() ; +// fPHOS->SetTreeAddress() ; - gAlice->TreeD()->GetEvent(0) ; - gAlice->TreeR()->GetEvent(0) ; +// gAlice->TreeD()->GetEvent(0) ; +// gAlice->TreeR()->GetEvent(0) ; - // Create and fill arrays of hits for each EMC module +// // Create and fill arrays of hits for each EMC module - Int_t nOfModules = fGeom->GetNModules(); - TClonesArray **hitsPerModule = new TClonesArray *[nOfModules]; - Int_t iModule; - for (iModule=0; iModule < nOfModules; iModule++) - hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100); - - AliPHOSCPVModule emcModule; - TClonesArray *emcHits; - Int_t nEMChits; - AliPHOSCPVHit *emcHit; - - // First go through all primary tracks and fill the arrays - // of hits per each EMC module - - for (Int_t itrack=0; itrackResetHits(); - gAlice->TreeH()->GetEvent(itrack); - for (Int_t iModule=0; iModule < nOfModules; iModule++) { - emcModule = fPHOS->GetEMCModule(iModule); - emcHits = emcModule.Hits(); - nEMChits = emcHits->GetEntriesFast(); - for (Int_t ihit=0; ihitUncheckedAt(ihit); - TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule]; - new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*emcHit); - } - emcModule.Clear(); - } - } +// Int_t nOfModules = fGeom->GetNModules(); +// TClonesArray **hitsPerModule = new TClonesArray *[nOfModules]; +// Int_t iModule; +// for (iModule=0; iModule < nOfModules; iModule++) +// hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100); - // Loop over reconstructed particles - - TClonesArray ** recParticleList = fPHOS->RecParticles() ; - AliPHOSRecParticle * recParticle ; - Int_t nEMCrecs = (*recParticleList)->GetEntries(); - if (nEMCrecs == 1) { - recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(0) ; - Float_t recE = recParticle->Energy(); - Int_t phosModule; - Double_t recX, recZ ; - fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), phosModule, recX, recZ) ; - - // for this rec.point take the hit list in the same PHOS module - - emcHits = hitsPerModule[phosModule-1]; - Int_t nEMChits = emcHits->GetEntriesFast(); - if (nEMChits == 1) { - Float_t genX, genZ, genE; - for (Int_t ihit=0; ihitUncheckedAt(ihit); - genX = emcHit->X(); - genZ = emcHit->Y(); - genE = emcHit->GetMomentum().E(); - } - Float_t dx = recX - genX; - Float_t dz = recZ - genZ; - Float_t de = recE - genE; - hDx1 ->Fill(dx); - hDz1 ->Fill(dz); - hDE1 ->Fill(de); - hDx2 ->Fill(genE,dx); - hDz2 ->Fill(genE,dz); - hDE2 ->Fill(genE,recE); - } - } - delete [] hitsPerModule; - } - // Save histograms +// AliPHOSCPVModule emcModule; +// TClonesArray *emcHits; +// Int_t nEMChits; +// AliPHOSCPVHit *emcHit; - Text_t outputname[80] ; - sprintf(outputname,"%s.analyzed",fRootFile->GetName()); - TFile output(outputname,"RECREATE"); - output.cd(); +// // First go through all primary tracks and fill the arrays +// // of hits per each EMC module - hDx1 ->Write() ; - hDz1 ->Write() ; - hDE1 ->Write() ; - hDx2 ->Write() ; - hDz2 ->Write() ; - hDE2 ->Write() ; - - // Plot histograms - - TCanvas *emcCanvas = new TCanvas("EMC","EMC analysis",20,20,700,300); - gStyle->SetOptStat(111111); - gStyle->SetOptFit(1); - gStyle->SetOptDate(1); - emcCanvas->Divide(3,1); - - emcCanvas->cd(1); - gPad->SetFillColor(10); - hDx1->SetFillColor(16); - hDx1->Draw(); - - emcCanvas->cd(2); - gPad->SetFillColor(10); - hDz1->SetFillColor(16); - hDz1->Draw(); - - emcCanvas->cd(3); - gPad->SetFillColor(10); - hDE1->SetFillColor(16); - hDE1->Draw(); - - emcCanvas->Print("EMC.ps"); +// for (Int_t itrack=0; itrackResetHits(); +// gAlice->TreeH()->GetEvent(itrack); +// for (Int_t iModule=0; iModule < nOfModules; iModule++) { +// emcModule = fPHOS->GetEMCModule(iModule); +// emcHits = emcModule.Hits(); +// nEMChits = emcHits->GetEntriesFast(); +// for (Int_t ihit=0; ihitUncheckedAt(ihit); +// TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule]; +// new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*emcHit); +// } +// emcModule.Clear(); +// } +// } + +// // Loop over reconstructed particles + +// TClonesArray ** recParticleList = fPHOS->RecParticles() ; +// AliPHOSRecParticle * recParticle ; +// Int_t nEMCrecs = (*recParticleList)->GetEntries(); +// if (nEMCrecs == 1) { +// recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(0) ; +// Float_t recE = recParticle->Energy(); +// Int_t phosModule; +// Double_t recX, recZ ; +// fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), phosModule, recX, recZ) ; + +// // for this rec.point take the hit list in the same PHOS module + +// emcHits = hitsPerModule[phosModule-1]; +// Int_t nEMChits = emcHits->GetEntriesFast(); +// if (nEMChits == 1) { +// Float_t genX, genZ, genE; +// for (Int_t ihit=0; ihitUncheckedAt(ihit); +// genX = emcHit->X(); +// genZ = emcHit->Y(); +// genE = emcHit->GetMomentum().E(); +// } +// Float_t dx = recX - genX; +// Float_t dz = recZ - genZ; +// Float_t de = recE - genE; +// hDx1 ->Fill(dx); +// hDz1 ->Fill(dz); +// hDE1 ->Fill(de); +// hDx2 ->Fill(genE,dx); +// hDz2 ->Fill(genE,dz); +// hDE2 ->Fill(genE,recE); +// } +// } +// delete [] hitsPerModule; +// } +// // Save histograms + +// Text_t outputname[80] ; +// sprintf(outputname,"%s.analyzed",fRootFile->GetName()); +// TFile output(outputname,"RECREATE"); +// output.cd(); + +// hDx1 ->Write() ; +// hDz1 ->Write() ; +// hDE1 ->Write() ; +// hDx2 ->Write() ; +// hDz2 ->Write() ; +// hDE2 ->Write() ; + +// // Plot histograms + +// TCanvas *emcCanvas = new TCanvas("EMC","EMC analysis",20,20,700,300); +// gStyle->SetOptStat(111111); +// gStyle->SetOptFit(1); +// gStyle->SetOptDate(1); +// emcCanvas->Divide(3,1); + +// emcCanvas->cd(1); +// gPad->SetFillColor(10); +// hDx1->SetFillColor(16); +// hDx1->Draw(); + +// emcCanvas->cd(2); +// gPad->SetFillColor(10); +// hDz1->SetFillColor(16); +// hDz1->Draw(); + +// emcCanvas->cd(3); +// gPad->SetFillColor(10); +// hDE1->SetFillColor(16); +// hDE1->Draw(); + +// emcCanvas->Print("EMC.ps"); } diff --git a/PHOS/AliPHOSDigit.cxx b/PHOS/AliPHOSDigit.cxx index c137b5896e7..7ad00684d19 100644 --- a/PHOS/AliPHOSDigit.cxx +++ b/PHOS/AliPHOSDigit.cxx @@ -95,7 +95,7 @@ AliPHOSDigit::~AliPHOSDigit() } //____________________________________________________________________________ -Int_t AliPHOSDigit::Compare(TObject * obj) +Int_t AliPHOSDigit::Compare(const TObject * obj) const { // Compares two digits with respect to its Id // to sort according increasing Id diff --git a/PHOS/AliPHOSDigit.h b/PHOS/AliPHOSDigit.h index 5d4c5ccb2cc..c8bafc6483b 100644 --- a/PHOS/AliPHOSDigit.h +++ b/PHOS/AliPHOSDigit.h @@ -39,7 +39,7 @@ class AliPHOSDigit : public AliDigitNew { Bool_t operator==(const AliPHOSDigit &rValue) const; AliPHOSDigit& operator+(AliPHOSDigit const &rValue) ; - Int_t Compare(TObject * obj) ; + Int_t Compare(const TObject * obj) const ; Int_t GetNprimary() const { // returns the number of primaries return fNprimary ; } diff --git a/PHOS/AliPHOSHit.cxx b/PHOS/AliPHOSHit.cxx index 762689547b0..669553ea341 100644 --- a/PHOS/AliPHOSHit.cxx +++ b/PHOS/AliPHOSHit.cxx @@ -54,27 +54,30 @@ AliPHOSHit::AliPHOSHit(const AliPHOSHit & hit) } //____________________________________________________________________________ -AliPHOSHit::AliPHOSHit(Int_t Shunt, Int_t primary, Int_t Track, Int_t id, Float_t *hits, Int_t pid) : AliHit(Shunt, Track) +AliPHOSHit::AliPHOSHit(Int_t shunt, Int_t primary, Int_t track, Int_t id, Float_t *hits, Int_t pid, TLorentzVector p, Float_t *xy): AliHit(shunt, track) { - // ctor - - fId = id ; - fX = hits[0] ; - fY = hits[1] ; - fZ = hits[2] ; - fELOS = hits[3] ; - fPrimary = primary ; - fPid = pid ; + // + // Create a CPV hit object + // + + fId = id ; + fELOS = hits[3] ; + fPrimary = primary ; + fPid = pid ; + fMomentum = p; + fX = xy[0]; //position of particle first entering cristall/pad + fY = xy[1]; + fZ = 9999.; //Fake Z to avoid FPE } //____________________________________________________________________________ Bool_t AliPHOSHit::operator==(AliPHOSHit const &rValue) const { - // Two hits are identical if they have the same Id and originate from the same primary + // Two hits are identical if they have the same Id and originate from the same primary Bool_t rv = kFALSE ; - if ( fId == rValue.GetId() && fPrimary == rValue.GetPrimary() ) + if ( (fId == rValue.GetId()) && ( fPrimary == rValue.GetPrimary() ) && (fPid*rValue.fPid ! = 0) ) rv = kTRUE; return rv; @@ -94,6 +97,8 @@ AliPHOSHit AliPHOSHit::operator+(const AliPHOSHit &rValue) const added.fELOS += rValue.GetEnergy() ; + if(fPid == 0) fPid = rValue.fPid ; + return added; } diff --git a/PHOS/AliPHOSHit.h b/PHOS/AliPHOSHit.h index 936da55fcc8..f2b0beddc92 100644 --- a/PHOS/AliPHOSHit.h +++ b/PHOS/AliPHOSHit.h @@ -20,6 +20,8 @@ #include +class TLorentzVector ; + class AliPHOSHit : public AliHit { friend ostream& operator << (ostream&, const AliPHOSHit&) ; @@ -30,7 +32,7 @@ class AliPHOSHit : public AliHit { // default ctor } AliPHOSHit(const AliPHOSHit & hit) ; - AliPHOSHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t id, Float_t *hits, Int_t pid) ; + AliPHOSHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t id, Float_t *hits, Int_t pid, TLorentzVector p, Float_t *xy); virtual ~AliPHOSHit(void) { // dtor } @@ -51,6 +53,8 @@ class AliPHOSHit : public AliHit { // returns the primary particle id at the origine of this hit return fPrimary ; } + TLorentzVector GetMomentum() { return fMomentum; } + // momentum of the particle which initiated this hit Bool_t operator == (AliPHOSHit const &rValue) const ; AliPHOSHit operator + (const AliPHOSHit& rValue) const ; @@ -58,10 +62,12 @@ class AliPHOSHit : public AliHit { private: - Int_t fId ; // Absolute Id number of PHOS Xtal or PPSD pad - Float_t fELOS ; // Energy deposited - Int_t fPid ; // type of the particle that initiates that hit - Int_t fPrimary ; // Primary particles at the origine of the hit + Int_t fId ; // Absolute Id number of PHOS Xtal or PPSD pad + Float_t fELOS ; // Energy deposited + Int_t fPid ; // type of the particle that initiates that hit + Int_t fPrimary ; // Primary particles at the origine of the hit + + TLorentzVector fMomentum; // 4-momentum of the particle ClassDef(AliPHOSHit,1) // Hit for PHOS diff --git a/PHOS/AliPHOSv1.cxx b/PHOS/AliPHOSv1.cxx index 30a4bdf9781..6e50ea0da96 100644 --- a/PHOS/AliPHOSv1.cxx +++ b/PHOS/AliPHOSv1.cxx @@ -67,16 +67,18 @@ AliPHOSv1::AliPHOSv1() fReconstructioner = 0; fTrackSegmentMaker = 0; + + fHits = new TClonesArray("AliPHOSHit",1000) ; - if ( 0==(fEMCModules=new TClonesArray("AliPHOSCPVModule",0)) ) { - Error("AliPHOSv1","Can not create array of EMC modules"); - exit(1); - } + // if ( 0==(fEMCModules=new TClonesArray("AliPHOSCPVModule",0)) ) { + // Error("AliPHOSv1","Can not create array of EMC modules"); + // exit(1); + // } - if ( 0==(fCPVModules=new TClonesArray("AliPHOSCPVModule",0)) ) { - Error("AliPHOSv1","Can not create array of CPV modules"); - exit(1); - } + // if ( 0==(fCPVModules=new TClonesArray("AliPHOSCPVModule",0)) ) { + // Error("AliPHOSv1","Can not create array of CPV modules"); + // exit(1); + // } } @@ -96,7 +98,10 @@ AliPHOSv0(name,title) // fPinElectronicNoise = 0.010 ; - fDigitThreshold = 0.1 ; // 1 GeV + fDigitThreshold = 0.01 ; // 1 GeV + fDigitizeA= 0. ; + fDigitizeB = 10000000. ; + // We do not want to save in TreeH the raw hits // But save the cumulated hits instead (need to create the branch myself) @@ -114,32 +119,32 @@ AliPHOSv0(name,title) // Create array of EMC modules of the size of PHOS modules number - if ( 0==(fEMCModules=new TClonesArray("AliPHOSCPVModule",fGeom->GetNModules())) ) { - Error("AliPHOSv1","Can not create array of EMC modules"); - exit(1); - } - TClonesArray &lemcmodule = *fEMCModules; - for (Int_t i=0; iGetNModules(); i++) new(lemcmodule[i]) AliPHOSCPVModule(); +// if ( 0==(fEMCModules=new TClonesArray("AliPHOSCPVModule",fGeom->GetNModules())) ) { +// Error("AliPHOSv1","Can not create array of EMC modules"); +// exit(1); +// } +// TClonesArray &lemcmodule = *fEMCModules; +// for (Int_t i=0; iGetNModules(); i++) new(lemcmodule[i]) AliPHOSCPVModule(); // Create array of CPV modules for the IHEP's version of CPV - if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) { - // Create array of CPV modules of the size of PHOS modules number +// if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) { +// // Create array of CPV modules of the size of PHOS modules number - if ( 0==(fCPVModules=new TClonesArray("AliPHOSCPVModule",fGeom->GetNCPVModules())) ) { - Error("AliPHOSv1","Can not create array of CPV modules"); - exit(1); - } - TClonesArray &lcpvmodule = *fCPVModules; - for (Int_t i=0; iGetNCPVModules(); i++) new(lcpvmodule[i]) AliPHOSCPVModule(); - } - else { - // Create an empty array of AliPHOSCPVModule to satisfy - // AliPHOSv1::Streamer when writing root file +// if ( 0==(fCPVModules=new TClonesArray("AliPHOSCPVModule",fGeom->GetNCPVModules())) ) { +// Error("AliPHOSv1","Can not create array of CPV modules"); +// exit(1); +// } +// TClonesArray &lcpvmodule = *fCPVModules; +// for (Int_t i=0; iGetNCPVModules(); i++) new(lcpvmodule[i]) AliPHOSCPVModule(); +// } +// else { +// // Create an empty array of AliPHOSCPVModule to satisfy +// // AliPHOSv1::Streamer when writing root file - fCPVModules=new TClonesArray("AliPHOSCPVModule",0); +// fCPVModules=new TClonesArray("AliPHOSCPVModule",0); - } +// } } //____________________________________________________________________________ @@ -194,6 +199,12 @@ AliPHOSv1::~AliPHOSv1() fHits = 0 ; } + if ( fSDigits) { + fSDigits->Delete() ; + delete fSDigits ; + fSDigits = 0 ; + } + if ( fDigits) { fDigits->Delete() ; delete fDigits ; @@ -221,7 +232,7 @@ AliPHOSv1::~AliPHOSv1() } //____________________________________________________________________________ -void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id, Float_t * hits, Int_t trackpid) +void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id, Float_t * hits, Int_t trackpid, TLorentzVector p, Float_t * lpos) { // Add a hit to the hit list. // A PHOS hit is the sum of all hits in a single crystal @@ -232,7 +243,7 @@ void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id, AliPHOSHit *curHit ; Bool_t deja = kFALSE ; - newHit = new AliPHOSHit(shunt, primary, tracknumber, Id, hits, trackpid) ; + newHit = new AliPHOSHit(shunt, primary, tracknumber, Id, hits, trackpid, p, lpos) ; for ( hitCounter = 0 ; hitCounter < fNhits && !deja ; hitCounter++ ) { curHit = (AliPHOSHit*) (*fHits)[hitCounter] ; @@ -250,49 +261,26 @@ void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id, delete newHit; } -//___________________________________________________________________________ -Int_t AliPHOSv1::Digitize(Float_t Energy) -{ - // Applies the energy calibration - - Float_t fB = 10000000. ; - Float_t fA = 0. ; - Int_t chan = Int_t(fA + Energy*fB ) ; - return chan ; -} - //____________________________________________________________________________ -void AliPHOSv1::Hit2Digit(Int_t ntracks){ - //Collects all hits in the same active volume into digits - - if(fDigits!= 0) - fDigits->Clear() ; - else - fDigits = new TClonesArray("AliPHOSDigit",1000) ; - - // Branch address for digit tree - char branchname[20]; - sprintf(branchname,"%s",GetName()); - gAlice->TreeD()->Branch(branchname,&fDigits,fBufferSize); - - gAlice->TreeD()->GetEvent(0); +void AliPHOSv1::Hits2SDigits(){ + //Collects all hits in the same active volume into digit - Int_t i ; - Int_t relid[4]; Int_t j ; AliPHOSHit * hit ; AliPHOSDigit * newdigit ; AliPHOSDigit * curdigit ; Bool_t deja = kFALSE ; + Int_t itrack ; - for (itrack=0; itrackGetNtrack(); itrack++){ + //=========== Get the Hits Tree for the Primary track itrack gAlice->ResetHits(); gAlice->TreeH()->GetEvent(itrack); + for ( i = 0 ; i < fHits->GetEntries() ; i++ ) { hit = (AliPHOSHit*)fHits->At(i) ; @@ -303,9 +291,8 @@ void AliPHOSv1::Hit2Digit(Int_t ntracks){ newdigit = new AliPHOSDigit( -1 , hit->GetId(), Digitize( hit->GetEnergy() ) ) ; deja =kFALSE ; - - for ( j = 0 ; j < fNdigits ; j++) { - curdigit = (AliPHOSDigit*) fDigits->At(j) ; + for ( j = 0 ; j < fnSdigits ; j++) { + curdigit = (AliPHOSDigit*) fSDigits->At(j) ; if ( *curdigit == *newdigit) { *curdigit = *curdigit + *newdigit ; deja = kTRUE ; @@ -313,44 +300,122 @@ void AliPHOSv1::Hit2Digit(Int_t ntracks){ } if ( !deja ) { - new((*fDigits)[fNdigits]) AliPHOSDigit(* newdigit) ; - fNdigits++ ; + new((*fSDigits)[fnSdigits]) AliPHOSDigit(* newdigit) ; + fnSdigits++ ; } delete newdigit ; } } // loop over tracks - - // Noise induced by the PIN diode of the PbWO crystals + + fSDigits->Sort() ; + + fnSdigits = fSDigits->GetEntries() ; + fSDigits->Expand(fnSdigits) ; + + for (i = 0 ; i < fnSdigits ; i++) { + AliPHOSDigit * digit = (AliPHOSDigit *) fSDigits->At(i) ; + digit->SetIndexInList(i) ; + } + + gAlice->TreeS()->Fill() ; - Float_t energyandnoise ; - for ( i = 0 ; i < fNdigits ; i++ ) { - newdigit = (AliPHOSDigit * ) fDigits->At(i) ; + TTree *ts = gAlice->TreeS(); + TDirectory *dir = ts->GetDirectory() ; + cout << "dir name is " << dir->GetName() << endl ; + // ts->Write(0,TObject::kOverwrite) ; + +gAlice->TreeS()->GetBranch("PHOS")->Fill(); +gAlice->TreeS()->GetBranch("PHOS")->Write(); + + +} +//____________________________________________________________________________ +void AliPHOSv1::SDigits2Digits(){ + //Adds noise to the summable digits and removes everething below thresholds + //Note, that sDigits should be SORTED in accordance with abs ID. + + + gAlice->TreeS()->GetEvent(0) ; + + cout << "fSdigits " << fSDigits << " " << fSDigits->GetEntries() << endl ; + + + // Noise induced by the PIN diode of the PbWO crystals + Int_t iCurSDigit = 0 ; + //we assume, that there is al least one EMC digit... + Int_t idCurSDigit = ((AliPHOSDigit *)fSDigits->At(0))->GetId() ; + + cout << "fDigits " << fDigits << " " <GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ(); absID++){ - fGeom->AbsToRelNumbering(newdigit->GetId(), relid) ; + cout << "absID " << absID << " " << idCurSDigit << endl ; + + Float_t noise = gRandom->Gaus(0., fPinElectronicNoise) ; + if(absID < idCurSDigit ){ + cout << "In < idC " << noise << " " << fDigitThreshold << endl ; + if(noise >fDigitThreshold ){ + cout << "noise " << absID << " " << noise << endl; + new((*fDigits)[fNdigits]) AliPHOSDigit( -1,absID,Digitize(noise) ) ; + cout << "FN " << fNdigits << endl ; + fNdigits++ ; + } + } + else{ //add noise and may be remove the true hit + cout << "correcting digit " << iCurSDigit << endl ; + Float_t signal = noise + Calibrate(((AliPHOSDigit *)fSDigits->At(iCurSDigit))->GetAmp()) ; + cout << "signal " << signal << endl ; + if( signal >fDigitThreshold ){ + cout << "signal " << signal << endl ; + AliPHOSDigit * digit = (AliPHOSDigit*) fSDigits->At(iCurSDigit) ; + new((*fDigits)[fNdigits]) AliPHOSDigit( *digit ) ; + ((AliPHOSDigit *)fDigits->At(fNdigits))->SetAmp(Digitize(signal)); + cout << "fNdigits " << fNdigits << endl ; + fNdigits++ ; + } + + if(iCurSDigit < fSDigits->GetEntries()-1){ + iCurSDigit++ ; + idCurSDigit = ((AliPHOSDigit*)fSDigits->At(iCurSDigit))->GetId() ; + } + else + idCurSDigit = 10000000; //no real hits left + } - if (relid[1]==0){ // Digits belong to EMC (PbW0_4 crystals) - energyandnoise = newdigit->GetAmp() + Digitize(gRandom->Gaus(0., fPinElectronicNoise)) ; - - if (energyandnoise < 0 ) - energyandnoise = 0 ; - - if ( newdigit->GetAmp() < fDigitThreshold ) // if threshold not surpassed, remove digit from list - fDigits->RemoveAt(i) ; + } + + //remove PPSD/CPV digits below thresholds + Int_t idigit ; + for(idigit = iCurSDigit; idigit < fSDigits->GetEntries() ; idigit++){ //loop over CPV/PPSD digits + + AliPHOSDigit * digit = (AliPHOSDigit *) fSDigits->At(idigit) ; + Float_t ene = Calibrate(digit->GetAmp()) ; + + Int_t relid[4] ; + fGeom->AbsToRelNumbering(digit->GetId(), relid) ; + if ( relid[0] > fGeom->GetNCPVModules() ){ //ppsd + if ( ( (relid[1] > 0) && (ene > fPpsdEnergyThreshold)) || //PPSD digit + ( (relid[1] < 0) && (ene > fCpvEnergyThreshold ) ) ) //CPV digit + new((*fDigits)[fNdigits]) AliPHOSDigit( *digit ) ; + fNdigits++ ; } - } + } fDigits->Compress() ; fNdigits = fDigits->GetEntries() ; fDigits->Expand(fNdigits) ; + Int_t i ; for (i = 0 ; i < fNdigits ; i++) { - newdigit = (AliPHOSDigit *) fDigits->At(i) ; - newdigit->SetIndexInList(i) ; + AliPHOSDigit * digit = (AliPHOSDigit *) fDigits->At(i) ; + digit->SetIndexInList(i) ; } + gAlice->TreeD()->Fill() ; gAlice->TreeD()->Write(0,TObject::kOverwrite) ; @@ -359,20 +424,81 @@ void AliPHOSv1::Hit2Digit(Int_t ntracks){ //___________________________________________________________________________ void AliPHOSv1::MakeBranch(Option_t* opt, char *file) -{ +{ + + + char *cH ; // Create new branche in the current Root Tree in the digit Tree AliDetector::MakeBranch(opt) ; - - // Create new branches EMC for hits in EMC modules - - for( Int_t i=0; iGetNModules(); i++ ) GetEMCModule(i).MakeBranch("EMC",i+1,file); - - // Create new branches CPV for hits in CPV modules for IHEP geometry - - if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) { - for( Int_t i=0; iGetNCPVModules(); i++ ) GetCPVModule(i).MakeBranch("CPV",i+1,file); + + + cH = strstr(opt,"S"); + //Create a branch for SDigits + if( cH ){ + char branchname[20]; + sprintf(branchname,"%s",GetName()); + if(fSDigits) + fSDigits->Clear(); + else + fSDigits = new TClonesArray("AliPHOSDigit",1000); + fnSdigits = 0 ; + cout << " AliPHOSv1::MakeBranch : " << file << endl ; + gAlice->MakeBranchInTree(gAlice->TreeS(),branchname,&fSDigits,fBufferSize,file); } + cH = strstr(opt,"D"); + //Create a branch for Digits + if( cH ){ + char branchname[20]; + sprintf(branchname,"%s",GetName()); + if(fSDigits) + fDigits->Clear(); + else + fDigits = new TClonesArray("AliPHOSDigit",1000); + fNdigits = 0 ; + + gAlice->MakeBranchInTree(gAlice->TreeD(),branchname,&fSDigits,fBufferSize,file); + } + + cH = strstr(opt,"R"); + //Create a branch for Reconstruction + if( cH ){ + char branchname[20]; + + Int_t splitlevel = 0 ; + + fEmcRecPoints->Delete() ; + + if ( fEmcRecPoints && gAlice->TreeR() ) { + sprintf(branchname,"%sEmcRP",GetName()) ; + gAlice->TreeR()->Branch(branchname, "TObjArray", &fEmcRecPoints, fBufferSize, splitlevel) ; + } + + fPpsdRecPoints->Delete() ; + + if ( fPpsdRecPoints && gAlice->TreeR() ) { + sprintf(branchname,"%sPpsdRP",GetName()) ; + gAlice->TreeR()->Branch(branchname, "TObjArray", &fPpsdRecPoints, fBufferSize, splitlevel) ; + } + + fTrackSegments->Delete() ; + + if ( fTrackSegments && gAlice->TreeR() ) { + sprintf(branchname,"%sTS",GetName()) ; + gAlice->TreeR()->Branch(branchname, &fTrackSegments, fBufferSize) ; + } + + fRecParticles->Delete() ; + + if ( fRecParticles && gAlice->TreeR() ) { + sprintf(branchname,"%sRP",GetName()) ; + gAlice->TreeR()->Branch(branchname, &fRecParticles, fBufferSize) ; + } + + } + + + } //_____________________________________________________________________________ @@ -385,41 +511,12 @@ void AliPHOSv1::Reconstruction(AliPHOSReconstructioner * Reconstructioner) // 5. Write the Tree to File fReconstructioner = Reconstructioner ; - - char branchname[10] ; - + // 1. // gAlice->MakeTree("R") ; - Int_t splitlevel = 0 ; - fEmcRecPoints->Delete() ; - - if ( fEmcRecPoints && gAlice->TreeR() ) { - sprintf(branchname,"%sEmcRP",GetName()) ; - gAlice->TreeR()->Branch(branchname, "TObjArray", &fEmcRecPoints, fBufferSize, splitlevel) ; - } - - fPpsdRecPoints->Delete() ; - - if ( fPpsdRecPoints && gAlice->TreeR() ) { - sprintf(branchname,"%sPpsdRP",GetName()) ; - gAlice->TreeR()->Branch(branchname, "TObjArray", &fPpsdRecPoints, fBufferSize, splitlevel) ; - } - - fTrackSegments->Delete() ; - - if ( fTrackSegments && gAlice->TreeR() ) { - sprintf(branchname,"%sTS",GetName()) ; - gAlice->TreeR()->Branch(branchname, &fTrackSegments, fBufferSize) ; - } - - fRecParticles->Delete() ; - - if ( fRecParticles && gAlice->TreeR() ) { - sprintf(branchname,"%sRP",GetName()) ; - gAlice->TreeR()->Branch(branchname, &fRecParticles, fBufferSize) ; - } + MakeBranch("R") ; // 3. @@ -462,10 +559,10 @@ void AliPHOSv1::ResetHits() // Yuri Kharlov, 28 September 2000 AliDetector::ResetHits(); - for (Int_t i=0; iGetNModules(); i++) ((AliPHOSCPVModule*)(*fEMCModules)[i]) -> Clear(); - if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) { - for (Int_t i=0; iGetNCPVModules(); i++) ((AliPHOSCPVModule*)(*fCPVModules)[i]) -> Clear(); - } + // for (Int_t i=0; iGetNModules(); i++) ((AliPHOSCPVModule*)(*fEMCModules)[i]) -> Clear(); + // if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) { + // for (Int_t i=0; iGetNCPVModules(); i++) ((AliPHOSCPVModule*)(*fCPVModules)[i]) -> Clear(); + // } } //____________________________________________________________________________ @@ -480,6 +577,12 @@ void AliPHOSv1::ResetReconstruction() } +//____________________________________________________________________________ +//void AliPHOSv1::SDigits2Digits() { +// // Adds the noise to the summable digits and keeps digits above a threshold +// // To make a digit +//} + //____________________________________________________________________________ void AliPHOSv1::SetTreeAddress() { @@ -496,14 +599,14 @@ void AliPHOSv1::SetTreeAddress() // Set branch address for the Hits Tree for hits in EMC modules // Yuri Kharlov, 23 November 2000. - for( Int_t i=0; iGetNModules(); i++ ) GetEMCModule(i).SetTreeAddress("EMC",i+1); + // for( Int_t i=0; iGetNModules(); i++ ) GetEMCModule(i).SetTreeAddress("EMC",i+1); // Set branch address for the Hits Tree for hits in CPV modules for IHEP geometry // Yuri Kharlov, 28 September 2000. - if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) { - for( Int_t i=0; iGetNCPVModules(); i++ ) GetCPVModule(i).SetTreeAddress("CPV",i+1); - } + // if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) { + // for( Int_t i=0; iGetNCPVModules(); i++ ) GetCPVModule(i).SetTreeAddress("CPV",i+1); + // } } @@ -522,13 +625,43 @@ void AliPHOSv1::StepManager(void) Int_t absid ; // absolute cell ID number Float_t xyze[4]={0,0,0,0} ; // position wrt MRS and energy deposited TLorentzVector pos ; // Lorentz vector of the track current position + TLorentzVector pmom ; //momentum of the particle initiated hit + Float_t xyd[2] ; //local posiiton of the entering + Bool_t entered = kFALSE ; Int_t copy ; - Int_t i ; Int_t tracknumber = gAlice->CurrentTrack() ; Int_t primary = gAlice->GetPrimary( gAlice->CurrentTrack() ); TString name = fGeom->GetName() ; - Int_t trackpid = gMC->TrackPid() ; + Int_t trackpid = 0 ; + + if( gMC->IsTrackEntering() ){ // create hit with position and momentum of new particle, + // but may be without energy deposition + + // Current position of the hit in the local ref. system + 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 + xyd[0] = xyzd[0]; + xyd[1] =-xyzd[2]; + + // Current momentum of the hit's track in the local ref. system + 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]; + + trackpid = gMC->TrackPid(); + entered = kTRUE ; // Mark to create hit even withou energy deposition + + } + + if ( name == "GPS2" || name == "MIXT" ) { // ======> CPV is a GPS' PPSD if( gMC->CurrentVolID(copy) == gMC->VolId("GCEL") ) // We are inside a gas cell @@ -539,7 +672,7 @@ void AliPHOSv1::StepManager(void) xyze[2] = pos[2] ; xyze[3] = gMC->Edep() ; - if ( xyze[3] != 0 ) { // there is deposited energy + if ( (xyze[3] != 0) || entered ) { // there is deposited energy or new particle entering PPSD gMC->CurrentVolOffID(5, relid[0]) ; // get the PHOS Module number if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(5),"PHO1") == 0 ){ relid[0] += fGeom->GetNModules() - fGeom->GetNPPSDModules(); @@ -555,7 +688,8 @@ void AliPHOSv1::StepManager(void) fGeom->RelToAbsNumbering(relid, absid) ; // add current hit to the hit list - AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid); + AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid, pmom, xyd); + } // there is deposited energy } // We are inside the gas of the CPV @@ -566,57 +700,18 @@ void AliPHOSv1::StepManager(void) // Yuri Kharlov, 28 September 2000 if( gMC->CurrentVolID(copy) == gMC->VolId("CPVQ") && - gMC->IsTrackEntering() && - gMC->TrackCharge() != 0) { - - // Charged track has just entered to the CPV sensitive plane + entered && + gMC->TrackCharge() != 0) { - AliPHOSv1 &phos = *(AliPHOSv1*)gAlice->GetModule("PHOS"); + // Digitize the current CPV hit: + + // 1. find pad response and Int_t moduleNumber; gMC->CurrentVolOffID(3,moduleNumber); moduleNumber--; - - // Current position of the hit in the CPV module ref. system - - gMC -> TrackPosition(pos); - Float_t xyzm[3], xyzd[3], xyd[3]={0,0,0}; - 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 - xyd[0] = xyzd[0]; - xyd[1] =-xyzd[2]; - - // Current momentum of the hit's track in the CPV module ref. system - - TLorentzVector pmom; - gMC -> TrackMomentum(pmom); - Float_t pm[3], pd[3]; - for (i=0; i<3; i++) pm[i] = pmom[i]; - gMC -> Gmtod (pm, pd, 2); // transform 3-momentum from master to daughter system - pmom[0] = pd[0]; - pmom[1] =-pd[1]; - pmom[2] =-pd[2]; - - // Current particle type of the hit's track - Int_t ipart = gMC->TrackPid(); - - // Add the current particle in the list of the CPV hits. - - phos.GetCPVModule(moduleNumber).AddHit(fIshunt,primary,pmom,xyd,ipart); - - if (fDebug == 1) { - printf("CPV hit added to module #%2d: p = (% .4f, % .4f, % .4f, % .4f) GeV,\n", - moduleNumber+1,pmom.Px(),pmom.Py(),pmom.Pz(),pmom.E()); - printf( " xy = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n", - xyd[0],xyd[1],ipart,primary); - } - - // Digitize the current CPV hit: - // 1. find pad response and - TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0); // array of digits for current hit CPVDigitize(pmom,xyd,moduleNumber,cpvDigits); @@ -664,7 +759,7 @@ void AliPHOSv1::StepManager(void) 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, trackpid); + AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid, pmom, xyd); if (cpvDigit->GetQpad() > 0.02) { xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5); @@ -676,6 +771,7 @@ void AliPHOSv1::StepManager(void) } } // end of IHEP configuration + if(gMC->CurrentVolID(copy) == gMC->VolId("PXTL") ) { // We are inside a PBWO crystal gMC->TrackPosition(pos) ; xyze[0] = pos[0] ; @@ -683,62 +779,11 @@ void AliPHOSv1::StepManager(void) xyze[2] = pos[2] ; xyze[3] = gMC->Edep() ; - // Track enters to the crystal from the top edge - - if (gMC->IsTrackEntering()) { - Float_t posloc[3]; - gMC -> Gmtod (xyze, posloc, 1); - if (posloc[1] > fGeom->GetCrystalSize(1)/2-0.01) { - Int_t row,cel; - Float_t xyd[3]={0,0,0}; - AliPHOSv1 &phos = *(AliPHOSv1*)gAlice->GetModule("PHOS"); - - Int_t moduleNumber; - gMC->CurrentVolOffID(10,moduleNumber); - - if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(10),"PHO1") == 0 ) - moduleNumber += fGeom->GetNModules() - fGeom->GetNPPSDModules(); - moduleNumber--; - - gMC->CurrentVolOffID(4, row) ; // get the row number inside the module - gMC->CurrentVolOffID(3, cel) ; // get the cell number inside the module - xyd[0] = -(posloc[2] + (cel-0.5-fGeom->GetNZ() /2) * - (fGeom->GetCrystalSize(2) + 2 * fGeom->GetGapBetweenCrystals())); - xyd[1] = posloc[0] + (row-0.5-fGeom->GetNPhi()/2) * - (fGeom->GetCrystalSize(0) + 2 * fGeom->GetGapBetweenCrystals()); - - // Current momentum of the hit's track in the CPV module ref. system - - TLorentzVector pmom; - gMC -> TrackMomentum(pmom); - Float_t pm[3], pd[3]; - for (i=0; i<3; i++) pm[i] = pmom[i]; - gMC -> Gmtod (pm, pd, 2); // transform 3-momentum from master to daughter system - pmom[0] = pd[0]; - pmom[1] =-pd[1]; - pmom[2] =-pd[2]; - - // Current particle type of the hit's track - - Int_t ipart = gMC->TrackPid(); - - // Add the current particle in the list of the EMC hits. - - phos.GetEMCModule(moduleNumber).AddHit(fIshunt,primary,pmom,xyd,ipart); - - if (fDebug == 1) { - printf("EMC hit added to module #%2d: p = (% .4f, % .4f, % .4f, % .4f) GeV,\n", - moduleNumber+1,pmom.Px(),pmom.Py(),pmom.Pz(),pmom.E()); - printf( " xy = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n", - xyd[0],xyd[1],ipart,primary); - } - } - } - - // Track is inside the crystal and deposits some energy + + if ( (xyze[3] != 0) || entered ) { // Track is inside the crystal and deposits some energy or just entered - if ( xyze[3] != 0 ) { gMC->CurrentVolOffID(10, relid[0]) ; // get the PHOS module number ; + if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(10),"PHO1") == 0 ) relid[0] += fGeom->GetNModules() - fGeom->GetNPPSDModules(); @@ -750,10 +795,13 @@ void AliPHOSv1::StepManager(void) fGeom->RelToAbsNumbering(relid, absid) ; // add current hit to the hit list - AddHit(fIshunt, primary,tracknumber, absid, xyze, trackpid); + AddHit(fIshunt, primary,tracknumber, absid, xyze, trackpid,pmom, xyd); + } // there is deposited energy } // we are inside a PHOS Xtal + + } //____________________________________________________________________________ diff --git a/PHOS/AliPHOSv1.h b/PHOS/AliPHOSv1.h index c5191285496..741ae58a1e1 100644 --- a/PHOS/AliPHOSv1.h +++ b/PHOS/AliPHOSv1.h @@ -39,23 +39,30 @@ public: } virtual ~AliPHOSv1(void) ; - virtual void AddHit( Int_t shunt, Int_t primary, Int_t track, Int_t id, Float_t *hits, Int_t pid) ; - Int_t Digitize(Float_t Energy); - virtual void Hit2Digit(Int_t event) ; - virtual Int_t IsVersion(void) const { - // Gives the version number - return 1 ; - } + virtual void AddHit( Int_t shunt, Int_t primary, Int_t track, Int_t id, Float_t *hits, Int_t pid, TLorentzVector p, Float_t *pos) ; + Float_t Calibrate(Int_t amp){ return (amp - fDigitizeA)/fDigitizeB ; } + Int_t Digitize(Float_t Energy){ return (Int_t ) (fDigitizeA + Energy*fDigitizeB); } + // virtual void Hit2Digit(Int_t event) ; + virtual void Hits2SDigits() ; virtual void MakeBranch(Option_t* opt, char *file=0 ) ; void Reconstruction(AliPHOSReconstructioner * Reconstructioner) ; void ResetClusters(){} ; virtual void ResetHits() ; + virtual void SDigits2Digits() ; + virtual Int_t IsVersion(void) const { + // Gives the version number + return 1 ; + } + virtual void ResetReconstruction() ; // Reset reconstructed objects void SetReconstructioner(AliPHOSReconstructioner& Reconstructioner) { // sets the reconstructionner object to be used fReconstructioner = &Reconstructioner ; } void SetDigitThreshold(Float_t th) { fDigitThreshold = th ; } + void SetPpsdEnergyThreshold(Float_t enth) { fPpsdEnergyThreshold = enth ; } + void SetCpvEnergyThreshold(Float_t enth) { fCpvEnergyThreshold = enth ; } + virtual void SetTreeAddress(); virtual void StepManager(void) ; virtual TString Version(void){ @@ -71,8 +78,8 @@ public: // IHEP's CPV specific functions - AliPHOSCPVModule &GetEMCModule(int n) { return *(AliPHOSCPVModule*)fEMCModules->operator[](n); } - AliPHOSCPVModule &GetCPVModule(int n) { return *(AliPHOSCPVModule*)fCPVModules->operator[](n); } + // AliPHOSCPVModule &GetEMCModule(int n) { return *(AliPHOSCPVModule*)fEMCModules->operator[](n); } + // AliPHOSCPVModule &GetCPVModule(int n) { return *(AliPHOSCPVModule*)fCPVModules->operator[](n); } void CPVDigitize (TLorentzVector p, Float_t *xy, Int_t moduleNumber, TClonesArray *digits) ; Float_t CPVPadResponseFunction(Float_t qhit, Float_t zg, Float_t xg) ; @@ -81,11 +88,16 @@ public: protected: Float_t fDigitThreshold ; // Threshold for the digit registration + Float_t fPpsdEnergyThreshold; //PPSD + Float_t fCpvEnergyThreshold; //CPV Float_t fPinElectronicNoise ; // Electronic Noise in the PIN - AliPHOSReconstructioner * fReconstructioner ; // Reconstrutioner of the PHOS event: Clusterization and subtracking procedures + Float_t fDigitizeA ; //Parameters of the + Float_t fDigitizeB ; //digitization + Int_t fnSdigits ; + AliPHOSReconstructioner * fReconstructioner ; // Clusterization and subtracking procedures AliPHOSTrackSegmentMaker * fTrackSegmentMaker ; // Reconstructioner of the PHOS track segment: 2 x PPSD + 1 x EMC - TClonesArray * fEMCModules; // Array of EMC modules - TClonesArray * fCPVModules; // Array of CPV modules for the IHEP's version of CPV + // TClonesArray * fEMCModules; // Array of EMC modules + // TClonesArray * fCPVModules; // Array of CPV modules for the IHEP's version of CPV ClassDef(AliPHOSv1,1) // Implementation of PHOS manager class for layout EMC+PPSD diff --git a/PHOS/AliPHOSv2.cxx b/PHOS/AliPHOSv2.cxx index 6de2a51fb44..838143c38be 100644 --- a/PHOS/AliPHOSv2.cxx +++ b/PHOS/AliPHOSv2.cxx @@ -69,13 +69,13 @@ AliPHOSv2::~AliPHOSv2() } //____________________________________________________________________________ -void AliPHOSv2::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id, Float_t * hits, Int_t pid) +void AliPHOSv2::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id, Float_t * hits, Int_t pid, TLorentzVector p, Float_t * lpos) { // Add a hit to the hit list. AliPHOSHit *newHit ; - newHit = new AliPHOSHit(shunt, primary, tracknumber, Id, hits, pid) ; + newHit = new AliPHOSHit(shunt, primary, tracknumber, Id, hits, pid, p, lpos) ; new((*fHits)[fNhits]) AliPHOSHit(*newHit) ; fNhits++ ; diff --git a/PHOS/AliPHOSv2.h b/PHOS/AliPHOSv2.h index 340887600f5..6e6e01ae260 100644 --- a/PHOS/AliPHOSv2.h +++ b/PHOS/AliPHOSv2.h @@ -24,7 +24,7 @@ public: AliPHOSv2(const char *name, const char *title="") ; virtual ~AliPHOSv2(void) ; - virtual void AddHit( Int_t shunt, Int_t primary, Int_t track, Int_t id, Float_t *hits, Int_t pid ) ; + virtual void AddHit( Int_t shunt, Int_t primary, Int_t track, Int_t id, Float_t *hits, Int_t pid, TLorentzVector p, Float_t *pos ) ; virtual Int_t IsVersion(void) const { // Gives the version number return 2 ; diff --git a/PHOS/AliPHOSv3.cxx b/PHOS/AliPHOSv3.cxx index 49877304ede..87a718c13fc 100644 --- a/PHOS/AliPHOSv3.cxx +++ b/PHOS/AliPHOSv3.cxx @@ -92,123 +92,225 @@ AliPHOSv3::AliPHOSv3(AliPHOSReconstructioner * Reconstructioner, const char *nam fRecalibrationFactor = 6.2 / fLightYieldMean ; fElectronsPerGeV = 2.77e+8 ; } - //____________________________________________________________________________ + void AliPHOSv3::StepManager(void) { // Accumulates hits as long as the track stays in a single crystal or PPSD gas Cell - // Adds the energy deposited in the PIN diode + +// if (gMC->IsTrackEntering()) +// cout << "Track enters the volume " << gMC->CurrentVolName() << endl; +// if (gMC->IsTrackExiting()) +// cout << "Track leaves the volume " << gMC->CurrentVolName() << endl; Int_t relid[4] ; // (box, layer, row, column) indices - Float_t xyze[4] ; // position wrt MRS and energy deposited - TLorentzVector pos ; - Int_t copy; - Float_t lightyield ; // Light Yield per GeV - Float_t nElectrons ; // Number of electrons in the PIN diode - TString name = fGeom->GetName() ; - Float_t global[3] ; - Float_t local[3] ; - Float_t lostenergy ; + Int_t absid ; // absolute cell ID number + Float_t xyze[4]={0,0,0,0} ; // position wrt MRS and energy deposited + TLorentzVector pos ; // Lorentz vector of the track current position + TLorentzVector pmom ; //momentum of the particle initiated hit + Float_t xyd[2] ; //local posiiton of the entering + Bool_t entered = kFALSE ; + Int_t copy ; Int_t tracknumber = gAlice->CurrentTrack() ; - Int_t primary = gAlice->GetPrimary( gAlice->CurrentTrack() ); - Int_t trackpid = gMC->TrackPid() ; + Int_t primary = gAlice->GetPrimary( gAlice->CurrentTrack() ); + TString name = fGeom->GetName() ; + Int_t trackpid = 0 ; + + if( gMC->IsTrackEntering() ){ // create hit with position and momentum of new particle, + // but may be without energy deposition + + // Current position of the hit in the local ref. system + 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 + xyd[0] = xyzd[0]; + xyd[1] =-xyzd[2]; + + // Current momentum of the hit's track in the local ref. system + 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]; + + trackpid = gMC->TrackPid(); + entered = kTRUE ; // Mark to create hit even withou energy deposition + + } + + + if ( name == "GPS2" || name == "MIXT" ) { // ======> CPV is a GPS' PPSD - if ( name == "GPS2" ) { // the CPV is a PPSD if( gMC->CurrentVolID(copy) == gMC->VolId("GCEL") ) // 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() ; - + xyze[3] = gMC->Edep() ; - if ( xyze[3] != 0 ) { // there is deposited energy + if ( (xyze[3] != 0) || entered ) { // there is deposited energy or new particle entering PPSD gMC->CurrentVolOffID(5, relid[0]) ; // get the PHOS Module number + if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(5),"PHO1") == 0 ){ + relid[0] += fGeom->GetNModules() - fGeom->GetNPPSDModules(); + } gMC->CurrentVolOffID(3, relid[1]) ; // get the Micromegas Module number - // 1-> Geom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() upper - // > fGeom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() lower + // 1-> fGeom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() upper + // > fGeom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() lower gMC->CurrentVolOffID(1, relid[2]) ; // get the row number of the cell gMC->CurrentVolID(relid[3]) ; // get the column number // get the absolute Id number - Int_t absid ; - fGeom->RelToAbsNumbering(relid,absid) ; - + fGeom->RelToAbsNumbering(relid, absid) ; + + // add current hit to the hit list + AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid, pmom, xyd); - AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid); } // there is deposited energy - } // We are inside the gas of the CPV - } // GPS2 configuration + } // 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("CPVQ") && + entered && + gMC->TrackCharge() != 0) { + + // Digitize the current CPV hit: + + // 1. find pad response and + + 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); + + 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 + + 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 + fGeom->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, trackpid, pmom, xyd); + + 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 - 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] ; - lostenergy = gMC->Edep() ; - xyze[3] = gMC->Edep() ; - - global[0] = pos[0] ; - global[1] = pos[1] ; - global[2] = pos[2] ; - - if ( xyze[3] != 0 ) { - gMC->CurrentVolOffID(10, relid[0]) ; // get the PHOS module number ; - relid[1] = 0 ; // means PW04 - 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 + 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() ; - Int_t absid ; - fGeom->RelToAbsNumbering(relid,absid) ; - gMC->Gmtod(global, local, 1) ; - - // calculating number of electrons in the PIN diode asociated to this hit - lightyield = gRandom->Poisson(fLightYieldMean) ; - nElectrons = lostenergy * lightyield * fIntrinsicPINEfficiency * - exp(-fLightYieldAttenuation * (local[1]+fGeom->GetCrystalSize(1)/2.0 ) ) ; + + if ( (xyze[3] != 0) || entered ) { // Track is inside the crystal and deposits some energy or just entered - xyze[3] = nElectrons * fRecalibrationFactor ; - // add current hit to the hit list - AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid); - - } // there is deposited energy - } // we are inside a PHOS Xtal - - if(gMC->CurrentVolID(copy) == gMC->VolId("PPIN") ) // We are inside de PIN diode - { - gMC->TrackPosition(pos) ; - xyze[0] = pos[0] ; - xyze[1] = pos[1] ; - xyze[2] = pos[2] ; - lostenergy = gMC->Edep() ; - xyze[3] = gMC->Edep() ; - - if ( xyze[3] != 0 ) { - gMC->CurrentVolOffID(11, relid[0]) ; // get the PHOS module number ; - relid[1] = 0 ; // means PW04 and PIN - gMC->CurrentVolOffID(5, relid[2]) ; // get the row number inside the module - gMC->CurrentVolOffID(4, relid[3]) ; // get the cell number inside the module + gMC->CurrentVolOffID(10, relid[0]) ; // get the PHOS module number ; + if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(10),"PHO1") == 0 ) + relid[0] += fGeom->GetNModules() - fGeom->GetNPPSDModules(); + + relid[1] = 0 ; // means PBW04 + gMC->CurrentVolOffID(4, relid[2]) ; // get the row number inside the module + gMC->CurrentVolOffID(3, relid[3]) ; // get the cell number inside the module + // get the absolute Id number + fGeom->RelToAbsNumbering(relid, absid) ; + + // add current hit to the hit list + AddHit(fIshunt, primary,tracknumber, absid, xyze, trackpid,pmom, xyd); - Int_t absid ; - fGeom->RelToAbsNumbering(relid,absid) ; - - // calculating number of electrons in the PIN diode asociated to this hit - nElectrons = lostenergy * fElectronsPerGeV ; - xyze[3] = nElectrons * fRecalibrationFactor ; + } // there is deposited energy + } // we are inside a PHOS Xtal + + if(gMC->CurrentVolID(copy) == gMC->VolId("PPIN") ) // We are inside de PIN diode + { + gMC->TrackPosition(pos) ; + xyze[0] = pos[0] ; + xyze[1] = pos[1] ; + xyze[2] = pos[2] ; + Float_t lostenergy = gMC->Edep() ; + xyze[3] = gMC->Edep() ; + + if ( xyze[3] != 0 ) { + gMC->CurrentVolOffID(11, relid[0]) ; // get the PHOS module number ; + relid[1] = 0 ; // means PW04 and PIN + gMC->CurrentVolOffID(5, relid[2]) ; // get the row number inside the module + gMC->CurrentVolOffID(4, relid[3]) ; // get the cell number inside the module + + // get the absolute Id number + + Int_t absid ; + fGeom->RelToAbsNumbering(relid,absid) ; + + // calculating number of electrons in the PIN diode asociated to this hit + Float_t nElectrons = lostenergy * fElectronsPerGeV ; + xyze[3] = nElectrons * fRecalibrationFactor ; + // add current hit to the hit list - AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid); + AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid,pmom,xyd); //printf("PIN volume is %d, %d, %d, %d \n",relid[0],relid[1],relid[2],relid[3]); //printf("Lost energy in the PIN is %f \n",lostenergy) ; - } // there is deposited energy + } // there is deposited energy } // we are inside a PHOS XtalPHOS PIN diode } -- 2.31.1