Now rings back again...

diff --git a/HMPID/AliHMPIDParam.h b/HMPID/AliHMPIDParam.h
index 85c0f8c55e6c7ffea3f78fdd7018fdc16978c213..e3c9ba55eac302718daaab099ba916f52bee3b21 100644 (file)
--- a/HMPID/AliHMPIDParam.h
+++ b/HMPID/AliHMPIDParam.h
@@ -31,7 +31,7 @@ public:
   void     Lors2Mars   (Int_t c,Float_t x,Float_t y,Double_t *m,Int_t pl=kPc)const{Double_t z=0; switch(pl){case kPc:z=8.0;break; case kAnod:z=7.806;break; case kRad:z=-1.25; break;}   Double_t l[3]={x-fX,y-fY,z};  fM[c]->LocalToMaster(l,m); }    
   TVector3 Lors2Mars   (Int_t c,Float_t x,Float_t y,            Int_t pl=kPc)const{Double_t m[3];Lors2Mars(c,x,y,m,pl); return TVector3(m);    }//MRS->LRS  
   void     Mars2Lors   (Int_t c,Double_t *m,Float_t &x,Float_t &y           )const{Double_t l[3];fM[c]->MasterToLocal(m,l);x=l[0]+fX;y=l[1]+fY;}//MRS->LRS
-  void     Mars2LorsVec(Int_t c,Double_t *m,Float_t &th,Float_t &ph         )const{Double_t l[3]; fM[c]->MasterToLocalVect(m,l); Float_t pt=TMath::Sqrt(l[0]*l[0]+l[1]*l[1]); th=TMath::ATan(l[2]/pt); ph=TMath::ATan(l[0]/pt);}    
+  void     Mars2LorsVec(Int_t c,Double_t *m,Float_t &th,Float_t &ph         )const{Double_t l[3]; fM[c]->MasterToLocalVect(m,l); Float_t pt=TMath::Sqrt(l[0]*l[0]+l[1]*l[1]); th=TMath::ATan(pt/l[2]); ph=TMath::ATan(l[1]/l[0]);}    
   TVector3 Norm        (Int_t c                                             )const{Double_t n[3]; Norm(c,n); return TVector3(n);               }//norm 
   void     Norm        (Int_t c,Double_t *n                                 )const{Double_t l[3]={0,0,1};fM[c]->LocalToMasterVect(l,n);        }//norm
   

diff --git a/HMPID/AliHMPIDRecon.h b/HMPID/AliHMPIDRecon.h
index 5c0cb76c0f97506d74b732e4b06fc29d5b51d2a0..b9ec332db7e00d8f7bbb63ed032a4e47e1e9df48 100644 (file)
--- a/HMPID/AliHMPIDRecon.h
+++ b/HMPID/AliHMPIDRecon.h
@@ -36,7 +36,7 @@ public :
   void     Propagate    (const TVector3 &dir,      TVector3 &pos,Double_t z                 )const;//propagate photon alogn the line  
   void     Refract      (      TVector3 &dir,                    Double_t n1,    Double_t n2)const;//refract photon on the boundary
   Double_t TracePhot    (Double_t ckovTh,Double_t ckovPh,TVector2 &pos                      )const;//trace photon created by track to PC 
-  void     SetTrack     (Double_t x,Double_t y,Double_t th,Double_t ph                      ){fTrkDir.SetMagThetaPhi(1,th,ph);  fTrkPos.Set(x,y);}//set track
+  void     SetTrack     (Double_t x,Double_t y,Double_t theta,Double_t phi                  ){fTrkDir.SetMagThetaPhi(1,theta,phi);  fTrkPos.Set(x,y);}//set track
   Double_t SigLoc       (Double_t ckovTh,Double_t ckovPh,Double_t beta                      )const;//error due to cathode segmetation
   Double_t SigGeom      (Double_t ckovTh,Double_t ckovPh,Double_t beta                      )const;//error due to unknown photon origin
   Double_t SigCrom      (Double_t ckovTh,Double_t ckovPh,Double_t beta                      )const;//error due to unknonw photon energy

