Minor changes ot protect against 1/0 values.
authordibari <dibari@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 8 Nov 2005 16:42:52 +0000 (16:42 +0000)
committerdibari <dibari@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 8 Nov 2005 16:42:52 +0000 (16:42 +0000)
RICH/AliRICHParam.cxx
RICH/AliRICHRecon.cxx
RICH/AliRICHReconstructor.cxx
RICH/AliRICHReconstructor.h
RICH/AliRICHTracker.cxx
RICH/RichMenu.C

index 3444026..bc62a72 100644 (file)
@@ -257,7 +257,7 @@ void AliRICHParam::ReadErrFiles()
   FILE *pChromErr, *pGeomErr, *pLocErr;  
 
   if(!count) {
-     AliInfoClass("reading RICH error parameters...");
+     AliInfoGeneral("ReadErrFiles","reading RICH error parameters...");
      pChromErr = fopen(Form("%s/RICH/RICHConfig/SigmaChromErr.txt",gSystem->Getenv("ALICE_ROOT")),"r");
      pGeomErr  = fopen(Form("%s/RICH/RICHConfig/SigmaGeomErr.txt",gSystem->Getenv("ALICE_ROOT")),"r");
      pLocErr   = fopen(Form("%s/RICH/RICHConfig/SigmaLocErr.txt",gSystem->Getenv("ALICE_ROOT")),"r");
@@ -279,7 +279,7 @@ void AliRICHParam::ReadErrFiles()
        fgErrLoc[2][i] = l2;
        fgErrLoc[3][i] = l3;    
      }
-     AliInfoClass("DONE successfully!");
+     AliInfoGeneral("ReadErrFiles","DONE successfully!");
      fclose(pChromErr);
      fclose(pGeomErr);
      fclose(pLocErr);
index 7c4cfcd..f73e178 100644 (file)
@@ -68,14 +68,14 @@ Double_t AliRICHRecon::ThetaCerenkov()
 // Pattern recognition method based on Hough transform
 // Return theta Cerenkov for a given track and list of clusters which are set in ctor  
 
-  if(fpClusters->GetEntries()==0) return -1;//no clusters at all for a given track
+  if(fpClusters->GetEntries()==0) return -10;//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;
 
-  SetThetaCerenkov(999.);
+  SetThetaCerenkov(-1);
   SetHoughPhotons(0);
   SetHoughPhotonsNorm(0);
 
@@ -116,8 +116,8 @@ Double_t AliRICHRecon::ThetaCerenkov()
  
   if(nPhotonHough < 1) 
     {
-      SetThetaCerenkov(999.);
-      SetHoughPhotonsNorm(0.);
+      SetThetaCerenkov(-1);
+      SetHoughPhotonsNorm(0);
       return -1;
     }
 
index d1b54a8..f5198be 100644 (file)
 #include <AliRawReader.h>         //Reconstruct(...)
 #include <AliRun.h>               //ConvertDigits uses gAlice
 #include <AliESD.h>            //RichAna()
+#include <AliStack.h>          //RichAna()
 #include <TFile.h>             //RichAna()
 #include <TMinuit.h>              //Dig2Clu()
+#include <TParticle.h>            //TParticle()
 ClassImp(AliRICHReconstructor)
 
 void AliRICHReconstructor::Reconstruct(AliRunLoader *pAL,AliRawReader* pRR)const
@@ -109,28 +111,63 @@ void  AliRICHReconstructor::FormCluster(AliRICHCluster *pClu,AliRICHDigit *pDig,
     if((pDig=UseDig(x[i],y[i],pDigList,pDigMap))) FormCluster(pClu,pDig,pDigList,pDigMap);    
 }//FormCluster()
 //__________________________________________________________________________________________________
-void AliRICHReconstructor::RichAna(Int_t nNevMax,Bool_t askPatRec)
+void AliRICHReconstructor::CheckPR()
+{
+//Pattern recognition with 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:Chamber:Particle");
+//  printf("\n\n");
+//  printf("Pattern Recognition done for event %5i",0);
+  AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH"));
+  AliMagF * magf = gAlice->Field();
+  AliTracker::SetFieldMap(magf,kTRUE);
+  for(Int_t iEvtN=0;iEvtN<pRich->GetLoader()->GetRunLoader()->GetNumberOfEvents();iEvtN++) {
+    pRich->GetLoader()->GetRunLoader()->GetEvent(iEvtN);
+    AliRICHTracker *tr = new AliRICHTracker();
+    tr->RecWithStack(hn);
+    AliInfoClass(Form("Pattern Recognition done for event %i \b",iEvtN));
+//    printf("\b\b\b\b\b%5i",iEvtN+1);
+  }
+  printf("\n\n");
+  pFile->Write();pFile->Close();
+}
+//__________________________________________________________________________________________________
+void AliRICHReconstructor::RichAna(Int_t iNevMin,Int_t iNevMax,Bool_t askPatRec)
 {
   TFile *pFile=TFile::Open("AliESDs.root","read");
-  if(!pFile || !pFile->IsOpen()) {return;}      //open AliESDs.root                                                                    
+  if(!pFile || !pFile->IsOpen()) {AliInfoClass("ESD file not open.");return;}      //open AliESDs.root                                                                    
   TTree *pTree = (TTree*) pFile->Get("esdTree");
-  if(!pTree){return;}                               //get ESD tree  
-
+  if(!pTree){AliInfoClass("ESD not found.");return;}                               //get ESD tree  
+  AliInfoClass("ESD found. Go ahead!");
+  
   AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH"));
+
+  AliMagF * magf = gAlice->Field();
+  AliTracker::SetFieldMap(magf,kTRUE);
+  pRich->GetLoader()->GetRunLoader()->LoadHeader();
+  pRich->GetLoader()->GetRunLoader()->LoadKinematics();
+  TString var1 = "Pmod:Charge:TrackTheta:TrackPhi:MinX:MinY:ThetaCerenkov:NPhotons:";
+  TString var2 = "ChargeMIP:Chamber:TOF:LengthTOF:prob1:prob2:prob3:";
+  TString var3 = "ErrPar1:ErrPar2:ErrPar3:Th1:Th2:Th3:nPhotBKG:pdgCode";
+  TString varList = var1+var2+var3;
   
-  AliMagF * magf = gAlice->Field();  AliTracker::SetFieldMap(magf,kTRUE); //  AliTracker::SetFieldMap(magf);
+  Double_t dx,dy;
+  Double_t hnvec[30];
 
-  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:ErrPar1:ErrPar2:ErrPar3:Th1:Th2:Th3:nPhotBKG");
-  if(nNevMax==0) nNevMax=999999;
-  if(pRich->GetLoader()->GetRunLoader()->GetNumberOfEvents()<nNevMax) nNevMax = pRich->GetLoader()->GetRunLoader()->GetNumberOfEvents();
+//  TFile *pFileRA = new TFile("$(HOME)/RichAna.root","RECREATE","RICH Pattern Recognition");
+  TFile *pFileRA = new TFile("./RichAna.root","RECREATE","RICH Pattern Recognition");
+  TNtupleD *hn = new TNtupleD("hn","ntuple",varList);
+  if(iNevMin<0) iNevMin=0;
+  if(iNevMin>iNevMax) {iNevMin=0;iNevMax=0;}  
+  if(iNevMax==0) iNevMax=999999;
+  if(pRich->GetLoader()->GetRunLoader()->GetNumberOfEvents()<iNevMax) iNevMax = pRich->GetLoader()->GetRunLoader()->GetNumberOfEvents();
   AliESD *pESD=new AliESD;  pTree->SetBranchAddress("ESD", &pESD);
-  for(Int_t iEvtN=0;iEvtN<nNevMax;iEvtN++) {
+  for(Int_t iEvtN=iNevMin;iEvtN<iNevMax;iEvtN++) {
     pTree->GetEvent(iEvtN);
     pRich->GetLoader()->GetRunLoader()->GetEvent(iEvtN);
     pRich->GetLoader()->LoadRecPoints();
     pRich->GetLoader()->TreeR()->GetEntry(0);
+    AliStack *pStack = pRich->GetLoader()->GetRunLoader()->Stack();
 //Pattern recognition started
     if(pESD->GetNumberOfTracks()) {
       Int_t iNtracks=pESD->GetNumberOfTracks();
@@ -140,11 +177,15 @@ void AliRICHReconstructor::RichAna(Int_t nNevMax,Bool_t askPatRec)
         AliRICHTracker *pTrRich = new AliRICHTracker();
         if(askPatRec==kTRUE) pTrRich->RecWithESD(pESD,pRich,iTrackN);
         AliESDtrack *pTrack = pESD->GetTrack(iTrackN);// get next reconstructed track
-        Double_t dx,dy;
-        Double_t hnvec[30];
+        
+        Int_t lab=TMath::Abs(pTrack->GetLabel());
+        TParticle *pPart=pStack->Particle(lab);
+        Int_t code=pPart->GetPdgCode();
+
         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]);
@@ -170,7 +211,8 @@ void AliRICHReconstructor::RichAna(Int_t nNevMax,Bool_t askPatRec)
           if(cosThetaTh>=1) continue;
           hnvec[18+i]= TMath::ACos(cosThetaTh);
         }
-        hnvec[21]=pTrRich->fnPhotBKG;
+        if(askPatRec==kTRUE) hnvec[21]=pTrRich->fnPhotBKG; else hnvec[21]=0;
+        hnvec[22]=code;
         hn->Fill(hnvec);
       }
     }
@@ -178,21 +220,5 @@ void AliRICHReconstructor::RichAna(Int_t nNevMax,Bool_t askPatRec)
   }
   pFileRA->Write();pFileRA->Close();// close RichAna.root
   delete pESD;  pFile->Close();//close AliESDs.root
-}//RichAna()
-//__________________________________________________________________________________________________
-void AliRICHReconstructor::CheckPR()
-{
-//Pattern recognition with 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:Chamber:Particle");
-  AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH"));
-  AliMagF * magf = gAlice->Field();
-  AliTracker::SetFieldMap(magf,kTRUE);
-  for(Int_t iEvtN=0;iEvtN<pRich->GetLoader()->GetRunLoader()->GetNumberOfEvents();iEvtN++) {
-    pRich->GetLoader()->GetRunLoader()->GetEvent(iEvtN);
-    AliRICHTracker *tr = new AliRICHTracker();
-    tr->RecWithStack(hn);
-    AliInfoClass(Form("Pattern Recognition done for event %i \b",iEvtN));
-  }
-  pFile->Write();pFile->Close();
 }
+//__________________________________________________________________________________________________
index 135989a..608669b 100644 (file)
@@ -25,7 +25,7 @@ public:
          void          FormCluster(AliRICHCluster *pClu,AliRICHDigit *pDig,TClonesArray *pDigList,TMatrixF *pDigMap)const;//form cluster recursive algorithm
   inline AliRICHDigit *UseDig     (Int_t padX,Int_t padY,TClonesArray *pDigList,TMatrixF *pDigMap                  )const;//use this pad's digit to form a cluster
   static void          CheckPR    (                                                                                );     //utility-> run staff for stack
-  static void          RichAna    (Int_t iNevMax=99999,Bool_t isPatRec=kFALSE                                      );     //utility-> create ntuples for analysis
+  static void          RichAna    (Int_t iNevMin=0, Int_t iNevMax=99999,Bool_t isPatRec=kFALSE                                      );     //utility-> create ntuples for analysis
   
 protected:
   ClassDef(AliRICHReconstructor, 0)   //class for the RICH reconstruction
index 345268a..7bc97e7 100644 (file)
@@ -48,15 +48,21 @@ void AliRICHTracker::RecWithESD(AliESD *pESD,AliRICH *pRich,Int_t iTrackN)
     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) return;//intersection with no chamber found
-//find MIP cluster candidate (cluster which is closest to track intersection point)    
+    if(!iChamber) {
+      pTrack->SetRICHsignal(-999); //to be improved by flags...
+      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
     Double_t chargeMip=0; //charge of the MIP
+    Bool_t kFound = kFALSE;
     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
+      if(pClus->Q()<AliRICHParam::QthMIP()) continue;
       Double_t distCurrent=pClus->DistTo(helix.PosPc());//distance between current cluster and helix intersection with PC
       if(distCurrent<distMip){
+        kFound = kTRUE;
         distMip=distCurrent;
         iMipId=iClusN;
         distX=pClus->DistX(helix.PosPc());
@@ -66,14 +72,21 @@ void AliRICHTracker::RecWithESD(AliESD *pESD,AliRICH *pRich,Int_t iTrackN)
       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
-    
+
+    if(!kFound) {
+      pTrack->SetRICHsignal(-999);
+      return;
+    }
     AliDebug(1,Form("Min distance cluster: %i dist is %f",iMipId,distMip));
+    pTrack->SetRICHcluster(((Int_t)chargeMip)+1000000*iChamber);
+    pTrack->SetRICHdxdy(distX,distY);
+    pTrack->SetRICHthetaPhi(helix.Ploc().Theta(),helix.Ploc().Phi());
 //
 // HERE CUTS ON GOLD RINGS....
 //
-    if(distMip>AliRICHParam::DmatchMIP()||chargeMip<AliRICHParam::QthMIP()) {
+    if(distMip>AliRICHParam::DmatchMIP()) {
       //track not accepted for pattern recognition
-      pTrack->SetRICHsignal(-999.); //to be improved by flags...
+      pTrack->SetRICHsignal(-990); //to be improved by flags...
       return;
     }
 //
@@ -81,9 +94,6 @@ void AliRICHTracker::RecWithESD(AliESD *pESD,AliRICH *pRich,Int_t iTrackN)
 
     Double_t thetaCerenkov=recon.ThetaCerenkov(); //search for mean Cerenkov angle for this track
     
-    pTrack->SetRICHcluster(((Int_t)chargeMip)+1000000*iChamber);
-    pTrack->SetRICHdxdy(distX,distY);
-    pTrack->SetRICHthetaPhi(helix.Ploc().Theta(),helix.Ploc().Phi());
     pTrack->SetRICHsignal(thetaCerenkov);
     pTrack->SetRICHnclusters(recon.GetHoughPhotons());
 
@@ -103,17 +113,15 @@ void AliRICHTracker::RecWithESD(AliESD *pESD,AliRICH *pRich,Int_t iTrackN)
           if(recon.GetPhotonFlag() == 2) {
             Double_t theta_g=recon.GetTrackTheta();
             Double_t phi_g=(recon.GetPhiPoint()-recon.GetTrackPhi());
-            Double_t sigma2 = AliRICHParam::SigmaSinglePhoton(iPart,pTrack->GetP(),theta_g,phi_g).Mag2();
-           if (sigma2>0)
-             sigmaPID[iPart] += 1/sigma2;
+            Double_t sigma2 = AliRICHParam::SigmaSinglePhoton(iPart,pTrack->GetP(),theta_g,phi_g).Mag2(); 
+            if(sigma2>0) sigmaPID[iPart] += 1/sigma2;
           }
         }
-       sigmaPID[iPart] *= (Double_t)(recon.GetHoughPhotons()-fnPhotBKG)/(Double_t)(recon.GetHoughPhotons()); // n total phots, m are background...the sigma are scaled..
        if (sigmaPID[iPart]>0)
-         sigmaPID[iPart] = 1/TMath::Sqrt(sigmaPID[iPart])*0.001; // sigma from parametrization are in mrad...
-       else
-         sigmaPID[iPart] = 0;
-       fErrPar[iPart]=sigmaPID[iPart];
+          sigmaPID[iPart] *= (Double_t)(recon.GetHoughPhotons()-fnPhotBKG)/(Double_t)(recon.GetHoughPhotons()); // n total phots, m are background...the sigma are scaled..
+         if(sigmaPID[iPart]>0) sigmaPID[iPart] = 1/TMath::Sqrt(sigmaPID[iPart])*0.001; // sigma from parametrization are in mrad...
+          else                  sigmaPID[iPart] = 0;
+          fErrPar[iPart]=sigmaPID[iPart];
         AliDebug(1,Form("sigma for %s is %f rad",AliPID::ParticleName(iPart),sigmaPID[iPart]));
       }
       CalcProb(thetaCerenkov,pTrack->GetP(),sigmaPID,richPID);
@@ -240,7 +248,8 @@ void AliRICHTracker::CalcProb(Double_t thetaCer,Double_t pmod, Double_t *sigmaPI
 //    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);
-    height[iPart] = TMath::Gaus(thetaCer,thetaTh[iPart],sigmaPID[iPart],kTRUE);
+    if(sigmaPID[iPart]>0) height[iPart] = TMath::Gaus(thetaCer,thetaTh[iPart],sigmaPID[iPart],kTRUE);
+    else                  height[iPart] = 0;
     totalHeight +=height[iPart];
     AliDebug(1,Form(" Particle %s with mass %f with height %f and thetaTH %f",AliPID::ParticleName(iPart),mass,height[iPart],thetaTh[iPart]));
     AliDebug(1,Form(" partial height %15.14f total height %15.14f",height[iPart],totalHeight));
@@ -248,7 +257,7 @@ void AliRICHTracker::CalcProb(Double_t thetaCer,Double_t pmod, Double_t *sigmaPI
   if(totalHeight<1e-5) {for(Int_t iPart=0;iPart<AliPID::kSPECIES;iPart++)richPID[iPart]=1.0/AliPID::kSPECIES;return;}
   for(Int_t iPart=0;iPart<AliPID::kSPECIES;iPart++) richPID[iPart] = height[iPart]/totalHeight;
   Int_t iPartNear = TMath::LocMax(AliPID::kSPECIES,richPID);
-  if(TMath::Abs(thetaCer-thetaTh[iPartNear]) > 5*sigmaPID[iPartNear]) for(Int_t iPart=0;iPart<AliPID::kSPECIES;iPart++)richPID[iPart]=1.0/AliPID::kSPECIES;
+  if(TMath::Abs(thetaCer-thetaTh[iPartNear])>5.*sigmaPID[iPartNear]) for(Int_t iPart=0;iPart<AliPID::kSPECIES;iPart++)richPID[iPart]=1.0/AliPID::kSPECIES;
   //last line is to check if the nearest thetacerenkov to the teorethical one is within 5 sigma, otherwise no response (equal prob to every particle
 
 }//CalcProb
index 03a6cfe..e40092d 100644 (file)
@@ -107,16 +107,16 @@ void RichMenu()
   if(AliceRead()){//it's from file, show some info
     pMenu->AddButton("Display single chambers"         ,"r->Display();"  , "Display Fast");
     pMenu->AddButton("Display ALL chambers"            ,"r->DisplayEvent(0,0);"  , "Display Fast");
-    pMenu->AddButton("Hits QA"                         ,"hqa()"    ,"????");
+    pMenu->AddButton("Hits QA"                         ,"hqa()"    ,"QA plots for hits: hqa()");
     pMenu->AddButton("Recon with stack"                ,"AliRICHReconstructor::CheckPR(        )","Create RSR.root with ntuple hn");    
-    pMenu->AddButton("RichAna no Recon"                ,"AliRICHReconstructor::RichAna(0,kFALSE)","Create RichAna.root with ntuple hn without PatRec");    
-    pMenu->AddButton("RichAna with Recon"              ,"AliRICHReconstructor::RichAna(0,kTRUE )","Create RichAna.root with ntuple hn with PatRec");    
-    pMenu->AddButton("Print hits"                      ,"h();"      ,"????");
-    pMenu->AddButton("Print sdigits"                   ,"s();"   ,"????");
-    pMenu->AddButton("Print digits"                    ,"d();"    ,"????");
-    pMenu->AddButton("Print clusters"                  ,"c();"  ,"????");  
-    pMenu->AddButton("Print occupancy"                 ,"r->OccupancyPrint(-1);" ,"????");  
-    pMenu->AddButton("Print event summary  "           ,"r->SummaryOfEvent();"   ,"????");  
+    pMenu->AddButton("RichAna no Recon"                ,"AliRICHReconstructor::RichAna(0,0,kFALSE)","Create RichAna.root with ntuple hn without PatRec");    
+    pMenu->AddButton("RichAna with Recon"              ,"AliRICHReconstructor::RichAna(0,0,kTRUE )","Create RichAna.root with ntuple hn with PatRec");    
+    pMenu->AddButton("Print hits"                      ,"h();"      ,"To print hits: h()");
+    pMenu->AddButton("Print sdigits"                   ,"s();"      ,"To print sdigits: s()");
+    pMenu->AddButton("Print digits"                    ,"d();"      ,"To print digits: d()");
+    pMenu->AddButton("Print clusters"                  ,"c();"      ,"To print clusters: c()");  
+    pMenu->AddButton("Print occupancy"                 ,"r->OccupancyPrint(-1);" ,"To print occupancy");  
+    pMenu->AddButton("Print event summary  "           ,"r->SummaryOfEvent();"   ,"To print a summary of the event");  
   }else{//it's aliroot, simulate
     pMenu->AddButton("Debug ON",     "DebugON();",   "Switch debug on-off");   
     pMenu->AddButton("Debug OFF",    "DebugOFF();",  "Switch debug on-off");