// --- Standard library ---
+#include <strstream.h>
// --- AliRoot header files ---
#include "AliPHOS.h"
//____________________________________________________________________________
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
}
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
(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);
//____________________________________________________________________________
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<Nevents; ievent++) {
+// //
+// // 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<Nevents; ievent++) {
- //========== 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 ;
+// //========== 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; itrack<ntracks; itrack++) {
- // Get the Hits Tree for the Primary track itrack
- gAlice->ResetHits();
- gAlice->TreeH()->GetEvent(itrack);
- for (Int_t iModule=0; iModule < nOfModules; iModule++) {
- cpvModule = fPHOS->GetCPVModule(iModule);
- cpvHits = cpvModule.Hits();
- nCPVhits = cpvHits->GetEntriesFast();
- for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
- cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(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; ihit<nCPVhits; ihit++) {
- cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(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; itrack<ntracks; itrack++) {
+// // Get the Hits Tree for the Primary track itrack
+// gAlice->ResetHits();
+// gAlice->TreeH()->GetEvent(itrack);
+// for (Int_t iModule=0; iModule < nOfModules; iModule++) {
+// cpvModule = fPHOS->GetCPVModule(iModule);
+// cpvHits = cpvModule.Hits();
+// nCPVhits = cpvHits->GetEntriesFast();
+// for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
+// cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(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; ihit<nCPVhits; ihit++) {
+// cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(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");
}
//____________________________________________________________________________
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; itrack<ntracks; itrack++) {
- //=========== Get the Hits Tree for the Primary track itrack
- gAlice->ResetHits();
- gAlice->TreeH()->GetEvent(itrack);
- Int_t iModule = 0 ;
- for (iModule=0; iModule < fGeom->GetNModules(); iModule++) {
- emcModule = fPHOS->GetEMCModule(iModule);
- emcHits = emcModule.Hits();
- nEMChits = emcHits->GetEntriesFast();
- for (Int_t ihit=0; ihit<nEMChits; ihit++) {
- nGenHits++;
- emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
- p = emcHit->GetMomentum();
- xgen = emcHit->X();
- zgen = emcHit->Y();
- ipart = emcHit->GetIpart();
- primary= emcHit->GetTrack();
- printf("EMC hit A: module %d, ",iModule+1);
- printf(" p = (%f .4, %f .4, %f .4, %f .4) GeV,\n",
- p.Px(),p.Py(),p.Pz(),p.Energy());
- printf(" (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
- xgen,zgen,ipart,primary);
- }
- }
- }
+// AliPHOSCPVModule emcModule;
+// TClonesArray *emcHits;
+// Int_t nEMChits;
+// AliPHOSCPVHit *emcHit;
+// TLorentzVector p;
+// Float_t xgen, zgen;
+// Int_t ipart, primary;
+// Int_t nGenHits = 0;
+// for (Int_t itrack=0; itrack<ntracks; itrack++) {
+// //=========== Get the Hits Tree for the Primary track itrack
+// gAlice->ResetHits();
+// gAlice->TreeH()->GetEvent(itrack);
+// Int_t iModule = 0 ;
+// for (iModule=0; iModule < fGeom->GetNModules(); iModule++) {
+// emcModule = fPHOS->GetEMCModule(iModule);
+// emcHits = emcModule.Hits();
+// nEMChits = emcHits->GetEntriesFast();
+// for (Int_t ihit=0; ihit<nEMChits; ihit++) {
+// nGenHits++;
+// emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
+// p = emcHit->GetMomentum();
+// xgen = emcHit->X();
+// zgen = emcHit->Y();
+// ipart = emcHit->GetIpart();
+// primary= emcHit->GetTrack();
+// printf("EMC hit A: module %d, ",iModule+1);
+// printf(" p = (%f .4, %f .4, %f .4, %f .4) GeV,\n",
+// p.Px(),p.Py(),p.Pz(),p.Energy());
+// printf(" (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
+// xgen,zgen,ipart,primary);
+// }
+// }
+// }
-// // Read and print EMC hits from PHOS branch
-
-// for (Int_t itrack=0; itrack<ntracks; itrack++) {
-// //=========== Get the Hits Tree for the Primary track itrack
-// gAlice->ResetHits();
-// gAlice->TreeH()->GetEvent(itrack);
-// TClonesArray *hits = fPHOS->Hits();
-// AliPHOSHit *hit ;
-// Int_t ihit;
-// for ( ihit = 0 ; ihit < hits->GetEntries() ; ihit++ ) {
-// hit = (AliPHOSHit*)hits->At(ihit) ;
-// Float_t hitXYZ[3];
-// hitXYZ[0] = hit->X();
-// hitXYZ[1] = hit->Y();
-// hitXYZ[2] = hit->Z();
-// ipart = hit->GetPid();
-// primary = hit->GetPrimary();
-// Int_t absId = hit->GetId();
-// Int_t relId[4];
-// fGeom->AbsToRelNumbering(absId, relId) ;
-// Int_t module = relId[0];
-// if (relId[1]==0 && !(hitXYZ[0]==0 && hitXYZ[2]==0))
-// printf("EMC hit B: module %d, (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
-// module,hitXYZ[0],hitXYZ[2],ipart,primary);
-// }
-// }
+// // // Read and print EMC hits from PHOS branch
+
+// // for (Int_t itrack=0; itrack<ntracks; itrack++) {
+// // //=========== Get the Hits Tree for the Primary track itrack
+// // gAlice->ResetHits();
+// // gAlice->TreeH()->GetEvent(itrack);
+// // TClonesArray *hits = fPHOS->Hits();
+// // AliPHOSHit *hit ;
+// // Int_t ihit;
+// // for ( ihit = 0 ; ihit < hits->GetEntries() ; ihit++ ) {
+// // hit = (AliPHOSHit*)hits->At(ihit) ;
+// // Float_t hitXYZ[3];
+// // hitXYZ[0] = hit->X();
+// // hitXYZ[1] = hit->Y();
+// // hitXYZ[2] = hit->Z();
+// // ipart = hit->GetPid();
+// // primary = hit->GetPrimary();
+// // Int_t absId = hit->GetId();
+// // Int_t relId[4];
+// // fGeom->AbsToRelNumbering(absId, relId) ;
+// // Int_t module = relId[0];
+// // if (relId[1]==0 && !(hitXYZ[0]==0 && hitXYZ[2]==0))
+// // printf("EMC hit B: module %d, (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
+// // module,hitXYZ[0],hitXYZ[2],ipart,primary);
+// // }
+// // }
- }
+// }
}
//____________________________________________________________________________
void AliPHOSAnalyze::AnalyzeEMC(Int_t Nevents)
{
- //
- // Read generated and reconstructed hits in EMC for Nevents events.
- // Plots the coordinate and energy resolution histograms.
- // Coordinate resolution is a difference between the reconstructed
- // coordinate and the exact coordinate on the face of the PHOS
- // Author: Yuri Kharlov
- // 27 November 2000
- //
-
- // Book histograms
-
- TH1F *hDx1 = new TH1F("hDx1" ,"EMC x-resolution", 100,-5. , 5.);
- TH1F *hDz1 = new TH1F("hDz1" ,"EMC z-resolution", 100,-5. , 5.);
- TH1F *hDE1 = new TH1F("hDE1" ,"EMC E-resolution", 100,-2. , 2.);
-
- TH2F *hDx2 = new TH2F("hDx2" ,"EMC x-resolution", 100, 0., 10., 100,-5. , 5.);
- TH2F *hDz2 = new TH2F("hDz2" ,"EMC z-resolution", 100, 0., 10., 100,-5. , 5.);
- TH2F *hDE2 = new TH2F("hDE2" ,"EMC E-resolution", 100, 0., 10., 100, 0. , 5.);
-
- cout << "Start EMC Analysis"<< endl ;
- for (Int_t ievent=0; ievent<Nevents; ievent++) {
+// //
+// // Read generated and reconstructed hits in EMC for Nevents events.
+// // Plots the coordinate and energy resolution histograms.
+// // Coordinate resolution is a difference between the reconstructed
+// // coordinate and the exact coordinate on the face of the PHOS
+// // Author: Yuri Kharlov
+// // 27 November 2000
+// //
+
+// // Book histograms
+
+// TH1F *hDx1 = new TH1F("hDx1" ,"EMC x-resolution", 100,-5. , 5.);
+// TH1F *hDz1 = new TH1F("hDz1" ,"EMC z-resolution", 100,-5. , 5.);
+// TH1F *hDE1 = new TH1F("hDE1" ,"EMC E-resolution", 100,-2. , 2.);
+
+// TH2F *hDx2 = new TH2F("hDx2" ,"EMC x-resolution", 100, 0., 10., 100,-5. , 5.);
+// TH2F *hDz2 = new TH2F("hDz2" ,"EMC z-resolution", 100, 0., 10., 100,-5. , 5.);
+// TH2F *hDE2 = new TH2F("hDE2" ,"EMC E-resolution", 100, 0., 10., 100, 0. , 5.);
+
+// cout << "Start EMC Analysis"<< endl ;
+// for (Int_t ievent=0; ievent<Nevents; ievent++) {
- //========== Event Number>
- if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
- cout << "==== AnalyzeEMC ====> Event is " << ievent+1 << endl ;
+// //========== 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; itrack<ntracks; itrack++) {
- // Get the Hits Tree for the Primary track itrack
- gAlice->ResetHits();
- gAlice->TreeH()->GetEvent(itrack);
- for (Int_t iModule=0; iModule < nOfModules; iModule++) {
- emcModule = fPHOS->GetEMCModule(iModule);
- emcHits = emcModule.Hits();
- nEMChits = emcHits->GetEntriesFast();
- for (Int_t ihit=0; ihit<nEMChits; ihit++) {
- emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
- TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
- new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*emcHit);
- }
- emcModule.Clear();
- }
- }
+// 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; ihit<nEMChits; ihit++) {
- emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
- genX = emcHit->X();
- genZ = emcHit->Y();
- genE = emcHit->GetMomentum().E();
- }
- Float_t dx = recX - genX;
- Float_t dz = recZ - genZ;
- Float_t de = recE - genE;
- hDx1 ->Fill(dx);
- hDz1 ->Fill(dz);
- hDE1 ->Fill(de);
- hDx2 ->Fill(genE,dx);
- hDz2 ->Fill(genE,dz);
- hDE2 ->Fill(genE,recE);
- }
- }
- delete [] hitsPerModule;
- }
- // Save histograms
+// 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; itrack<ntracks; itrack++) {
+// // Get the Hits Tree for the Primary track itrack
+// gAlice->ResetHits();
+// gAlice->TreeH()->GetEvent(itrack);
+// for (Int_t iModule=0; iModule < nOfModules; iModule++) {
+// emcModule = fPHOS->GetEMCModule(iModule);
+// emcHits = emcModule.Hits();
+// nEMChits = emcHits->GetEntriesFast();
+// for (Int_t ihit=0; ihit<nEMChits; ihit++) {
+// emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
+// TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
+// new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*emcHit);
+// }
+// emcModule.Clear();
+// }
+// }
+
+// // Loop over reconstructed particles
+
+// TClonesArray ** recParticleList = fPHOS->RecParticles() ;
+// AliPHOSRecParticle * recParticle ;
+// Int_t nEMCrecs = (*recParticleList)->GetEntries();
+// if (nEMCrecs == 1) {
+// recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(0) ;
+// Float_t recE = recParticle->Energy();
+// Int_t phosModule;
+// Double_t recX, recZ ;
+// fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), phosModule, recX, recZ) ;
+
+// // for this rec.point take the hit list in the same PHOS module
+
+// emcHits = hitsPerModule[phosModule-1];
+// Int_t nEMChits = emcHits->GetEntriesFast();
+// if (nEMChits == 1) {
+// Float_t genX, genZ, genE;
+// for (Int_t ihit=0; ihit<nEMChits; ihit++) {
+// emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
+// genX = emcHit->X();
+// genZ = emcHit->Y();
+// genE = emcHit->GetMomentum().E();
+// }
+// Float_t dx = recX - genX;
+// Float_t dz = recZ - genZ;
+// Float_t de = recE - genE;
+// hDx1 ->Fill(dx);
+// hDz1 ->Fill(dz);
+// hDE1 ->Fill(de);
+// hDx2 ->Fill(genE,dx);
+// hDz2 ->Fill(genE,dz);
+// hDE2 ->Fill(genE,recE);
+// }
+// }
+// delete [] hitsPerModule;
+// }
+// // Save histograms
+
+// Text_t outputname[80] ;
+// sprintf(outputname,"%s.analyzed",fRootFile->GetName());
+// TFile output(outputname,"RECREATE");
+// output.cd();
+
+// hDx1 ->Write() ;
+// hDz1 ->Write() ;
+// hDE1 ->Write() ;
+// hDx2 ->Write() ;
+// hDz2 ->Write() ;
+// hDE2 ->Write() ;
+
+// // Plot histograms
+
+// TCanvas *emcCanvas = new TCanvas("EMC","EMC analysis",20,20,700,300);
+// gStyle->SetOptStat(111111);
+// gStyle->SetOptFit(1);
+// gStyle->SetOptDate(1);
+// emcCanvas->Divide(3,1);
+
+// emcCanvas->cd(1);
+// gPad->SetFillColor(10);
+// hDx1->SetFillColor(16);
+// hDx1->Draw();
+
+// emcCanvas->cd(2);
+// gPad->SetFillColor(10);
+// hDz1->SetFillColor(16);
+// hDz1->Draw();
+
+// emcCanvas->cd(3);
+// gPad->SetFillColor(10);
+// hDE1->SetFillColor(16);
+// hDE1->Draw();
+
+// emcCanvas->Print("EMC.ps");
}
}
//____________________________________________________________________________
-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
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 ; }
}
//____________________________________________________________________________
-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;
added.fELOS += rValue.GetEnergy() ;
+ if(fPid == 0) fPid = rValue.fPid ;
+
return added;
}
#include <iostream.h>
+class TLorentzVector ;
+
class AliPHOSHit : public AliHit {
friend ostream& operator << (ostream&, const AliPHOSHit&) ;
// 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
}
// 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 ;
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
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);
+ // }
}
//
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)
// Create array of EMC modules of the size of PHOS modules number
- if ( 0==(fEMCModules=new TClonesArray("AliPHOSCPVModule",fGeom->GetNModules())) ) {
- Error("AliPHOSv1","Can not create array of EMC modules");
- exit(1);
- }
- TClonesArray &lemcmodule = *fEMCModules;
- for (Int_t i=0; i<fGeom->GetNModules(); i++) new(lemcmodule[i]) AliPHOSCPVModule();
+// if ( 0==(fEMCModules=new TClonesArray("AliPHOSCPVModule",fGeom->GetNModules())) ) {
+// Error("AliPHOSv1","Can not create array of EMC modules");
+// exit(1);
+// }
+// TClonesArray &lemcmodule = *fEMCModules;
+// for (Int_t i=0; i<fGeom->GetNModules(); i++) new(lemcmodule[i]) AliPHOSCPVModule();
// Create array of CPV modules for the IHEP's version of CPV
- if ( strcmp(fGeom->GetName(),"IHEP") == 0 || 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; i<fGeom->GetNCPVModules(); 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; i<fGeom->GetNCPVModules(); 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);
- }
+// }
}
//____________________________________________________________________________
fHits = 0 ;
}
+ if ( fSDigits) {
+ fSDigits->Delete() ;
+ delete fSDigits ;
+ fSDigits = 0 ;
+ }
+
if ( fDigits) {
fDigits->Delete() ;
delete fDigits ;
}
//____________________________________________________________________________
-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
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] ;
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; itrack<ntracks; itrack++){
-
+ for (itrack=0; itrack<gAlice->GetNtrack(); 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) ;
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 ;
}
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 << " " <<idCurSDigit<< endl ;
+
+ Int_t absID ;
+ for(absID = 1; absID < fGeom->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) ;
//___________________________________________________________________________
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<i> for hits in EMC modules
-
- for( Int_t i=0; i<fGeom->GetNModules(); i++ ) GetEMCModule(i).MakeBranch("EMC",i+1,file);
-
- // Create new branches CPV<i> for hits in CPV modules for IHEP geometry
-
- if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
- for( Int_t i=0; i<fGeom->GetNCPVModules(); i++ ) GetCPVModule(i).MakeBranch("CPV",i+1,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) ;
+ }
+
+ }
+
+
+
}
//_____________________________________________________________________________
// 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.
// Yuri Kharlov, 28 September 2000
AliDetector::ResetHits();
- for (Int_t i=0; i<fGeom->GetNModules(); i++) ((AliPHOSCPVModule*)(*fEMCModules)[i]) -> Clear();
- if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
- for (Int_t i=0; i<fGeom->GetNCPVModules(); i++) ((AliPHOSCPVModule*)(*fCPVModules)[i]) -> Clear();
- }
+ // for (Int_t i=0; i<fGeom->GetNModules(); i++) ((AliPHOSCPVModule*)(*fEMCModules)[i]) -> Clear();
+ // if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
+ // for (Int_t i=0; i<fGeom->GetNCPVModules(); i++) ((AliPHOSCPVModule*)(*fCPVModules)[i]) -> Clear();
+ // }
}
//____________________________________________________________________________
}
+//____________________________________________________________________________
+//void AliPHOSv1::SDigits2Digits() {
+// // Adds the noise to the summable digits and keeps digits above a threshold
+// // To make a digit
+//}
+
//____________________________________________________________________________
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; i<fGeom->GetNModules(); i++ ) GetEMCModule(i).SetTreeAddress("EMC",i+1);
+ // for( Int_t i=0; i<fGeom->GetNModules(); i++ ) GetEMCModule(i).SetTreeAddress("EMC",i+1);
// Set branch address for the Hits Tree for hits in CPV modules for IHEP geometry
// Yuri Kharlov, 28 September 2000.
- if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
- for( Int_t i=0; i<fGeom->GetNCPVModules(); i++ ) GetCPVModule(i).SetTreeAddress("CPV",i+1);
- }
+ // if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
+ // for( Int_t i=0; i<fGeom->GetNCPVModules(); i++ ) GetCPVModule(i).SetTreeAddress("CPV",i+1);
+ // }
}
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
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();
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
// 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);
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);
}
} // 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[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();
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
+
+
}
//____________________________________________________________________________
}
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){
// 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) ;
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
}
//____________________________________________________________________________
-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++ ;
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 ;
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; idigit<ndigits-1; idigit++) {
+ AliPHOSCPVDigit *cpvDigit1 = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(idigit);
+ Float_t x1 = cpvDigit1->GetXpad() ;
+ Float_t z1 = cpvDigit1->GetYpad() ;
+ for (Int_t jdigit=idigit+1; jdigit<ndigits; jdigit++) {
+ AliPHOSCPVDigit *cpvDigit2 = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(jdigit);
+ Float_t x2 = cpvDigit2->GetXpad() ;
+ Float_t z2 = cpvDigit2->GetYpad() ;
+ if (x1==x2 && z1==z2) {
+ Float_t qsum = cpvDigit1->GetQpad() + cpvDigit2->GetQpad() ;
+ cpvDigit2->SetQpad(qsum) ;
+ cpvDigits->RemoveAt(idigit) ;
+ }
+ }
+ }
+ cpvDigits->Compress() ;
+
+ // 3. add digits to temporary hit list fTmpHits
+
+ ndigits = cpvDigits->GetEntriesFast();
+ for (idigit=0; idigit<ndigits; idigit++) {
+ AliPHOSCPVDigit *cpvDigit = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(idigit);
+ relid[0] = moduleNumber + 1 ; // CPV (or PHOS) module number
+ relid[1] =-1 ; // means CPV
+ relid[2] = cpvDigit->GetXpad() ; // column number of a pad
+ relid[3] = cpvDigit->GetYpad() ; // row number of a pad
+
+ // get the absolute Id number
+ 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
}