Added functionalities in display the RICH event and pattern recognition. Minor bugs...
authordibari <dibari@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 20 Oct 2004 19:50:41 +0000 (19:50 +0000)
committerdibari <dibari@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 20 Oct 2004 19:50:41 +0000 (19:50 +0000)
RICH/AliRICH.cxx
RICH/AliRICH.h
RICH/AliRICHChamber.cxx
RICH/AliRICHClusterFinder.cxx
RICH/AliRICHDisplFast.cxx
RICH/AliRICHDisplFast.h
RICH/AliRICHRecon.cxx
RICH/AliRICHTracker.cxx
RICH/AliRICHTracker.h
RICH/menu.C

index 842bd21..c818cef 100644 (file)
@@ -134,12 +134,21 @@ void AliRICH::Hits2SDigits()
         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 
index 6ab61c6..bbfba2e 100644 (file)
@@ -114,6 +114,8 @@ public:
           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     
index c945d3a..ac5c42a 100644 (file)
@@ -29,7 +29,7 @@ AliRICHChamber::AliRICHChamber(Int_t iChamber):TNamed()
 //  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.
index a0ae4fd..4e1a7d9 100644 (file)
@@ -109,7 +109,7 @@ void AliRICHClusterFinder::FindClusterContribs(AliRICHcluster *pCluster)
 //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;}
index 54a253f..7137e0c 100644 (file)
 #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 *)
 {
@@ -50,9 +116,6 @@ 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);
index f46ccd6..78feb00 100644 (file)
@@ -7,7 +7,7 @@
 #include <TTask.h>
 #include "AliRICH.h"
 #include <AliRun.h>
-
+#include "TH2.h"
 class AliRICH;
 
 class AliRICHDisplFast : public TTask 
@@ -15,10 +15,12 @@ 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
index 5ceada6..ccce9a7 100644 (file)
@@ -65,6 +65,8 @@ Double_t AliRICHRecon::ThetaCerenkov()
 
   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;
 
@@ -88,7 +90,7 @@ Double_t AliRICHRecon::ThetaCerenkov()
     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
@@ -1009,7 +1011,7 @@ void AliRICHRecon::FlagPhotons()
   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);
index bffdd84..bdd3b74 100644 (file)
@@ -17,11 +17,11 @@ 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
-  
+  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;
@@ -43,8 +43,8 @@ void AliRICHTracker::RecWithESD(AliESD *pESD)
     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);   
@@ -70,60 +70,87 @@ void AliRICHTracker::RecWithESD(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);         
+//    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
 
index 8a1f98d..c0054a8 100644 (file)
@@ -3,6 +3,7 @@
 
 #include <AliTracker.h>
 #include <AliLog.h>
+#include "TNtupleD.h"
 
 class AliCluster;
 class AliESD;
@@ -19,7 +20,7 @@ public:
   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 
 
index 10e3d18..c677746 100644 (file)
@@ -512,3 +512,17 @@ void ParticleContribs()
     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();
+}