From: dibari Date: Tue, 23 Aug 2005 09:56:04 +0000 (+0000) Subject: Possibility to run standalone PatRec without run again ESD and creation of ntuples... X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=910735ae7857719da2aa193d4969def28753d72a;p=u%2Fmrichter%2FAliRoot.git Possibility to run standalone PatRec without run again ESD and creation of ntuples for PatRec debugging + minor changes --- diff --git a/RICH/AliRICH.cxx b/RICH/AliRICH.cxx index a73a87ab34b..b9981d22022 100644 --- a/RICH/AliRICH.cxx +++ b/RICH/AliRICH.cxx @@ -618,13 +618,13 @@ void AliRICH::OccupancyPrint(Int_t iEvtNreq)const Int_t totPadsPerChamber = AliRICHParam::NpadsX()*AliRICHParam::NpadsY(); - Int_t nDigCh[kNchambers]={0,0,0,0,0,0,0}; - Int_t iChHits[kNchambers]={0,0,0,0,0,0,0}; - Int_t nPrim[kNchambers]={0,0,0,0,0,0,0}; - Int_t nSec[kNchambers]={0,0,0,0,0,0,0}; - for(Int_t iEvtN=iEvtNmin;iEvtNGetRunLoader()->GetEvent(iEvtN)) return; AliStack *pStack = GetLoader()->GetRunLoader()->Stack(); for(Int_t iPrimN=0;iPrimNTreeH()->GetEntries();iPrimN++){//prims loop @@ -638,13 +638,13 @@ void AliRICH::OccupancyPrint(Int_t iEvtNreq)const } } GetLoader()->TreeD()->GetEntry(0); - for(Int_t iChamber=1;iChamber<=kNchambers;iChamber++) nDigCh[iChamber-1]= Digits(iChamber)->GetEntries(); - } - for(Int_t iChamber=1;iChamber<=kNchambers;iChamber++){ - Double_t occupancy = (Double_t)nDigCh[iChamber-1]/(Double_t)totPadsPerChamber; - Info("Occupancy","for chamber %i = %4.2f %% and charged prim tracks %i and sec. tracks %i with total %i", + for(Int_t iChamber=1;iChamber<=kNchambers;iChamber++) { + nDigCh[iChamber-1]= Digits(iChamber)->GetEntries(); + Double_t occupancy = (Double_t)nDigCh[iChamber-1]/(Double_t)totPadsPerChamber; + Info("Occupancy","for chamber %i = %4.2f %% and charged prim tracks %i and sec. tracks %i with total %i", iChamber,occupancy*100.,nPrim[iChamber-1],nSec[iChamber-1],iChHits[iChamber-1]); - } + } + } GetLoader()->UnloadHits(); GetLoader()->UnloadDigits(); GetLoader()->GetRunLoader()->UnloadHeader(); @@ -879,6 +879,67 @@ void AliRICH::CheckPR()const pFile->Write();pFile->Close(); } //__________________________________________________________________________________________________ +void AliRICH::RichAna() +{ + TFile *pFile=TFile::Open("AliESDs.root","read"); + if(!pFile || !pFile->IsOpen()) {AliInfo("ESD file not open.");return;} //open AliESDs.root + TTree *pTree = (TTree*) pFile->Get("esdTree"); + if(!pTree){AliInfo("ESD not found.");return;} //get ESD tree + AliInfo("ESD found. Go ahead!"); + + AliMagF * magf = gAlice->Field(); + AliTracker::SetFieldMap(magf,kTRUE); +// AliTracker::SetFieldMap(magf); + + TFile *pFileRA = new TFile("$(HOME)/RichAna.root","RECREATE","RICH Pattern Recognition"); + TNtupleD *hn = new TNtupleD("hn","ntuple","Pmod:Charge:TrackTheta:TrackPhi:MinX:MinY:ThetaCerenkov:NPhotons:ChargeMIP:Chamber:TOF:LengthTOF:prob1:prob2:prob3"); + + AliESD *pESD=new AliESD; pTree->SetBranchAddress("ESD", &pESD); +// for(Int_t iEvtN=0;iEvtNGetRunLoader()->GetNumberOfEvents();iEvtN++) { + for(Int_t iEvtN=0;iEvtN<2;iEvtN++) { + pTree->GetEvent(iEvtN); + AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH")); + pRich->GetLoader()->GetRunLoader()->GetEvent(iEvtN); + pRich->GetLoader()->LoadRecPoints(); + pRich->GetLoader()->TreeR()->GetEntry(0); +//Pattern recognition started + if(pESD->GetNumberOfTracks()) { + Int_t iNtracks=pESD->GetNumberOfTracks(); + Info("RichAna",Form("Start with %i tracks",iNtracks)); + for(Int_t iTrackN=0;iTrackNRecWithESD(pESD,pRich,iTrackN); + AliESDtrack *pTrack = pESD->GetTrack(iTrackN);// get next reconstructed track + Double_t dx,dy; + Double_t hnvec[20]; + pTrack->GetRICHdxdy(dx,dy); + hnvec[0]=pTrack->GetP(); + hnvec[1]=pTrack->GetSign(); +// cout << " Track momentum " << pTrack->GetP() << " charge " << pTrack->GetSign() << endl; + + pTrack->GetRICHthetaPhi(hnvec[2],hnvec[3]); + pTrack->GetRICHdxdy(hnvec[4],hnvec[5]); + hnvec[6]=pTrack->GetRICHsignal(); + hnvec[7]=pTrack->GetRICHnclusters(); + hnvec[9]=pTrack->GetRICHcluster()/1000000; + hnvec[8]=pTrack->GetRICHcluster()-hnvec[9]*1000000; + hnvec[10]=pTrack->GetTOFsignal(); + hnvec[11]=pTrack->GetIntegratedLength(); + Double_t prob[5]; + pTrack->GetRICHpid(prob); + hnvec[12]=prob[0]+prob[1]+prob[2]; + hnvec[13]=prob[3]; + hnvec[14]=prob[4]; + hn->Fill(hnvec); + } + } + Info("RichAna","Pattern Recognition done for event %i",iEvtN); + } + pFileRA->Write();pFileRA->Close();// close RichAna.root + delete pESD; pFile->Close();//close AliESDs.root +} +//__________________________________________________________________________________________________ void AliRICH::DisplayEvent(Int_t iEvtNmin,Int_t iEvtNmax)const { TH2F *pDigitsH2[8]; @@ -980,7 +1041,7 @@ void AliRICH::Display()const pHitsH2->SetMarkerColor(kRed); pHitsH2->SetMarkerStyle(29); pHitsH2->SetMarkerSize(0.4); ReadESD(iEventN,iChamber); pHitsH2->Draw(); - ReadESD(iEventN,iChamber); +// ReadESD(iEventN,iChamber); AliRICHParam::DrawSectors(); TLatex l; l.SetNDC(); l.SetTextSize(0.02); if(!isHits) {l.SetTextColor(kRed) ;l.DrawLatex(0.1,0.01,"No Hits" );} diff --git a/RICH/AliRICH.h b/RICH/AliRICH.h index 10cc5842029..be8a3f3081f 100644 --- a/RICH/AliRICH.h +++ b/RICH/AliRICH.h @@ -84,6 +84,7 @@ public: void PrintTracks (Int_t iEvent=0); //prints a list of tracks for a given event void CheckPR ()const; //utility-> run staff for stack ?????? void ReadESD(Int_t iEventN, Int_t iChamber)const; + void RichAna (); //utility-> reconstruct ESD rings outside ESD void DrawRing(TVector3 entrance,TVector3 vectorTrack,Double_t thetaCer)const; protected: diff --git a/RICH/AliRICHRecon.cxx b/RICH/AliRICHRecon.cxx index b3021fda029..4466b3c18be 100644 --- a/RICH/AliRICHRecon.cxx +++ b/RICH/AliRICHRecon.cxx @@ -25,6 +25,7 @@ #include #include #include +#include #include "AliRICH.h" #include "AliRICHParam.h" @@ -40,16 +41,18 @@ AliRICHRecon::AliRICHRecon(AliRICHHelix *pHelix,TClonesArray *pClusters,Int_t iM { // main ctor SetFreonScaleFactor(1); - fIsWEIGHT = kFALSE; +// fIsWEIGHT = kFALSE; + fIsWEIGHT = kTRUE; fThetaBin=750; fThetaMin = 0.0; fThetaMax = 0.75; fDTheta = 0.001; fWindowWidth = 0.045; + fMinNumPhots = 3; fRadiatorWidth = AliRICHParam::Zfreon(); fQuartzWidth = AliRICHParam::Zwin(); fGapWidth = AliRICHParam::Freon2Pc() - fRadiatorWidth - fQuartzWidth; - fXmin = -AliRICHParam::PcSizeX()/2.; - fXmax = AliRICHParam::PcSizeX()/2.; - fYmin = -AliRICHParam::PcSizeY()/2.; - fYmax = AliRICHParam::PcSizeY()/2.; + fXmin = 0.; + fXmax = AliRICHParam::PcSizeX(); + fYmin = 0.; + fYmax = AliRICHParam::PcSizeY(); SetTrackTheta(pHelix->Ploc().Theta()); SetTrackPhi(pHelix->Ploc().Phi()); SetMipIndex(iMipId); @@ -476,6 +479,7 @@ void AliRICHRecon::FindAreaAndPortionOfRing() Int_t zone = CheckDetectorAcceptance(); + AliDebug(1,Form("acceptance to detector zone -> %d",zone)); if (zone != 0) { @@ -792,176 +796,54 @@ Float_t AliRICHRecon::SnellAngle(Float_t n1, Float_t n2, Float_t theta1) //__________________________________________________________________________________________________ void AliRICHRecon::HoughResponse() { - //Hough response - -// Implement Hough response pat. rec. method - - Float_t *hCSspace; - - int bin=0; - int bin1=0; - int bin2=0; - int i, j, k, nCorrBand; - float hcs[750],hcsw[750]; - float angle, weight; - float lowerlimit,upperlimit; - - float etaPeak[100]; - - int nBin; - - float etaPeakPos = -1; - - Int_t etaPeakCount = -1; - - Float_t thetaCerenkov = 0.; - - nBin = (int)(0.5+fThetaMax/(fDTheta)); - nCorrBand = (int)(0.5+ fWindowWidth/(2 * fDTheta)); - - memset ((void *)hcs, 0, fThetaBin*sizeof(float)); - memset ((void *)hcsw, 0, fThetaBin*sizeof(float)); - + TH1F *phots = new TH1F("phots","phots",750,0.,0.75); + TH1F *photsw = new TH1F("photsw","photsw",750,0.,0.75); + TH1F *resultw = new TH1F("resultw","resultw",750,0.,0.75); Int_t nPhotons = GetPhotonsNumber(); - - Int_t weightFlag = 0; - - for (k=0; k< nPhotons; k++) { - - SetPhotonIndex(k); - - angle = GetPhotonEta(); - + Int_t nBin = (Int_t)(fThetaMax/fDTheta); + Int_t nCorrBand = (Int_t)(fWindowWidth/(2*fDTheta)); + AliDebug(1,Form("Ring reconstruction for track with theta %f",GetTrackTheta()*TMath::RadToDeg())); + for (Int_t kPhot=0; kPhot< nPhotons; kPhot++){ + SetPhotonIndex(kPhot); + Double_t angle = GetPhotonEta(); if(angle == -999.) continue; - - if (angle>=fThetaMin && angle<= fThetaMax) - - { - - bin = (int)(0.5+angle/(fDTheta)); - - bin1= bin-nCorrBand; - bin2= bin+nCorrBand; - - // calculate weights - - if(fIsWEIGHT) - { - lowerlimit = ((Float_t)bin1)*fDTheta + 0.5*fDTheta; - SetThetaOfRing(lowerlimit); - FindAreaAndPortionOfRing(); - Float_t area1 = GetAreaOfRing(); - - upperlimit = ((Float_t)bin2)*fDTheta + 0.5*fDTheta; - SetThetaOfRing(upperlimit); - FindAreaAndPortionOfRing(); - Float_t area2 = GetAreaOfRing(); - - // cout << "lowerlimit" << lowerlimit << "upperlimit " << upperlimit << endl; - Float_t diffarea = area2 - area1; - - if(diffarea>0) - { - weight = 1./(area2-area1); - } - else - { - weightFlag = 1; - weight = 1.; - } - - // cout <<" low "<< lowerlimit << " up " << upperlimit << - // " area1 " << area1 << " area2 " << area2 << " weight " << weight << endl; - - } - else - { - weight = 1.; - } - - SetPhotonWeight(weight); - - // cout << "weight..." << weight << endl; - - - if (bin1<0) bin1=0; - if (bin2>nBin) bin2=nBin; - - for (j=bin1; j etaPeakPos) { - etaPeakCount = 0; - etaPeakPos = hCSspace[bin]; - etaPeak[0]=angle; - } - else { - if (hCSspace[bin] == etaPeakPos) { - etaPeak[++etaPeakCount] = angle; - } - } - } - - for (i=0; i=0) { - thetaCerenkov /= etaPeakCount+1; - fThetaPeakPos = etaPeakPos; - } - - SetThetaCerenkov(thetaCerenkov); -} -//__________________________________________________________________________________________________ -void AliRICHRecon::HoughFiltering(float hcs[]) -{ - // filter for Hough - -// hough filtering - - float hcsFilt[750]; - float k[5] = {0.05, 0.25, 0.4, 0.25, 0.05}; - int nx, i, nxDx; - int sizeHCS; - int nBin; - - nBin = (int)(1+fThetaMax/fDTheta); - sizeHCS = fThetaBin*sizeof(float); - - memset ((void *)hcsFilt, 0, sizeHCS); - - for (nx = 0; nx < nBin; nx++) { - for (i = 0; i < 5; i++) { - nxDx = nx + (i-2); - if (nxDx> -1 && nxDxFill(angle); + Int_t bin = (Int_t)(0.5+angle/(fDTheta)); + Double_t weight=1.; + if(fIsWEIGHT){ + Double_t lowerlimit = ((Float_t)bin)*fDTheta - 0.5*fDTheta; + SetThetaOfRing(lowerlimit); + FindAreaAndPortionOfRing(); + Float_t area1 = GetAreaOfRing(); + Double_t upperlimit = ((Float_t)bin)*fDTheta + 0.5*fDTheta; + SetThetaOfRing(upperlimit); + FindAreaAndPortionOfRing(); + Float_t area2 = GetAreaOfRing(); + AliDebug(1,Form("lowerlimit %f area %f ; upperlimit %f area %f",lowerlimit,area1,upperlimit,area2)); + Float_t diffarea = area2 - area1; + if(diffarea>0){weight = 1./(area2-area1);}else{weight = 1.;} + } + photsw->Fill(angle,weight); + SetPhotonWeight(weight); + } + for (Int_t i=1; i<=nBin;i++){ + Int_t bin1= i-nCorrBand; + Int_t bin2= i+nCorrBand; + if(bin1<1) bin1=1; + if(bin2>nBin)bin2=nBin; + Double_t sumPhots=phots->Integral(bin1,bin2); + if(sumPhotsIntegral(bin1,bin2); + resultw->Fill((Float_t)(i*fDTheta),sumPhotsw); +} +// evaluate the "BEST" thetacerenkov.... + Float_t *pVec = resultw->GetArray(); + Int_t locMax = TMath::LocMax(nBin,pVec); + SetThetaCerenkov((Double_t)(locMax*fDTheta+0.5*fDTheta)); + +// Reset and delete objects + phots->Delete();photsw->Delete();resultw->Delete(); +}//HoughResponse //__________________________________________________________________________________________________ void AliRICHRecon::FindThetaCerenkov() { @@ -1039,4 +921,3 @@ void AliRICHRecon::FlagPhotons() } SetHoughPhotons(nPhotonHough); } - diff --git a/RICH/AliRICHRecon.h b/RICH/AliRICHRecon.h index 0c520b5ef36..e1d383b1bef 100644 --- a/RICH/AliRICHRecon.h +++ b/RICH/AliRICHRecon.h @@ -160,6 +160,7 @@ protected: Int_t fCandidatePhotonsNumber; // number of candidate photons Int_t fHoughPhotons; // n. photons after Hough Int_t fFittedHoughPhotons; // n. photons after Hough and after minimization + Int_t fMinNumPhots; // minimum number of photons for a given ring Float_t fTrackTheta; // Theta of track at RICH Float_t fTrackPhi; // Phi of track at RICH diff --git a/RICH/AliRICHTracker.cxx b/RICH/AliRICHTracker.cxx index 8d03a9bd64d..6ee209e26e7 100644 --- a/RICH/AliRICHTracker.cxx +++ b/RICH/AliRICHTracker.cxx @@ -19,36 +19,35 @@ Int_t AliRICHTracker::PropagateBack(AliESD *pESD) // 1. AliESD - reconstructed tracks are used // 2. RICH private ntuple for debug- stack particles used instead of reconstructed tracks AliDebug(1,"Start pattern recognition"); - if(pESD->GetNumberOfTracks()) - RecWithESD(pESD); - else - RecWithStack(0); + if(pESD->GetNumberOfTracks()) { + Int_t iNtracks=pESD->GetNumberOfTracks(); + AliDebug(1,Form("Start with %i tracks",iNtracks)); + AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH")); + for(Int_t iTrackN=0;iTrackNGetNumberOfTracks(); - Double_t b=GetFieldMap()->SolenoidField()/10;// magnetic field in Tesla - AliDebug(1,Form("Start with %i tracks in %f Tesla field",iNtracks,b)); - - AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH")); - - for(Int_t iTrackN=0;iTrackNSolenoidField()/10;// magnetic field in Tesla AliESDtrack *pTrack = pESD->GetTrack(iTrackN);// get next reconstructed track // if((pTrack->GetStatus()&AliESDtrack::kTOFout)==0) continue; //ignore tracks not recontructed by TOF // pTrack->GetXYZ(xb); // pTrack->GetPxPyPz(pb); Int_t status=pTrack->GetStatus()&AliESDtrack::kTOFout;//get running track parameters - Int_t charge = (Int_t)(-TMath::Sign(1.,pTrack->GetSign()*b)); + Int_t charge = (Int_t)(-TMath::Sign(1.,pTrack->GetSign()*fField)); AliDebug(1,Form("Track %i pmod=%f charge=%i stat=%i",iTrackN,pTrack->GetP(),charge,status)); - AliRICHHelix helix(pTrack->X3(),pTrack->P3(),charge,b); + AliRICHHelix helix(pTrack->X3(),pTrack->P3(),charge,fField); Int_t iChamber=helix.RichIntersect(pRich->P()); AliDebug(1,Form("intersection with %i chamber found",iChamber)); - if(!iChamber) continue;//intersection with no chamber found + if(!iChamber) return;//intersection with no chamber found //find MIP cluster candidate (cluster which is closest to track intersection point) Double_t distMip=9999,distX=0,distY=0; //min distance between clusters and track position on PC Int_t iMipId=0; //index of that min distance cluster @@ -74,14 +73,14 @@ void AliRICHTracker::RecWithESD(AliESD *pESD) if(distMip>1||chargeMip<100) { //track not accepted for pattern recognition pTrack->SetRICHsignal(-999.); //to be improved by flags... - continue; + return; } // AliRICHRecon recon(&helix,pRich->Clusters(iChamber),iMipId); //actual job is done there Double_t thetaCerenkov=recon.ThetaCerenkov(); //search for mean Cerenkov angle for this track - pTrack->SetRICHcluster(iMipId+100000*iChamber); + pTrack->SetRICHcluster(((Int_t)chargeMip)+1000000*iChamber); pTrack->SetRICHdxdy(distX,distY); pTrack->SetRICHthetaPhi(helix.Ploc().Theta(),helix.Ploc().Phi()); pTrack->SetRICHsignal(thetaCerenkov); @@ -95,13 +94,11 @@ void AliRICHTracker::RecWithESD(AliESD *pESD) Double_t richPID[AliPID::kSPECIES]; for (Int_t iPart=0;iPartClusters(iChamber)->GetEntries();iphot++) { recon.SetPhotonIndex(iphot); if(recon.GetPhotonFlag() == 2) { - Double_t sigma2 = AliRICHParam::SigmaSinglePhoton(iPart,pTrack->GetP(),recon.GetTrackTheta(),recon.GetPhiPoint()-recon.GetTrackPhi()).Mag2(); - if (sigma2 >0) - sigmaPID[iPart] += 1/sigma2; + Double_t sigma = AliRICHParam::SigmaSinglePhoton(iPart,pTrack->GetP(),recon.GetTrackTheta(),recon.GetPhiPoint()-recon.GetTrackPhi()).Mag(); + sigmaPID[iPart] += 1/(sigma*sigma); } } if (sigmaPID[iPart]>0) @@ -111,10 +108,7 @@ void AliRICHTracker::RecWithESD(AliESD *pESD) CalcProb(thetaCerenkov,pTrack->GetP(),sigmaPID,richPID); pTrack->SetRICHpid(richPID); AliDebug(1,Form("PROBABILITIES ---> %f - %f - %f - %f - %f",richPID[0],richPID[1],richPID[2],richPID[3],richPID[4])); - } - - - }//ESD tracks loop + } AliDebug(1,"Stop."); } //RecWithESD //__________________________________________________________________________________________________ diff --git a/RICH/AliRICHTracker.h b/RICH/AliRICHTracker.h index 52ae65fc3bd..ec6e99a636f 100644 --- a/RICH/AliRICHTracker.h +++ b/RICH/AliRICHTracker.h @@ -7,6 +7,7 @@ class AliCluster; class AliESD; +class AliRICH; class TTree; class AliRICHTracker : public AliTracker @@ -19,13 +20,15 @@ public: void UnloadClusters() {AliDebug(1,"Start.");} //pure virtual from AliTracker AliCluster *GetCluster(Int_t )const {AliDebug(1,"Start.");return 0;} //pure virtual from AliTracker Int_t PropagateBack(AliESD *); //pure virtual from AliTracker - void RecWithESD(AliESD *); //recon with ESD + void RecWithESD(AliESD *,AliRICH *,Int_t iTrackN); //recon with ESD void RecWithStack(TNtupleD *hn); //recon from Stack in case ESD empty void CalcProb(Double_t thetaCer,Double_t pmod,Double_t *sigmaPID, Double_t *richPID); // calculate pid for RICH Int_t LoadClusters(TTree *); //pure virtual from AliTracker protected: + Double_t fField; // magnetic field stored + ClassDef(AliRICHTracker,0) };//class AliRICHTracker diff --git a/RICH/RichMenu.C b/RICH/RichMenu.C index f01946ce29b..66a239ce26a 100644 --- a/RICH/RichMenu.C +++ b/RICH/RichMenu.C @@ -106,7 +106,8 @@ void MenuRich() pMenu->AddButton("Print occupancy" ,"r->OccupancyPrint(-1);" ,"????"); pMenu->AddButton("Event Summary " ,"r->SummaryOfEvent();" ,"????"); pMenu->AddButton("Hits plots" ,"r->ControlPlots()" ,"????"); - pMenu->AddButton("Recon with stack" ,"r->CheckPR()" , "Create RSR.root with ntuple hn"); + pMenu->AddButton("Recon with stack" ,"r->CheckPR()", "Create RSR.root with ntuple hn"); + pMenu->AddButton("RichAna" ,"r->RichAna()", "Create RichAna.root with ntuple hn"); pMenu->Show(); }//TestMenu() //__________________________________________________________________________________________________