AliDebug(1,"Start cluster finder.");AliRICHClusterFinder clus(GetRICH(pAL)); clus.Exec();
}
//__________________________________________________________________________________________________
-void AliRICHReconstructor::FillESD(AliRunLoader* /*pAL*/, AliESD* pESD) const
+void AliRICHReconstructor::FillESD(AliRunLoader* /*pAL*/, AliESD* /*pESD*/) const
{
-//This methode fills AliESDtrack with information from RICH
- AliDebug(1,Form("Start with %i tracks",pESD->GetNumberOfTracks()));
-/* const Double_t masses[5]={0.000511,0.105658,0.139567,0.493677,0.93828};//electron,muon,pion,kaon,proton
- const Double_t refIndex = 1.29052;
-
- Double_t thetaTh[5];
- Double_t sinThetaThNorm;
- Double_t sigmaThetaTh[5];
- Double_t height[5];
- Double_t totalHeight=0;
- Double_t x[3],p[3]; //tmp storage for track parameters
-
- for(Int_t iTrackN=0;iTrackN<pESD->GetNumberOfTracks();iTrackN++){//ESD tracks loop
- AliESDtrack *pTrack = pESD->GetTrack(iTrackN);// get next reconstructed track
-// if((pTrack->GetStatus()&AliESDtrack::kTOFout)==0) continue; //ignore tracks not recontructed by TOF
- pTrack->GetXYZ(x); pTrack->GetPxPyPz(p); Double_t pmod=pTrack->GetP();//get running track parameters
- continue;
-// AliRICHHelix helix(x[0],x[1],x[2],p[0],p[1],p[2]); //construct helix from running track parameters
-
-// TVector rad(1,5); TVector pc(1,5);
-// helix.RichIntersect(GetRICH(pAL)->P(),rad,pc); //returns cross point of track with RICH PC in LRS
-
- for(Int_t iPart=4;iPart>=0;iPart--){
- Double_t cosThetaTh = TMath::Sqrt(masses[iPart]*masses[iPart]+pmod*pmod)/(refIndex*pmod);
- if(cosThetaTh>=1) {break;}
- thetaTh[iPart] = TMath::ACos(cosThetaTh);
- sinThetaThNorm = TMath::Sin(thetaTh[iPart])/TMath::Sqrt(1-1/(refIndex*refIndex));
- sigmaThetaTh[iPart] = (0.014*(1/sinThetaThNorm-1) + 0.0043)*1.25;
- height[iPart] = TMath::Gaus(thetaExp,thetaTh[iPart],sigmaThetaTh[iPart]);
- totalHeight +=height[iPart];
- }
-
- Double_t richPID[5];
- for(Int_t iPart=0;iPart<5;iPart++) richPID[iPart] = height[iPart]/totalHeight;
- pTrack->SetRICHpid(richPID);
- }//ESD tracks loop
- */
+//All is done in AliRICHTracker. Nothing to do here.
}//FillESD
//__________________________________________________________________________________________________
AliRICH* AliRICHReconstructor::GetRICH(AliRunLoader* pAL) const
#include "AliRICHRecon.h"
#include <AliStack.h>
#include <TParticle.h>
-
+#include <TDatabasePDG.h>
+#include <TMath.h>
ClassImp(AliRICHTracker)
-
-
+//__________________________________________________________________________________________________
Int_t AliRICHTracker::PropagateBack(AliESD *pESD)
{
+ //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
+
+ AliDebug(1,"Start pattern recognition");
+ if(pESD->GetNumberOfTracks()) RecWithESD(pESD);
+ else
+ RecWithStack();
+ AliDebug(1,"Stop pattern recognition");
+
+ return 0; // error code: 0=no error;
+} //pure virtual from AliTracker
+//__________________________________________________________________________________________________
+void AliRICHTracker::RecWithESD(AliESD *pESD)
+{
+ //recontruction from ESD
+ //
Int_t iNtracks=pESD->GetNumberOfTracks();
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"));
-// pRich->GetLoader()->GetRunLoader()->LoadHeader();pRich->GetLoader()->GetRunLoader()->LoadKinematics();
-// AliStack *pStack = pRich->GetLoader()->GetRunLoader()->Stack();
-// TParticle *pParticle = pStack->Particle(0);
-// p0.SetMagThetaPhi(pParticle->P(),pParticle->Theta(),pParticle->Phi());
-
for(Int_t iTrackN=0;iTrackN<iNtracks;iTrackN++){//ESD tracks loop
AliESDtrack *pTrack = pESD->GetTrack(iTrackN);// get next reconstructed track
// if((pTrack->GetStatus()&AliESDtrack::kTOFout)==0) continue; //ignore tracks not recontructed by TOF
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));
-// x.Print();p.Print();
x0.SetXYZ(xb[0],xb[1],xb[2]); p0.SetXYZ(xb[0],xb[1],xb[2]);
AliRICHHelix helix(x0,p0,pTrack->GetSign(),b);
Int_t iChamber=helix.RichIntersect(pRich->P());
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);
}//ESD tracks loop
AliDebug(1,"Stop.");
- return 0;
-} //pure virtual from AliTracker
+} //RecWithESD
+//__________________________________________________________________________________________________
+void AliRICHTracker::RecWithStack()
+{
+ // reconstruction for particles from STACK
+ //
+ AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH"));
+
+ pRich->GetLoader()->GetRunLoader()->LoadHeader();pRich->GetLoader()->GetRunLoader()->LoadKinematics();
+ AliStack *pStack = pRich->GetLoader()->GetRunLoader()->Stack();
+ Int_t iNtracks=pStack->GetNtrack();
+
+ 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
+
+
+ for(Int_t iTrackN=0;iTrackN<iNtracks;iTrackN++){//ESD tracks loop
+ TParticle *pParticle = pStack->Particle(iTrackN);
+ 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
+ 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
+
+ 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));
+
+ AliRICHRecon recon(&helix,pRich->Clusters(iChamber),iMipId);
+ Double_t thetaCerenkov=recon.ThetaCerenkov(); //search for mean Cerenkov angle for this track
+ 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);
+
+ }//ESD tracks loop
+ AliDebug(1,"Stop.");
+} //RecWithStack
Int_t AliRICHTracker::LoadClusters(TTree *pTree)
{
// Load clusters for RICH
AliDebug(1,"Start."); pTree->GetEntry(0); AliDebug(1,"Stop."); return 0;
}
+
+//__________________________________________________________________________________________________
+void AliRICHTracker::CalcProb(Double_t thetaCer,Double_t pmod, Double_t *richPID)
+{
+//
+ Double_t height[5];Double_t totalHeight=0;
+ Int_t code[5]={kElectron,kMuonPlus,kPiPlus,kKPlus,kProton};
+ TDatabasePDG *db = new TDatabasePDG();
+ for(Int_t iPart=0;iPart<4;iPart++){
+ Double_t mass = db->GetParticle(code[iPart])->Mass();
+ Double_t refIndex=AliRICHParam::IndOfRefC6F14(6.755);
+ Double_t cosThetaTh = TMath::Sqrt(mass*mass+pmod*pmod)/(refIndex*pmod);
+ if(cosThetaTh>=1) {break;}
+ Double_t thetaTh = TMath::ACos(cosThetaTh);
+ Double_t sinThetaThNorm = TMath::Sin(thetaTh)/TMath::Sqrt(1-1/(refIndex*refIndex));
+ Double_t sigmaThetaTh = (0.014*(1/sinThetaThNorm-1) + 0.0043)*1.25;
+ height[iPart] = TMath::Gaus(thetaCer,thetaTh,sigmaThetaTh);
+ totalHeight +=height[iPart];
+ }
+ for(Int_t iPart=0;iPart<5;iPart++) richPID[iPart] = height[iPart]/totalHeight;
+}//CalcProb
fprintf(fp," pGen->Init();\n");
break;
case kGun1:
- fprintf(fp," AliGenFixed *pGen=new AliGenFixed(1);\n");
+ fprintf(fp," AliGenBox *pGen=new AliGenBox(1);\n");
fprintf(fp," pGen->SetPart(%i); pGen->SetOrigin(0,0,0);\n",fGenPartIdCombo->GetSelected());
- fprintf(fp," pGen->SetMomentumRange(%3.1f,%3.1f); pGen->SetTheta(pRICH->C(%i)->ThetaD()-2); pGen->SetPhi(pRICH->C(%i)->PhiD()); \n",
- float(fGenMinMomCombo->GetSelected())/10, float(fGenMaxMomCombo->GetSelected())/10,
- fGenChamberCombo->GetSelected(),fGenChamberCombo->GetSelected());
+ fprintf(fp," pGen->SetMomentumRange(%3.1f,%3.1f); \n",float(fGenMinMomCombo->GetSelected())/10, float(fGenMaxMomCombo->GetSelected())/10);
+ fprintf(fp," pGen->SetThetaRange(pRICH->C(%i)->ThetaD()-3,pRICH->C(%i)->ThetaD()-1); \n",fGenChamberCombo->GetSelected(),fGenChamberCombo->GetSelected());
+ fprintf(fp," pGen->SetPhiRange(pRICH->C(%i)->PhiD()-1,pRICH->C(%i)->PhiD()+1); \n",fGenChamberCombo->GetSelected(),fGenChamberCombo->GetSelected());
fprintf(fp," pGen->Init();\n");
break;
case kGun7:
fprintf(fp," AliGenCocktail *pCocktail=new AliGenCocktail();\n");
fprintf(fp," for(int i=1;i<=7;i++){\n");
fprintf(fp," AliGenFixed *pFixed=new AliGenFixed(1);\n");
- fprintf(fp," pFixed->SetPart(%i); pFixed->SetMomentum(2.5+i*0.4); pFixed->SetOrigin(0,0,0);\n",fGenPartIdCombo->GetSelected());
+ fprintf(fp," pFixed->SetPart(%i); pFixed->SetMomentum(1.0+i*0.5); pFixed->SetOrigin(0,0,0);\n",fGenPartIdCombo->GetSelected());
fprintf(fp," pFixed->SetPhi(pRICH->C(i)->PhiD()); pFixed->SetTheta(pRICH->C(i)->ThetaD()-2);\n");
fprintf(fp," pCocktail->AddGenerator(pFixed,Form(\"Fixed %i\",i),1);\n }\n");
fprintf(fp," pCocktail->Init();\n");
fprintf(fp," TStopwatch sw;TDatime time;\n\n");
fprintf(fp," AliSimulation *pSim=new AliSimulation;\n");
- if(fGenTypeCombo->GetSelected()==kGun1)
- fprintf(fp," pSim->SetRegionOfInterest(kTRUE); pSim->SetMakeDigitsFromHits(\"ITS TPC TRD TOF\");\n");
+ if(fGenTypeCombo->GetSelected()==kGun1||fGenTypeCombo->GetSelected()==kGun7) {
+ fprintf(fp," pSim->SetRegionOfInterest(kTRUE);\n");
+ fprintf(fp," pSim->SetMakeSDigits(\"TOF RICH\");\n");
+ fprintf(fp," pSim->SetMakeDigitsFromHits(\"ITS TPC TRD\");\n");
+ }
fprintf(fp," pSim->Run(iNevents); delete pSim;\n\n");
fprintf(fp," AliReconstruction *pRec=new AliReconstruction;\n");
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff