#include "TTree.h"
#include "TMath.h"
#include "TCanvas.h"
+#include "TStyle.h"
// --- Standard library ---
#include "AliPHOSTrackSegment.h"
#include "AliPHOSRecParticle.h"
#include "AliPHOSIndexToObject.h"
+#include "AliPHOSCPVHit.h"
ClassImp(AliPHOSAnalyze)
fEvt = -999 ;
}
+ fDebugLevel = 0;
fClu = 0 ;
fPID = 0 ;
fTrs = 0 ;
fRootFile->Close() ;
if(fRootFile)
delete fRootFile ;
-
- if(fPHOS)
- delete fPHOS ;
-
- if(fClu)
- delete fClu ;
-
- if(fPID)
- delete fPID ;
-
- if(fRec)
- delete fRec ;
-
- if(fTrs)
- delete fTrs ;
+ if(fRootFile) {delete fRootFile ; fRootFile=0 ;}
+ if(fPHOS) {delete fPHOS ; fPHOS =0 ;}
+ if(fClu) {delete fClu ; fClu =0 ;}
+ if(fPID) {delete fPID ; fPID =0 ;}
+ if(fRec) {delete fRec ; fRec =0 ;}
+ if(fTrs) {delete fTrs ; fTrs =0 ;}
}
//=========== Gets the Kine TTree
gAlice->TreeK()->GetEvent(0) ;
-
+
//=========== Get the Digit Tree
gAlice->TreeD()->GetEvent(0) ;
//____________________________________________________________________________
void AliPHOSAnalyze::Reconstruct(Int_t Nevents,Int_t FirstEvent )
{
- Int_t ievent ;
- for ( ievent=FirstEvent; ievent<Nevents; ievent++)
- {
- if (ievent==FirstEvent)
- {
- cout << "Analyze > Starting Reconstructing " << endl ;
- //========== Create the Clusterizer
- fClu = new AliPHOSClusterizerv1() ;
- fClu->SetEmcEnergyThreshold(0.03) ;
- fClu->SetEmcClusteringThreshold(0.20) ;
- fClu->SetPpsdEnergyThreshold (0.0000001) ;
- fClu->SetPpsdClusteringThreshold(0.0000001) ;
- fClu->SetLocalMaxCut(0.03) ;
- fClu->SetCalibrationParameters(0., 0.00000001) ;
-
- //========== Creates the track segment maker
- fTrs = new AliPHOSTrackSegmentMakerv1() ;
- // fTrs->UnsetUnfoldFlag() ;
-
- //========== Creates the particle identifier
- fPID = new AliPHOSPIDv1() ;
- fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ;
- fPID->SetDispersionCutOff(2.0) ;
- fPID->SetRelativeDistanceCut(3.) ;
-
- //========== Creates the Reconstructioner
- fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
- // fRec -> SetDebugReconstruction(kTRUE);
- }
+ // Perform reconstruction of EMC and CPV (GPS2 or IHEP) for <Nevents> events
+ // Yuri Kharlov. 19 October 2000
+
+ Int_t ievent ;
+ for ( ievent=FirstEvent; ievent<Nevents; ievent++) {
+ if (ievent==FirstEvent) {
+ cout << "Analyze > Starting Reconstructing " << endl ;
+ //========== Create the Clusterizer
+ fClu = new AliPHOSClusterizerv1() ;
+ fClu->SetEmcEnergyThreshold(0.05) ;
+ fClu->SetEmcClusteringThreshold(0.20) ;
+ fClu->SetLocalMaxCut(0.03) ;
+ if (strcmp(fGeom->GetName(),"GPS2") == 0) {
+ fClu->SetPpsdEnergyThreshold (0.0000002) ;
+ fClu->SetPpsdClusteringThreshold(0.0000001) ;
+ }
+ else if (strcmp(fGeom->GetName(),"IHEP") == 0) {
+ fClu->SetLocalMaxCutCPV(0.03) ;
+ fClu->SetLogWeightCutCPV(4.0) ;
+ fClu->SetPpsdEnergyThreshold (0.09) ;
+ }
+ fClu->SetCalibrationParameters(0., 0.00000001) ;
- //========== Event Number>
- // if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
- cout << "Reconstruct > Event is " << ievent << endl ;
+ //========== Creates the track segment maker
+ fTrs = new AliPHOSTrackSegmentMakerv1() ;
+ // fTrs->UnsetUnfoldFlag() ;
+
+ //========== Creates the particle identifier for GPS2 only
+ if (strcmp(fGeom->GetName(),"GPS2") == 0) {
+ fPID = new AliPHOSPIDv1() ;
+ fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ;
+ }
- //=========== Connects the various Tree's for evt
- gAlice->GetEvent(ievent);
+ //========== Creates the Reconstructioner
+ fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
+ if (fDebugLevel != 0) fRec -> SetDebugReconstruction(kTRUE);
+ }
+
+ if (fDebugLevel != 0 ||
+ (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
+ gAlice->GetEvent(ievent);
+
+ //=========== Gets the Digit TTree
+ gAlice->TreeD()->GetEvent(0) ;
+
+ //=========== Do the reconstruction
+ fPHOS->Reconstruction(fRec);
+ }
- //=========== Gets the Digit TTree
- gAlice->TreeD()->GetEvent(0) ;
+ if(fClu) {delete fClu ; fClu =0 ;}
+ if(fPID) {delete fPID ; fPID =0 ;}
+ if(fRec) {delete fRec ; fRec =0 ;}
+ if(fTrs) {delete fTrs ; fTrs =0 ;}
+
+}
- //=========== Do the reconstruction
- fPHOS->Reconstruction(fRec);
+//-------------------------------------------------------------------------------------
+void AliPHOSAnalyze::ReadAndPrintCPV(Int_t Nevents)
+{
+ //
+ // Read and print generated and reconstructed hits in CPV
+ // Author: Yuri Kharlov
+ // 12 October 2000
+ //
+
+ cout << "Start CPV Analysis"<< endl ;
+ for ( Int_t ievent=0; ievent<Nevents; ievent++) {
+
+ //========== Event Number>
+ cout << endl << "==== ReadAndPrintCPV ====> Event is " << ievent+1 << endl ;
+
+ //=========== Connects the various Tree's for evt
+ gAlice->GetEvent(ievent);
+
+ //=========== Get the Hits Tree
+ gAlice->ResetHits();
+ gAlice->TreeH()->GetEvent(0);
+
+ //========== Creating branches ===================================
+ AliPHOSRecPoint::RecPointsList ** EmcRecPoints = fPHOS->EmcRecPoints() ;
+ gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , EmcRecPoints ) ;
+
+ AliPHOSRecPoint::RecPointsList ** CpvRecPoints = fPHOS->PpsdRecPoints() ;
+ gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", CpvRecPoints ) ;
+
+ //=========== Gets the Reconstruction TTree
+ gAlice->TreeR()->GetEvent(0) ;
+
+ // Read and print CPV hits
+
+ TClonesArray *CPVhits;
+ for (Int_t iModule=0; iModule < fGeom->GetNModules(); iModule++) {
+ AliPHOSCPVModule cpvModule = fPHOS->GetCPVModule(iModule);
+ CPVhits = cpvModule.Hits();
+ Int_t nCPVhits = CPVhits->GetEntriesFast();
+ for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
+ AliPHOSCPVHit *cpvHit = (AliPHOSCPVHit*)CPVhits->UncheckedAt(ihit);
+ TLorentzVector p = cpvHit->GetMomentum();
+ Float_t xgen = cpvHit->GetX();
+ Float_t zgen = cpvHit->GetY();
+ Int_t ipart = cpvHit->GetIpart();
+ printf("CPV hit in module %d: ",iModule+1);
+ printf(" p = (%f, %f, %f, %f) GeV,\n",
+ p.Px(),p.Py(),p.Pz(),p.Energy());
+ printf(" xy = (%8.4f, %8.4f) cm, ipart = %d\n",
+ xgen,zgen,ipart);
+ }
+ }
+
+ // Read and print CPV reconstructed points
+
+ TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
+ AliPHOSPpsdRecPoint *cpvRecPoint ;
+ while( ( cpvRecPoint = (AliPHOSPpsdRecPoint *)nextRP() ) ) {
+ TVector3 locpos;
+ cpvRecPoint->GetLocalPosition(locpos);
+ Int_t PHOSModule = cpvRecPoint->GetPHOSMod();
+ printf("CPV recpoint in module %d: (X,Y,Z) = (%f,%f,%f) cm\n",
+ PHOSModule,locpos.X(),locpos.Y(),locpos.Z());
+ }
+ }
+}
+
+//____________________________________________________________________________
+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.);
+ 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 ;
+
+ //=========== Connects the various Tree's for evt
+ gAlice->GetEvent(ievent);
+ //=========== Get the Hits Tree
+ gAlice->ResetHits();
+ gAlice->TreeH()->GetEvent(0);
+
+ //========== Creating branches ===================================
+ AliPHOSRecPoint::RecPointsList ** EmcRecPoints = fPHOS->EmcRecPoints() ;
+ gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , EmcRecPoints ) ;
+
+ AliPHOSRecPoint::RecPointsList ** CpvRecPoints = fPHOS->PpsdRecPoints() ;
+ gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", CpvRecPoints ) ;
+ //=========== Gets the Reconstruction TTree
+ gAlice->TreeR()->GetEvent(0) ;
+ TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
+ AliPHOSCpvRecPoint *cpvRecPoint ;
+ AliPHOSCPVModule cpvModule;
+ TClonesArray *CPVhits;
+ 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;
+
+ cpvModule = fPHOS->GetCPVModule(PHOSModule-1);
+ CPVhits = cpvModule.Hits();
+ Int_t nCPVhits = CPVhits->GetEntriesFast();
+ for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
+ AliPHOSCPVHit *cpvHit = (AliPHOSCPVHit*)CPVhits->UncheckedAt(ihit);
+ Float_t xgen = cpvHit->GetX();
+ Float_t zgen = cpvHit->GetY();
+ if ( TMath::Abs(xgen-xrec) < TMath::Abs(dxmin) ) dxmin = xgen-xrec;
+ if ( TMath::Abs(zgen-zrec) < TMath::Abs(dzmin) ) dzmin = zgen-zrec;
+ }
+ cpvModule.Clear();
+ hDx ->Fill(dxmin);
+ hDz ->Fill(dzmin);
+ hNrp ->Fill(rpMult);
+ hNrpX->Fill(rpMultX);
+ hNrpZ->Fill(rpMultZ);
}
+ }
+ // Save histograms
+
+ Text_t outputname[80] ;
+ sprintf(outputname,"%s.analyzed",fRootFile->GetName());
+ TFile output(outputname,"RECREATE");
+ output.cd();
- fClu->Delete();
- fClu=0 ;
- fTrs->Delete();
- fTrs = 0 ;
- fPID->Delete();
- fPID = 0 ;
- fRec->Delete();
- fRec = 0 ;
+ hDx ->Write() ;
+ hDz ->Write() ;
+ hNrp ->Write() ;
+ hNrpX->Write() ;
+ hNrpZ->Write() ;
+
+ // Plot histograms
+
+ TCanvas *CPVcanvas = new TCanvas("CPV","CPV analysis",20,20,600,600);
+ gStyle->SetOptStat(111111);
+ gStyle->SetOptFit(1);
+ gStyle->SetOptDate(1);
+ CPVcanvas->Divide(3,3);
+
+ 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->Print("CPV.ps");
}
+
//____________________________________________________________________________
void AliPHOSAnalyze::InvariantMass(Int_t Nevents )
{
NRecParticles[ievent] = iRecPhot-1 ;
}
-
//Now calculate invariant mass:
Int_t irp1,irp2 ;
Int_t NCurEvent = 0 ;
} //loop over second rp
}//loop over first rp
-
-
AllRecParticleList->Delete() ;
} //Loop over events
TotalRPwithPrim++;
Int_t PrimaryType = ((TParticle *)PrimaryList->At(ClosestPrimary))->GetPdgCode() ;
- TParticlePDG* PDGparticle = ((TParticle *)PrimaryList->At(ClosestPrimary))->GetPDG();
- Double_t charge = PDGparticle->Charge() ;
+// TParticlePDG* PDGparticle = ((TParticle *)PrimaryList->At(ClosestPrimary))->GetPDG();
+// Double_t charge = PDGparticle->Charge() ;
// if(charge)
// cout <<"Primary " <<PrimaryType << " E " << ((TParticle *)PrimaryList->At(ClosestPrimary))->Energy() << endl ;
Int_t PrimaryCode ;
fhAllPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
fhVeto->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
fhPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
-
Counter[5][PrimaryCode]++ ;
break ;
case AliPHOSFastRecParticle::kABSURDEM:
//========== Creates the Reconstructioner
fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
- fRec -> SetDebugReconstruction(kFALSE);
+// fRec -> SetDebugReconstruction(kFALSE);
+ fRec -> SetDebugReconstruction(kTRUE);
//=========== Connect the various Tree's for evt
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff