Forseen the possibility to have RICH reconstruction from Stack (useful for debug...
authordibari <dibari@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 18 Oct 2004 15:00:49 +0000 (15:00 +0000)
committerdibari <dibari@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 18 Oct 2004 15:00:49 +0000 (15:00 +0000)
RICH/AliRICHReconstructor.cxx
RICH/AliRICHTracker.cxx
RICH/AliRICHTracker.h
RICH/CreateConfig.C

index 47424a9..db2dcfc 100644 (file)
@@ -37,45 +37,9 @@ void AliRICHReconstructor::Reconstruct(AliRunLoader* pAL) const
   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
index 0208b58..bffdd84 100644 (file)
@@ -8,12 +8,29 @@
 #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));
@@ -22,11 +39,6 @@ Int_t AliRICHTracker::PropagateBack(AliESD *pESD)
   
   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
@@ -34,7 +46,6 @@ Int_t AliRICHTracker::PropagateBack(AliESD *pESD)
       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());        
@@ -59,14 +70,86 @@ Int_t AliRICHTracker::PropagateBack(AliESD *pESD)
     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
index 7fd064d..8a1f98d 100644 (file)
@@ -12,12 +12,15 @@ class AliRICHTracker : public AliTracker
 {
 public:
            AliRICHTracker() :AliTracker()      {AliDebug(1,"Start.");}
-  virtual ~AliRICHTracker()                    {AliDebug(1,"Start.");}
+  virtual ~AliRICHTracker()                    {AliDebug(1,"Stop.");}
   Int_t Clusters2Tracks(AliESD *)              {AliDebug(1,"Start.");return 0;} //pure virtual from AliTracker 
   Int_t RefitInward(AliESD *)                  {AliDebug(1,"Start.");return 0;} //pure virtual from AliTracker 
   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 RecWithStack();                                                          //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 
 
 protected:
index f59d737..7608e51 100644 (file)
@@ -224,18 +224,18 @@ void RichConfig::CreateConfigFile()
       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");
@@ -320,8 +320,11 @@ void RichConfig::CreateRichBatch()
   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");