diff --git a/HMPID/AliHMPIDTracker.cxx b/HMPID/AliHMPIDTracker.cxx
index f1d97dc622e4fbd835eba7946b76c8920ce1beab..2ed03e7a23bd85cf819cacda9ff9ea4486e3894f 100644 (file)
--- a/HMPID/AliHMPIDTracker.cxx
+++ b/HMPID/AliHMPIDTracker.cxx
@@ -82,20 +82,20 @@ Int_t AliHMPIDTracker::IntTrkCha(AliESDtrack *pTrk,Float_t &xRa,Float_t &yRa)
 //   Returns: intersected chamber ID or -1
   
   AliHMPIDParam *pParam=AliHMPIDParam::Instance();
-  Float_t xPc=0,yPc=0,th=0,ph=0;                                                                //track intersection point and angles, LORS  
-  for(Int_t i=0;i<7;i++){                                                                       //chambers loop
+  Float_t xPc=0,yPc=0,theta=0,phi=0;                                                            //track intersection point and angles, LORS  
+  for(Int_t i=AliHMPIDDigit::kMinCh;i<=AliHMPIDDigit::kMaxCh;i++){                                                                       //chambers loop
     Double_t p1[3],n1[3]; pParam->Norm(i,n1); pParam->Lors2Mars(i,0,0,p1,AliHMPIDParam::kRad);  //point & norm  for RAD
     Double_t p2[3],n2[3]; pParam->Norm(i,n2); pParam->Lors2Mars(i,0,0,p2,AliHMPIDParam::kPc);   //point & norm  for PC
       
     if(pTrk->Intersect(p1,n1,-GetBz())==kFALSE) continue;                                       //try to intersect track with the middle of radiator
     if(pTrk->Intersect(p2,n2,-GetBz())==kFALSE) continue;                                       //try to intersect track with PC
       
-    pParam->Mars2LorsVec(i,n1,th,ph);                                                           //track angles
+    pParam->Mars2LorsVec(i,n1,theta,phi);                                                       //track angles
     pParam->Mars2Lors   (i,p1,xRa,yRa);                                                         //TRKxRAD position
     pParam->Mars2Lors   (i,p2,xPc,yPc);                                                         //TRKxPC position
       
     if(AliHMPIDDigit::IsInside(xPc,yPc)==kFALSE) continue;                                      //not in active area  
-    pTrk->SetHMPIDtrk      (xPc,yPc,th,ph);                                                     //store track intersection info
+    pTrk->SetHMPIDtrk      (xPc,yPc,theta,phi);                                                 //store track intersection info
     pTrk->SetHMPIDcluIdx   (i,0);
     return i;
   }                                                                                             //chambers loop

diff --git a/HMPID/Hdisp.C b/HMPID/Hdisp.C
index 0a41475d082d18b88c099894394d36abb5710c8d..8f8619774f126c1c6d27b028afcd3a8021adbcad 100644 (file)
--- a/HMPID/Hdisp.C
+++ b/HMPID/Hdisp.C
@@ -4,7 +4,8 @@ Int_t gEvt=0; Int_t gMaxEvt=0;
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 void Hdisp()
 {//display events from files if any in current directory or simulated events
-  pAll=new TCanvas("all","",1000,900); pAll->Divide(3,3,0,0);pAll->ToggleEditor();
+  pAll=new TCanvas("all","",1300,900); pAll->Divide(3,3,0,0);
+//  pAll->ToggleEditor();
   pAll->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)",0,"","DoZoom(Int_t,Int_t,Int_t,TObject*)");
   
   if(gSystem->IsFileInIncludePath("galice.root")){// tries to open session
@@ -48,13 +49,14 @@ void SimEvt()
   TObjArray    digs(7); for(Int_t i=0;i<7;i++) digs.AddAt(new TClonesArray("AliHMPIDDigit"),i);
   TObjArray    clus(7); for(Int_t i=0;i<7;i++) clus.AddAt(new TClonesArray("AliHMPIDCluster"),i);
   AliESD esd;
-    
+  AliHMPIDDigit::fSigmas=3;
+  AliHMPIDDigitizer::DoNoise(kTRUE);
   SimEsd(&esd);
   SimHits(&esd,&hits);
              AliHMPIDv1::Hit2Sdi(&hits,&sdig);                               
       AliHMPIDDigitizer::Sdi2Dig(&sdig,&digs);     
-  AliHMPIDReconstructor::Dig2Clu(&digs,&clus);
-        AliHMPIDTracker::Recon(&esd,&clus);
+      AliHMPIDReconstructor::Dig2Clu(&digs,&clus);
+      AliHMPIDTracker::Recon(&esd,&clus,OpenCalib());
   
   pAll->cd(3);  gPad->Clear(); TLatex txt;txt.DrawLatex(0.2,0.2,Form("Simulated event %i",gEvt));
   DrawEvt(&hits,&digs,&clus,&esd);  
@@ -64,12 +66,12 @@ void SimEvt()
 void SimEsd(AliESD *pEsd)
 {
   TParticle part; TLorentzVector mom;
-  for(Int_t iTrk=0;iTrk<25;iTrk++){//stack loop
+  for(Int_t iTrk=0;iTrk<100;iTrk++){//stack loop
     part.SetPdgCode(kProton);
-    part.SetProductionVertex(0,0,0,0);  
-    Double_t eta= -0.4+gRandom->Rndm()*0.8; //rapidity is random [-0.4,+0.4]
-    Double_t phi= gRandom->Rndm()*1.4;      //phi is random      [ 0  , 80 ] degrees    
-    mom.SetPtEtaPhiM(2,eta,phi,part.GetMass());
+    part.SetProductionVertex(0,0,0,0);
+    Double_t eta= -0.2+gRandom->Rndm()*0.4; //rapidity is random [-0.2,+0.2]
+    Double_t phi= gRandom->Rndm()*60.*TMath::DegToRad();   //phi is random      [ 0  , 60 ] degrees    
+    mom.SetPtEtaPhiM(5,eta,phi,part.GetMass());
     part.SetMomentum(mom);
     AliESDtrack trk(&part);
     pEsd->AddTrack(&trk);
@@ -78,8 +80,8 @@ void SimEsd(AliESD *pEsd)
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 void SimHits(AliESD *pEsd, TClonesArray *pHits)
 {
-  const Int_t kCerenkov=AliHMPIDParam::CkovCode();
-  const Int_t kFeedback=AliHMPIDParam::FeedCode();
+  const Int_t kCerenkov=50000050;
+  const Int_t kFeedback=50000051;
   
   AliHMPIDRecon rec;
   
@@ -89,15 +91,24 @@ void SimHits(AliESD *pEsd, TClonesArray *pHits)
     Float_t xRa,yRa;
     Int_t ch=AliHMPIDTracker::IntTrkCha(pTrk,xRa,yRa);
     if(ch<0) continue; //this track does not hit HMPID
-    Float_t ckov=0.63;
+    Float_t beta = pTrk->GetP()/(TMath::Sqrt(pTrk->GetP()*pTrk->GetP()+0.938*0.938));
+    Float_t ckov=TMath::ACos(1./(beta*1.292));
 
-    Float_t th,ph,xPc,yPc,; pTrk->GetHMPIDtrk(xPc,yPc,th,ph); rec.SetTrack(xRa,yRa,th,ph); 
+    Float_t theta,phi,xPc,yPc,; 
+    pTrk->GetHMPIDtrk(xPc,yPc,theta,phi); 
+    rec.SetTrack(xRa,yRa,theta,phi); 
     
     if(!AliHMPIDDigit::IsInDead(xPc,yPc)) new((*pHits)[hc++]) AliHMPIDHit(ch,200e-9,kProton  ,iTrk,xPc,yPc);                 //mip hit
-    for(int i=0;i<4;i++)  new((*pHits)[hc++]) AliHMPIDHit(ch,7.5e-9,kFeedback,iTrk,gRandom->Rndm()*130,gRandom->Rndm()*126); //bkg hits 4 per track
-    for(int i=0;i<16;i++){
+    /*
+    for(int i=0;i<4;i++) { 
+      Float_t x=gRandom->Rndm()*130;Float_t y=gRandom->Rndm()*126;
+      if(!AliHMPIDDigit::IsInDead(x,y)) new((*pHits)[hc++]) AliHMPIDHit(ch,7.5e-9,kFeedback,iTrk,x,y); //bkg hits 4 per track
+    }
+    */
+    Int_t nPhots = (Int_t)(20.*TMath::Power(TMath::Sin(ckov),2)/TMath::Power(TMath::Sin(TMath::ACos(1./1.292)),2));
+    for(int i=0;i<nPhots;i++){
       rec.TracePhot(ckov,gRandom->Rndm()*TMath::TwoPi(),pos);
-      new((*pHits)[hc++]) AliHMPIDHit(ch,7.5e-9,kCerenkov,iTrk,pos.X(),pos.Y());
+      if(!AliHMPIDDigit::IsInDead(pos.X(),pos.Y())) new((*pHits)[hc++]) AliHMPIDHit(ch,7.5e-9,kCerenkov,iTrk,pos.X(),pos.Y());
     }                      //photon hits  
   }//tracks loop    
 }//SimHits()
@@ -156,7 +167,11 @@ void DrawEvt(TClonesArray *pHitLst,TObjArray *pDigLst,TObjArray *pCluLst,AliESD
     Float_t ckov=pTrk->GetHMPIDsignal();  Float_t err=TMath::Sqrt(pTrk->GetHMPIDchi2());
     if(ckov>0){
       rec.SetTrack(xPc,yPc,th,ph);
-     TVector2 pos;  for(int j=0;j<100;j++){rec.TracePhot(ckov,j*0.0628,pos); pRin[ch]->SetNextPoint(pos.X(),pos.Y());}      
+      TVector2 pos;  
+      for(int j=0;j<100;j++){
+       rec.TracePhot(ckov,j*0.0628,pos); 
+       if(!AliHMPIDDigit::IsInDead(pos.X(),pos.Y())) pRin[ch]->SetNextPoint(pos.X(),pos.Y());
+      }      
     }
   }//tracks loop
       
@@ -182,7 +197,7 @@ void DrawEvt(TClonesArray *pHitLst,TObjArray *pDigLst,TObjArray *pCluLst,AliESD
     ((TClonesArray*)pDigLst->At(iCh))->Draw();  //draw digits
     ((TClonesArray*)pCluLst->At(iCh))->Draw();  //draw clusters
                             pTxC[iCh]->Draw();  //draw intersections
-                            pRin[iCh]->Draw();  //draw rings
+                            pRin[iCh]->Draw("CLP");  //draw rings
     gPad->SetEditable(kFALSE);
   }//chambers loop
   
@@ -207,7 +222,7 @@ void DoZoom(Int_t evt, Int_t px, Int_t py, TObject *obj)
   if(evt==6&&zoom==minZoom) return; 
   
  // if(!obj->IsA()->InheritsFrom("TPad")) return;  //current object is not pad
-  TVirtualPad *pPad=gROOT->GetSelectedPad();
+  TVirtualPad *pPad=gPad->GetSelectedPad();
   if(pPad->GetNumber()==3 || pPad->GetNumber()==7) return; //current pad is wrong
 
  // Printf("evt=%i (%i,%i) %s",evt,px,py,obj->GetName());
@@ -262,3 +277,17 @@ void PrintClus()
   
   Printf("totally %i clusters for event %i",iCluCnt,gEvt);
 }
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+TObjArray* OpenCalib()
+{
+  AliCDBManager* pCDB = AliCDBManager::Instance();
+  pCDB->SetDefaultStorage("local://$HOME");
+  AliCDBEntry *pQthreEnt=pCDB->Get("HMPID/Calib/Qthre",0);
+  AliCDBEntry *pNmeanEnt=pCDB->Get("HMPID/Calib/Nmean",0);
+  
+  if(!pQthreEnt || ! pNmeanEnt) return;
+  
+  TObjArray *pNmean=(TObjArray*)pNmeanEnt->GetObject(); 
+  TObjArray *pQthre=(TObjArray*)pQthreEnt->GetObject(); 
+  return pNmean;
+}