TVector2 x2 = C(pHit->C())->Mrs2Pc(pHit->OutX3());//hit position in photocathode plane
Int_t iTotQdc=P()->TotQdc(x2,pHit->Eloss());//total charge produced by hit, 0 if hit in dead zone
if(iTotQdc==0) continue;
- TVector area=P()->Loc2Area(x2);//determine affected pads, dead zones analysed inside
- AliDebug(1,Form("hit(%6.2f,%6.2f)->area(%3.0f,%3.0f)-(%3.0f,%3.0f) QDC=%4i",x2.X(),x2.Y(),area[0],area[1],area[2],area[3],iTotQdc));
+ //
+ //need to quantize the anod....
+ TVector padHit=AliRICHParam::Loc2Pad(x2);
+ TVector2 padHitXY=AliRICHParam::Pad2Loc(padHit);
+ TVector2 anod;
+ if((x2.Y()-padHitXY.Y())>0) anod.Set(x2.X(),padHitXY.Y()+AliRICHParam::PitchAnod()/2);
+ else anod.Set(x2.X(),padHitXY.Y()-AliRICHParam::PitchAnod()/2);
+ //end to quantize anod
+ //
+ TVector area=P()->Loc2Area(anod);//determine affected pads, dead zones analysed inside
+ AliDebug(1,Form("hitanod(%6.2f,%6.2f)->area(%3.0f,%3.0f)-(%3.0f,%3.0f) QDC=%4i",anod.X(),anod.Y(),area[0],area[1],area[2],area[3],iTotQdc));
TVector pad(2);
for(pad[1]=area[1];pad[1]<=area[3];pad[1]++)//affected pads loop
for(pad[0]=area[0];pad[0]<=area[2];pad[0]++){
- Double_t padQdc=iTotQdc*P()->FracQdc(x2,pad);
+ Double_t padQdc=iTotQdc*P()->FracQdc(anod,pad);
AliDebug(1,Form("current pad(%3.0f,%3.0f) with QDC =%6.2f",pad[0],pad[1],padQdc));
if(padQdc>0.1) AddSDigit(pHit->C(),pad,padQdc,GetLoader()->GetRunLoader()->Stack()->Particle(pHit->GetTrack())->GetPdgCode(),pHit->GetTrack());
}//affected pads loop
void Fill(AliRICHcluster *pRaw,Double_t x,Double_t y,Double_t q,Int_t cfm) //form new resolved cluster from raw one
{fCFM=cfm;fChamber=pRaw->Fchamber();fSize=pRaw->Fsize();fQdc=(Int_t)(q*pRaw->Q());fX=x;fY=y;fStatus=kResolved;}
Double_t DistTo(TVector2 x) const{return TMath::Sqrt((x.X()-fX)*(x.X()-fX)+(x.Y()-fY)*(x.Y()-fY));} //distance to given point
+ Double_t DistX(TVector2 x) const{return (x.X()-fX);} //distance in x to given point
+ Double_t DistY(TVector2 x) const{return (x.Y()-fY);} //distance to given point
protected:
Int_t fCFM; //1000000*Ncerenkovs+1000*Nfeedbacks+Nmips
Int_t fSize; //10000*(how many digits belong to this cluster) + nLocalMaxima
// horizontal angle between chambers 19.5 grad
// vertical angle between chambers 20 grad
RotY(90);//rotate around y
- fCenterX3.SetXYZ(490,0,0);fPcX3.SetXYZ(490+8,0,0);fRadX3.SetXYZ(490-2,0,0); //shift center along x by 490 cm
+ fCenterX3.SetXYZ(490,0,0);fPcX3.SetXYZ(490+8-0.4,0,0);fRadX3.SetXYZ(490-2,0,0); //shift center along x by 490 cm
switch(iChamber){
case 0: //special test beam configuration without rotation.
//Finds cerenkov-feedback-mip mixture for a given cluster
AliDebug(1,"Start.");ToAliDebug(1,pCluster->Print());
- R()->GetLoader()->GetRunLoader()->LoadHeader();
+// R()->GetLoader()->GetRunLoader()->LoadHeader();
AliStack *pStack = R()->GetLoader()->GetRunLoader()->Stack();
if(!pStack)
{AliInfo("No Stack found!!! No contrib to cluster found.");return;}
#include <TPolyLine.h>
#include <TParticle.h>
#include <TStyle.h>
-#include <TH2.h>
#include <TMath.h>
#include <TLatex.h>
ClassImp(AliRICHDisplFast)
+//__________________________________________________________________________________________________
+void AliRICHDisplFast::ShowEvent(Int_t iEvtNmin,Int_t iEvtNmax)
+{
+ TH2F *pDigitsH2[8];
+ char titobj[11],titdisp[20];
+
+ AliRICH *pRich = (AliRICH*)gAlice->GetDetector("RICH");
+ Bool_t isDigits =!pRich->GetLoader()->LoadDigits();
+ if(!isDigits){Error("ShoEvent","No digits. Nothing to display.");return;}
+
+ TCanvas *canvas = new TCanvas("RICHDisplay","RICH Display",0,0,1226,900);
+// TCanvas *canvas = (TCanvas*)gROOT->FindObjectAny("RICHDisplay");
+
+ gStyle->SetPalette(1);
+ canvas->Divide(3,3);
+
+ for(Int_t iChamber=1;iChamber<=7;iChamber++) {
+ sprintf(titobj,"pDigitsH2_%i",iChamber);
+ sprintf(titdisp,"Chamber %i",iChamber);
+ pDigitsH2[iChamber] = new TH2F(titobj,titdisp,165,0,AliRICHParam::PcSizeX(),144,0,AliRICHParam::PcSizeY());
+ pDigitsH2[iChamber]->SetMarkerColor(kGreen);
+ pDigitsH2[iChamber]->SetMarkerStyle(29);
+ pDigitsH2[iChamber]->SetMarkerSize(0.4);
+ pDigitsH2[iChamber]->SetStats(kFALSE);
+ }
+
+ if(iEvtNmax>gAlice->GetEventsPerRun()) iEvtNmax=gAlice->GetEventsPerRun();
+
+ TText *t = new TText();
+ t->SetTextSize(0.10);
+
+ TText *tit = new TText(0.1,0.6,"RICH Display");
+ tit->SetTextSize(0.10);
+
+ for(Int_t iEventN=iEvtNmin;iEventN<=iEvtNmax;iEventN++) {
+
+ canvas->cd(1);
+ sprintf(titdisp,"Event Number %i",iEventN);
+ t->SetText(0.2,0.4,titdisp);
+ t->Draw();
+ tit->Draw();
+ pRich->GetLoader()->GetRunLoader()->GetEvent(iEventN);
+ pRich->GetLoader()->TreeD()->GetEntry(0);
+ for(Int_t iChamber=1;iChamber<=7;iChamber++) {
+ pDigitsH2[iChamber]->Reset();
+ for(Int_t j=0;j<pRich->Digits(iChamber)->GetEntries();j++) {//digits loop
+ AliRICHdigit *pDig = (AliRICHdigit*)pRich->Digits(iChamber)->At(j);
+ TVector2 x2=AliRICHParam::Pad2Loc(pDig->Pad());
+ pDigitsH2[iChamber]->Fill(x2.X(),x2.Y());
+ }
+ if(iChamber==1) canvas->cd(7);
+ if(iChamber==2) canvas->cd(8);
+ if(iChamber==3) canvas->cd(4);
+ if(iChamber==4) canvas->cd(5);
+ if(iChamber==5) canvas->cd(6);
+ if(iChamber==6) canvas->cd(2);
+ if(iChamber==7) canvas->cd(3);
+ pDigitsH2[iChamber]->Draw();
+ }
+ canvas->Update();
+ canvas->Modified();
+
+ if(iEvtNmin<iEvtNmax) gPad->WaitPrimitive();
+// canvas->cd(3);
+// gPad->Clear();
+ }
+}
//__________________________________________________________________________________________________
void AliRICHDisplFast::Exec(Option_t *)
{
144,0,AliRICHParam::PcSizeY());
if(isClusters) pClustersH2 = new TH2F("pClustersH2","Event Display",165,0,AliRICHParam::PcSizeX(),
144,0,AliRICHParam::PcSizeY());
-
-
-
for(Int_t iEventN=0;iEventN<gAlice->GetEventsPerRun();iEventN++){//events Loop
pRich->GetLoader()->GetRunLoader()->GetEvent(iEventN);
#include <TTask.h>
#include "AliRICH.h"
#include <AliRun.h>
-
+#include "TH2.h"
class AliRICH;
class AliRICHDisplFast : public TTask
public :
AliRICHDisplFast() {;}
virtual ~AliRICHDisplFast() {;}
- static void DrawSectors(); //Draw sectors in plot
- virtual void Exec(Option_t *opt=0); //virtual do the main job
+ static void DrawSectors(); //Draw sectors in plot
+ void ShowEvent(Int_t iEvtNmin,Int_t iEvtNmax); // Display looping on events
+ void ShowEvent(Int_t iEvent) {ShowEvent(iEvent,iEvent);} // Display only one event
+ virtual void Exec(Option_t *opt=0); //virtual do the main job
protected:
- ClassDef(AliRICHDisplFast,0) //Utility class to draw the current event topology
+ ClassDef(AliRICHDisplFast,0) //Utility class to draw the current event topology
};
#endif // #ifdef AliRICHDisplFast_cxx
if(fpClusters->GetEntries()==0) return kBad;//no clusters at all for a given track
Bool_t kPatRec = kFALSE;
+
+ AliDebug(1,Form("---Track Parameters--- Theta: %f , Phi: %f ",GetTrackTheta()*TMath::RadToDeg(),GetTrackPhi()*TMath::RadToDeg()));
Int_t candidatePhotons = 0;
SetPhotonFlag(1);
FindThetaPhotonCerenkov();
Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov();
- AliDebug(1,Form("THETA CERENKOV ---> %i",thetaPhotonCerenkov));
+ AliDebug(1,Form("THETA CERENKOV ---> %f",thetaPhotonCerenkov));
SetPhotonEta(thetaPhotonCerenkov);
candidatePhotons++;
}//clusters loop
Int_t nPhotonHough = 0;
Float_t thetaCerenkov = GetThetaCerenkov();
- AliDebug(1,Form(" fThetaCerenkov ",thetaCerenkov));
+ AliDebug(1,Form(" fThetaCerenkov %f ",thetaCerenkov));
Float_t thetaDist= thetaCerenkov - fThetaMin;
Int_t steps = (Int_t)(thetaDist / fDTheta);
//invoked by AliRecontruction for RICH
//if ESD doesn't contain tracks, try to reconstruct with particle from STACK
//(this case is just to forsee a standalone RICH simulation
-
+ TNtupleD *hn=0;
AliDebug(1,"Start pattern recognition");
if(pESD->GetNumberOfTracks()) RecWithESD(pESD);
else
- RecWithStack();
+ RecWithStack(hn);
AliDebug(1,"Stop pattern recognition");
return 0; // error code: 0=no error;
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
+ pTrack->GetPxPyPz(pb);
+ Int_t status=pTrack->GetStatus()&AliESDtrack::kTOFout;//get running track parameters
AliDebug(1,Form("Track %i pmod=%f mass=%f stat=%i",iTrackN,pTrack->GetP(),pTrack->GetMass(),status));
x0.SetXYZ(xb[0],xb[1],xb[2]); p0.SetXYZ(xb[0],xb[1],xb[2]);
AliRICHHelix helix(x0,p0,pTrack->GetSign(),b);
AliDebug(1,Form("FINAL Theta Cerenkov=%f",thetaCerenkov));
pTrack->SetRICHsignal(thetaCerenkov);
- Double_t richPID[5]={0.2,0.2,0.2,0.2,0.2}; //start with equal probs for (e,mu,pi,k,p)
- CalcProb(thetaCerenkov,p0.Mag(),richPID);
- pTrack->SetRICHpid(richPID);
+// Double_t richPID[5]={0.2,0.2,0.2,0.2,0.2}; //start with equal probs for (e,mu,pi,k,p)
+// CalcProb(thetaCerenkov,p0.Mag(),richPID);
+// pTrack->SetRICHpid(richPID);
}//ESD tracks loop
AliDebug(1,"Stop.");
} //RecWithESD
//__________________________________________________________________________________________________
-void AliRICHTracker::RecWithStack()
+void AliRICHTracker::RecWithStack(TNtupleD *hn)
{
// reconstruction for particles from STACK
//
AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH"));
- pRich->GetLoader()->GetRunLoader()->LoadHeader();pRich->GetLoader()->GetRunLoader()->LoadKinematics();
+// pRich->GetLoader()->GetRunLoader()->LoadHeader();
+ pRich->GetLoader()->GetRunLoader()->LoadKinematics();
AliStack *pStack = pRich->GetLoader()->GetRunLoader()->Stack();
Int_t iNtracks=pStack->GetNtrack();
-
+ AliDebug(1,Form(" Start reconstruction with %i track(s) from Stack",iNtracks));
+
+ Double_t hnvec[12];
+
Double_t b=GetFieldMap()->SolenoidField()/10;// magnetic field in Tesla
AliDebug(1,Form("Start with simulated %i tracks in %f Tesla field",iNtracks,b));
TVector3 x0(0,0,0); TVector3 p0(0,0,0);//tmp storage for AliRICHHelix
+ if(pRich->GetLoader()->LoadRecPoints()) {AliDebug(1,Form("No clusters found in RICH"));return;}
+ pRich->GetLoader()->TreeR()->GetEntry(0);
+
for(Int_t iTrackN=0;iTrackN<iNtracks;iTrackN++){//ESD tracks loop
TParticle *pParticle = pStack->Particle(iTrackN);
+ if(pParticle->GetMother(0)!=-1) continue; //consider only primaries
+ if(pParticle->GetPDG()->Charge()==0) continue; //to avoid photons from stack...
+ hnvec[0]=pParticle->P();
+ hnvec[1]=pParticle->GetPDG()->Charge();
+ hnvec[2]=pParticle->Theta();
+ hnvec[3]=pParticle->Phi();
p0.SetMagThetaPhi(pParticle->P(),pParticle->Theta(),pParticle->Phi());
x0.SetXYZ(pParticle->Vx(),pParticle->Vy(),pParticle->Vz());
AliRICHHelix helix(x0,p0,Int_t(pParticle->GetPDG()->Charge()/3+0.1),b);
Int_t iChamber=helix.RichIntersect(pRich->P());
AliDebug(1,Form("intersection with %i chamber found",iChamber));
- if(!iChamber) continue;//intersection with no chamber found
-
- Double_t distMip=9999; //min distance between clusters and track position on PC
- Int_t iMipId=0; //index of that min distance cluster
+ if(!iChamber) continue;// no intersection with RICH found
+ hnvec[4]=helix.PosPc().X();
+ hnvec[5]=helix.PosPc().Y();
+ Double_t distMip=9999; //min distance between clusters and track position on PC
+ Double_t mipX=kBad; //min distance between clusters and track position on PC
+ Double_t mipY=kBad; //min distance between clusters and track position on PC
+ Double_t chargeMip=kBad; // charge MIP to find
+ Int_t iMipId=kBad; //index of that min distance cluster
for(Int_t iClusN=0;iClusN<pRich->Clusters(iChamber)->GetEntries();iClusN++){//clusters loop for intersected chamber
AliRICHcluster *pClus=(AliRICHcluster*)pRich->Clusters(iChamber)->UncheckedAt(iClusN);//get pointer to current cluster
Double_t distCurrent=pClus->DistTo(helix.PosPc());//ditance between current cluster and helix intersection with PC
- if(distCurrent<distMip){distMip=distCurrent;iMipId=iClusN;}//find cluster nearest to the track
+ if(distCurrent<distMip){distMip=distCurrent;mipX=pClus->X();
+ mipY=pClus->Y();
+ chargeMip=pClus->Q();iMipId=1000000*iChamber+iClusN;}//find cluster nearest to the track
AliDebug(1,Form("Ploc (%f,%f,%f) dist= %f",helix.Ploc().Mag(),helix.Ploc().Theta()*TMath::RadToDeg(),
helix.Ploc().Phi()*TMath::RadToDeg(),pClus->DistTo(helix.PosPc())));
}////clusters loop for intersected chamber
AliDebug(1,Form("Min distance cluster: %i dist is %f",iMipId,distMip));
-
+ hnvec[6]=mipX;hnvec[7]=mipY;
+ hnvec[8]=chargeMip;
AliRICHRecon recon(&helix,pRich->Clusters(iChamber),iMipId);
Double_t thetaCerenkov=recon.ThetaCerenkov(); //search for mean Cerenkov angle for this track
+ hnvec[9]=thetaCerenkov;
+ hnvec[10]=recon.GetHoughPhotons();
+ hnvec[11]=(Double_t)iMipId;
+ if(hn) hn->Fill(hnvec);
AliDebug(1,Form("FINAL Theta Cerenkov=%f",thetaCerenkov));
// pTrack->SetRICHsignal(thetaCerenkov);
- Double_t richPID[5]={0.2,0.2,0.2,0.2,0.2}; //start with equal probs for (e,mu,pi,k,p)
- CalcProb(thetaCerenkov,p0.Mag(),richPID);
+// Double_t richPID[5]={0.2,0.2,0.2,0.2,0.2}; //start with equal probs for (e,mu,pi,k,p)
+// CalcProb(thetaCerenkov,p0.Mag(),richPID);
}//ESD tracks loop
+
+ pRich->GetLoader()->UnloadRecPoints();
+
AliDebug(1,"Stop.");
} //RecWithStack
#include <AliTracker.h>
#include <AliLog.h>
+#include "TNtupleD.h"
class AliCluster;
class AliESD;
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 RecWithStack(); //recon from Stack in case ESD empty
+ void RecWithStack(TNtupleD *hn); //recon from Stack in case ESD empty
void CalcProb(Double_t thetaCer,Double_t pmod,Double_t *richPID); // calculate pid for RICH
Int_t LoadClusters(TTree *); //pure virtual from AliTracker
distH1->Draw();
}//ParticleContribs()
+//__________________________________________________________________________________________________
+void CheckPR()
+{
+//Pattern recognition wirh Stack particles
+ TFile *pFile = new TFile("$(HOME)/RPR.root","RECREATE","RICH Pattern Recognition");
+ TNtupleD *hn = new TNtupleD("hn","ntuple","Pmod:Charge:TrackTheta:TrackPhi:TrackX:TrackY:MinX:MinY:ChargeMIP:ThetaCerenkov:NPhotons:MipIndex");
+ for(Int_t iEvtN=0;iEvtN<R()->GetLoader()->GetRunLoader()->GetNumberOfEvents();iEvtN++) {
+ R()->GetLoader()->GetRunLoader()->GetEvent(iEvtN);
+ AliRICHTracker *tr = new AliRICHTracker();
+ tr->RecWithStack(hn);
+ Info("CheckPR","Pattern Recognition done for event %i",iEvtN);
+ }
+ pFile->Write();pFile->Close();
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff