Possibility to run standalone PatRec without run again ESD and creation of ntuples...
authordibari <dibari@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 23 Aug 2005 09:56:04 +0000 (09:56 +0000)
committerdibari <dibari@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 23 Aug 2005 09:56:04 +0000 (09:56 +0000)
RICH/AliRICH.cxx
RICH/AliRICH.h
RICH/AliRICHRecon.cxx
RICH/AliRICHRecon.h
RICH/AliRICHTracker.cxx
RICH/AliRICHTracker.h
RICH/RichMenu.C

index a73a87a..b9981d2 100644 (file)
@@ -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;iEvtN<iEvtNmax;iEvtN++){
-    if(iEvtN%10==0) AliInfo(Form("events processed %i",iEvtN));
+  for(Int_t iEvtN=iEvtNmin;iEvtN<iEvtNmax;iEvtN++){    
+    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};
+    AliInfo(Form("events processed %i",iEvtN));
     if(GetLoader()->GetRunLoader()->GetEvent(iEvtN)) return;    
     AliStack *pStack = GetLoader()->GetRunLoader()->Stack();
     for(Int_t iPrimN=0;iPrimN<GetLoader()->TreeH()->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;iEvtN<GetLoader()->GetRunLoader()->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;iTrackN<iNtracks;iTrackN++){//ESD tracks loop
+        if(iTrackN%100==0)Info("RichAna",Form("Track %i to be processed",iTrackN));
+        AliRICHTracker *pTrRich = new AliRICHTracker();
+        pTrRich->RecWithESD(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"    );}
index 10cc584..be8a3f3 100644 (file)
@@ -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:  
index b3021fd..4466b3c 100644 (file)
@@ -25,6 +25,7 @@
 #include <TMath.h>
 #include <TRotation.h>
 #include <TVector3.h>
+#include <TH1F.h>
 
 #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<bin2; j++) 
-         {
-           hcs[j] += 1; 
-           hcsw[j] += weight;
-         }
-      }
-  }
-  
-
-  if(weightFlag == 0) 
-    {
-      hCSspace = hcsw;
-    }
-  else
-    {
-      hCSspace = hcs;
-      //      cout << " probems with weight...normal procedure adopted " << endl;
-    }
-
-  HoughFiltering(hCSspace);
-
-  for (bin=0; bin <nBin; bin++) {
-    angle = (bin+0.5) * (fDTheta);
-    if (hCSspace[bin] && hCSspace[bin] > etaPeakPos) {
-      etaPeakCount = 0;
-      etaPeakPos = hCSspace[bin];
-      etaPeak[0]=angle;
-    }
-    else { 
-      if (hCSspace[bin] == etaPeakPos) {
-       etaPeak[++etaPeakCount] = angle;
-      }
-    }
-  } 
-
-  for (i=0; i<etaPeakCount+1; i++) {
-    thetaCerenkov += etaPeak[i];
-  }
-  if (etaPeakCount>=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 && nxDx<nBin)
-             hcsFilt[nx] +=  hcs[nxDx] * k[i];
-      }      
-   }
-     
-   for (nx = 0; nx < nBin; nx++) {
-     hcs[nx] = hcsFilt[nx];
-   }
-}
+    phots->Fill(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(sumPhots<fMinNumPhots) continue; // cut on minimum n. of photons per ring
+    Double_t sumPhotsw=photsw->Integral(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);
 }
-
index 0c520b5..e1d383b 100644 (file)
@@ -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
index 8d03a9b..6ee209e 100644 (file)
@@ -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;iTrackN<iNtracks;iTrackN++){//ESD tracks loop
+      RecWithESD(pESD,pRich,iTrackN);
+    }
+  }
+  else RecWithStack(0);
   AliDebug(1,"Stop pattern recognition");
   return 0; // error code: 0=no error;
 }//PropagateBack()
 //__________________________________________________________________________________________________
-void AliRICHTracker::RecWithESD(AliESD *pESD)
+void AliRICHTracker::RecWithESD(AliESD *pESD,AliRICH *pRich,Int_t iTrackN)
 {
 //recontruction from ESD- primary way to reconstruct particle ID signal from tracks provided by core detectors
 
-  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"));
-  
-  for(Int_t iTrackN=0;iTrackN<iNtracks;iTrackN++){//ESD tracks loop
+    Double_t fField=GetFieldMap()->SolenoidField()/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;iPart<AliPID::kSPECIES;iPart++) {
         sigmaPID[iPart] = 0;
-       richPID[iPart] = 0;
         for(Int_t iphot=0;iphot<pRich->Clusters(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
 //__________________________________________________________________________________________________
index 52ae65f..ec6e99a 100644 (file)
@@ -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
 
index f01946c..66a239c 100644 (file)
@@ -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()
 //__________________________________________________________________________________________________