new 2d event display + optical properties restored in v1
authorkir <kir@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 23 Feb 2007 11:32:03 +0000 (11:32 +0000)
committerkir <kir@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 23 Feb 2007 11:32:03 +0000 (11:32 +0000)
HMPID/AliHMPIDRecon.cxx
HMPID/AliHMPIDRecon.h
HMPID/AliHMPIDTracker.cxx
HMPID/AliHMPIDTracker.h
HMPID/AliHMPIDv1.cxx
HMPID/Hmenu.C

index 99905b5..865dd8c 100644 (file)
 #include <TRotation.h>     //TracePhoton()
 #include <TH1D.h>          //HoughResponse()
 #include <TClonesArray.h>  //CkovAngle()
-
-#include <TTree.h>         //Display()
-#include <TFile.h>         //Display()
-#include <AliESD.h>        //Display()
-#include <TPolyMarker.h>   //Display()
-#include <TLatex.h>        //Display()
-#include <TCanvas.h>       //Display()
-
+#include <AliESDtrack.h>   //CkovAngle()
 
 const Double_t AliHMPIDRecon::fgkRadThick=1.5;
 const Double_t AliHMPIDRecon::fgkWinThick=0.5;
@@ -61,29 +54,46 @@ AliHMPIDRecon::AliHMPIDRecon():TTask("RichRec","RichPat"),
   }
 }
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Double_t AliHMPIDRecon::CkovAngle(TClonesArray *pCluLst,Int_t &iNaccepted)
+void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst)
 {
 // Pattern recognition method based on Hough transform
-// Arguments: pCluLst  - list of clusters for this chamber   
-//   Returns:          - track ckov angle, [rad], 
+// Arguments:   pTrk     - track for which Ckov angle is to be found
+//              pCluLst  - list of clusters for this chamber   
+//   Returns:            - track ckov angle, [rad], 
   
   if(pCluLst->GetEntries()>200) fIsWEIGHT = kTRUE; // offset to take into account bkg in reconstruction
   else                          fIsWEIGHT = kFALSE;
 
   // Photon Flag:  Flag = 0 initial set; Flag = 1 good candidate (charge compatible with photon); Flag = 2 photon used for the ring;
+  Float_t xPc,yPc,th,ph;      pTrk->GetHMPIDtrk(xPc,yPc,th,ph);  SetTrack(xPc,yPc,th,ph); //initialize this track            
 
+  
+  
+  Float_t dMin=999,mipX=-1,mipY=-1;Int_t chId=-1,mipId=-1,mipQ=-1;                                                                           
   fPhotCnt=0;                                                      
   for (Int_t iClu=0; iClu<pCluLst->GetEntriesFast();iClu++){//clusters loop
     AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu);                       //get pointer to current cluster    
-    if(pClu->Q()>100) continue;                                                             //avoid MIP clusters from bkg
-    
-    fPhotCkov[fPhotCnt]=FindPhotCkov(pClu->X(),pClu->Y());                                  //find ckov angle for this  photon candidate
-    fPhotCnt++;         //increment counter of photon candidates
+    chId=pClu->Ch();
+    if(pClu->Q()>100){                                                                        //charge compartible with MIP clusters      
+      Float_t dX=xPc-pClu->X(),dY=yPc-pClu->Y(),d =TMath::Sqrt(dX*dX+dY*dY);                  //distance between current cluster and intersection point
+      if( d < dMin) {mipId=iClu; dMin=d;mipX=pClu->X();mipY=pClu->Y();mipQ=(Int_t)pClu->Q();}//current cluster is closer, overwrite data for min cluster
+    }else{                                                                                    //charge compartible with photon cluster
+      fPhotCkov[fPhotCnt]=FindPhotCkov(pClu->X(),pClu->Y());                                  //find ckov angle for this  photon candidate
+      fPhotCnt++;         //increment counter of photon candidates
+    }
   }//clusters loop
+  Int_t iNacc=FlagPhot(HoughResponse());                                   //flag photons according to individual theta ckov with respect to most probable
+
+                 pTrk->SetHMPIDmip      (mipX,mipY,mipQ,iNacc);                 //store mip info 
+
+  if(mipId==-1) {pTrk->SetHMPIDsignal   (kMipQdcCut);  return;}                 //no clusters with QDC more the threshold at all
+  if(dMin>1)    {pTrk->SetHMPIDsignal   (kMipDistCut); return;}                 //closest cluster with enough charge is still too far from intersection
+                 pTrk->SetHMPIDcluIdx(chId,mipId); 
+  if(iNacc<1)    pTrk->SetHMPIDsignal(kNoPhotAccept);                         //no photon candidates is accepted  
+  else           pTrk->SetHMPIDsignal(FindRingCkov(pCluLst->GetEntries()));   //find best Theta ckov for ring i.e. track
+  
+                 pTrk->SetHMPIDchi2(fCkovSigma2);                              //error squared 
 
-  iNaccepted=FlagPhot(HoughResponse()); //flag photons according to individual theta ckov with respect to most probable track theta ckov
-  if(iNaccepted<1) return -11; 
-  else             return FindRingCkov(pCluLst->GetEntries());  //find best Theta ckov for ring i.e. track
 }//ThetaCerenkov()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 Double_t AliHMPIDRecon::FindPhotCkov(Double_t cluX,Double_t cluY)
index f8bddad..67eef64 100644 (file)
@@ -16,7 +16,9 @@
 #include <TTask.h>        //base class
 #include <TVector3.h>     //fields 
 
-class TClonesArray;
+class TClonesArray; //CkovAngle()
+class AliESDtrack;  //CkovAngle()
+
 class AliHMPIDRecon : public TTask 
 {
 public : 
@@ -24,8 +26,7 @@ public :
     virtual ~AliHMPIDRecon()                                                          {}
 
   
-  Double_t CkovAngle    (TClonesArray *pCluLst,Int_t &iNaccepted);                                                         //reconstructed Theta Cerenkov
-  Double_t CkovSigma2   (                                                                   )const{ return fCkovSigma2;} //track ckov angle error squared
+  void     CkovAngle    (AliESDtrack *pTrk,TClonesArray *pCluLst                            );                           //reconstructed Theta Cerenkov
   Double_t FindPhotCkov (Double_t cluX,Double_t cluY                                        );     //find ckov angle for single photon candidate
   Double_t FindPhotPhi  (Double_t cluX,Double_t cluY                                        );     //find phi angle for single photon candidate
   Double_t FindRingCkov (Int_t iNclus                                                       );     //best ckov for ring formed by found photon candidates
@@ -35,11 +36,12 @@ 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 th,Double_t ph,Double_t x,Double_t y                      ){fTrkDir.SetMagThetaPhi(1,th,ph);  fTrkPos.Set(x,y);}//set track
+  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
   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
   Double_t Sigma2       (Double_t ckovTh,Double_t ckovPh                                    )const;//photon candidate sigma
+  enum ETrackingFlags {kMipDistCut=-9,kMipQdcCut=-5,kNoPhotAccept=-11};
 protected:
   static const Double_t fgkRadThick;                      //radiator thickness
   static const Double_t fgkWinThick;                      //window thickness
index 064eebc..1b3e4f5 100644 (file)
@@ -41,62 +41,43 @@ Int_t AliHMPIDTracker::Recon(AliESD *pESD,TObjArray *pCluAll)
 // Arguments: pESD - pointer to Event Summary Data class instance which contains a list of tracks
 //   Returns: error code, 0 if no errors   
   Int_t iNtracks=pESD->GetNumberOfTracks();  AliDebugClass(1,Form("Start with %i tracks",iNtracks));
-  AliHMPIDRecon recon;                                                        //instance of reconstruction class, nothing important in ctor
-   
-  AliHMPIDParam *pParam=AliHMPIDParam::Instance();
-
+  
+  AliHMPIDRecon recon;                                                                           //instance of reconstruction class, nothing important in ctor
+  Float_t xRa,yRa;
   for(Int_t iTrk=0;iTrk<iNtracks;iTrk++){                                                        //ESD tracks loop
     AliESDtrack *pTrk = pESD->GetTrack(iTrk);                                                    //get next reconstructed track    
-        
-    Float_t xRa=0,yRa=0,xPc=0,yPc=0,th=0,ph=0;                                                   //track intersection point and angles, LORS  
-    Int_t iCh=-1;                                                                                //intersected chamber 
-    for(Int_t i=0;i<7;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->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  
-      iCh=i;
-      break;
-    }//chambers loop      
-    
-    if(iCh==-1) continue;                                                                  //no intersection at all, go after next track
-    
-    TClonesArray *pCluLst=(TClonesArray *)pCluAll->At(iCh);                                //get clusters list for intersected chamber
-    
-    Double_t    dMin=999;                                                                  //distance between track-PC intersection point and current cluster
-    Int_t   iMip=-1;                                                                       //index of cluster nearest to intersection point
-    for(Int_t iClu=0;iClu<pCluLst->GetEntries();iClu++){                                   //clusters loop for intersected chamber
-      AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->At(iClu);                             //get pointer to current cluster
-      if(pClu->Q()<100) continue;                                                          //QDC is incompartible with mip, go after another one
-      
-      Float_t dX=xPc-pClu->X();                                                            //distance between current cluster and intersection point
-      Float_t dY=yPc-pClu->Y();
-      Float_t d =TMath::Sqrt(dX*dX+dY*dY);
-      
-      if( d < dMin) {iMip=iClu; dMin=d;}                                                   //current cluster is closer, overwrite data for min cluster
-    }//clusters loop for intersected chamber    
-    
-                   pTrk->SetHMPIDtrk      (xPc,yPc,th,ph);                                  //store track info
-    if(iMip==-1)  {pTrk->SetHMPIDsignal   (kMipQdcCut);  continue;}                         //no clusters with QDC more the threshold at all
-    
-    AliHMPIDCluster *pMipClu=(AliHMPIDCluster*)pCluLst->At(iMip);                            //take mip cluster 
-    
-                   pTrk->SetHMPIDmip      ((Float_t)pMipClu->X(),(Float_t)pMipClu->Y(),(Int_t)pMipClu->Q());          //store mip info 
-    if(dMin>1)    {pTrk->SetHMPIDsignal   (kMipDistCut); continue;}                          //closest cluster with enough charge is still too far 
-                   pTrk->SetHMPIDcluIdx   (iCh,iMip);                                        //set mip cluster index
-  recon.SetTrack(th,ph,xRa,yRa); Int_t iNphot=0;                                            //initialize track parameters  
-                   pTrk->SetHMPIDsignal   (recon.CkovAngle(pCluLst,iNphot));                 //search for Cerenkov angle for this track
-                   pTrk->SetHMPIDchi2     (recon.CkovSigma2());                              //error squared 
-                   pTrk->SetHMPIDmip      ((Float_t)pMipClu->X(),(Float_t)pMipClu->Y(),(Int_t)pMipClu->Q(),iMip);     //info on mip cluster + n. phot.
- }//ESD tracks loop
+    Int_t cham=IntTrkCha(pTrk,xRa,yRa);                                                          //get chamber intersected by thie track 
+    if(cham<0) continue;                                                                         //no intersection at all, go after next track
+    recon.CkovAngle(pTrk,(TClonesArray *)pCluAll->At(cham));                                     //search for Cerenkov angle for this track
+  }                                                                                              //ESD tracks loop
   AliDebugClass(1,"Stop pattern recognition");
   return 0; // error code: 0=no error;
 }//PropagateBack()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Int_t AliHMPIDTracker::IntTrkCha(AliESDtrack *pTrk,Float_t &xRa,Float_t &yRa)
+{
+// Static method to find intersection in between given track and HMPID chambers
+// Arguments: pTrk    - ESD track 
+//            xRa,yRa - track intersection with the middle of radiator, LORS
+//   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
+    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->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->SetHMPIDcluIdx   (i,0);
+    return i;
+  }                                                                                             //chambers loop
+  return -1; //no intersection with HMPID chambers
+}
index e1f836f..d1bf4f3 100644 (file)
@@ -5,8 +5,8 @@
 #include "AliHMPID.h"   //Recon()
 #include <AliRun.h>     //Recon()
 
-class AliESD;           
-
+class AliESD;      //Recon()     
+class AliESDtrack; //IntTrkCha()
 class AliHMPIDTracker : public AliTracker
 {
 public:
@@ -21,8 +21,8 @@ public:
          Int_t       RefitInward    (AliESD *                   )       {return 0;} //pure virtual from AliTracker 
          void        UnloadClusters (                           )       {         } //pure virtual from AliTracker 
 //private part  
-  enum ETrackingFlags {kMipDistCut=-9,kMipQdcCut=-5};
-  static Int_t       Recon(AliESD *pEsd,TObjArray *pCluAll);                        //do actual job 
+  static Int_t       IntTrkCha(AliESDtrack *pTrk,Float_t &x,Float_t &y);                    //find track-chamber intersection, retuns chamber ID
+  static Int_t       Recon    (AliESD *pEsd,TObjArray *pCluAll        );                    //do actual job, returns status code  
 protected:
   ClassDef(AliHMPIDTracker,0)
 };//class AliHMPIDTracker
index 2e679f2..1df7678 100644 (file)
@@ -91,6 +91,9 @@ void AliHMPIDv1::CreateMaterials()
     AliMaterial(++matId,"W"   ,aW   ,zW   ,dW   ,radW   ,absW   );  AliMedium(kW   ,"W"   , matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
     AliMaterial(++matId,"Al"  ,aAl  ,zAl  ,dAl  ,radAl  ,absAl  );  AliMedium(kAl  ,"Al"  , matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
     
+    
+    DefineOpticalProperties();
+    
     AliDebug(1,"Stop v1 HMPID.");
 
   TString ttl=GetTitle(); if(!ttl.Contains("ShowOptics")) return;              //user didn't aks to plot optical curves
@@ -277,16 +280,16 @@ void AliHMPIDv1::DefineOpticalProperties()
   Float_t aIdxRad[kNbins], aIdxWin[kNbins], aIdxGap[kNbins], aIdxMet[kNbins], aIdxPc[kNbins]; 
   Float_t                                                    aQeAll [kNbins], aQePc [kNbins];
 
-  TF2 *pRaIF=new TF2("RidxRad","sqrt(1+0.554*(1239.84/x)^2/((1239.84/x)^2-5796)-0.0005*(y-20))"                                       ,emin,emax,0,50); //DiMauro mail temp 0-50 degrees C
-  TF1 *pWiIF=new TF1("RidxWin","sqrt(1+46.411/(10.666*10.666-x*x)+228.71/(18.125*18.125-x*x))"                                        ,emin,emax);      //SiO2 idx TDR p.35
-  TF1 *pGaIF=new TF1("RidxGap","1+0.12489e-6/(2.62e-4 - x*x/1239.84/1239.84)"                                                         ,emin,emax);      //?????? from where  
+  TF2 *pRaIF=new TF2("HidxRad","sqrt(1+0.554*(1239.84/x)^2/((1239.84/x)^2-5796)-0.0005*(y-20))"                                       ,emin,emax,0,50); //DiMauro mail temp 0-50 degrees C
+  TF1 *pWiIF=new TF1("HidxWin","sqrt(1+46.411/(10.666*10.666-x*x)+228.71/(18.125*18.125-x*x))"                                        ,emin,emax);      //SiO2 idx TDR p.35
+  TF1 *pGaIF=new TF1("HidxGap","1+0.12489e-6/(2.62e-4 - x*x/1239.84/1239.84)"                                                         ,emin,emax);      //?????? from where  
 
-  TF1 *pRaAF=new TF1("RabsRad","(x<7.8)*(gaus+gaus(3))+(x>=7.8)*0.0001"                                                               ,emin,emax);  //fit from DiMauro data 28.10.03 
+  TF1 *pRaAF=new TF1("HabsRad","(x<7.8)*(gaus+gaus(3))+(x>=7.8)*0.0001"                                                               ,emin,emax);  //fit from DiMauro data 28.10.03 
   pRaAF->SetParameters(3.20491e16,-0.00917890,0.742402,3035.37,4.81171,0.626309);
-  TF1 *pWiAF=new TF1("RabsWin","(x<8.2)*(818.8638-301.0436*x+36.89642*x*x-1.507555*x*x*x)+(x>=8.2)*0.0001"                            ,emin,emax);  //fit from DiMauro data 28.10.03 
-  TF1 *pGaAF=new TF1("RabsGap","(x<7.75)*6512.399+(x>=7.75)*3.90743e-2/(-1.655279e-1+6.307392e-2*x-8.011441e-3*x*x+3.392126e-4*x*x*x)",emin,emax);  //????? from where  
+  TF1 *pWiAF=new TF1("HabsWin","(x<8.2)*(818.8638-301.0436*x+36.89642*x*x-1.507555*x*x*x)+(x>=8.2)*0.0001"                            ,emin,emax);  //fit from DiMauro data 28.10.03 
+  TF1 *pGaAF=new TF1("HabsGap","(x<7.75)*6512.399+(x>=7.75)*3.90743e-2/(-1.655279e-1+6.307392e-2*x-8.011441e-3*x*x+3.392126e-4*x*x*x)",emin,emax);  //????? from where  
   
-  TF1 *pQeF =new TF1("Qe"    ,"0+(x>6.07267)*0.344811*(1-exp(-1.29730*(x-6.07267)))"                                                  ,emin,emax);  //fit from DiMauro data 28.10.03  
+  TF1 *pQeF =new TF1("Hqe"    ,"0+(x>6.07267)*0.344811*(1-exp(-1.29730*(x-6.07267)))"                                                 ,emin,emax);  //fit from DiMauro data 28.10.03  
                             
   for(Int_t i=0;i<kNbins;i++){
     Float_t eV=emin+0.1*i;  //Ckov energy in eV
index 9e6c56c..f695400 100644 (file)
@@ -1,15 +1,15 @@
-AliRun *a;     AliRunLoader *al;   TGeoManager *g; //globals for easy manual manipulations
-AliHMPID   *r; AliLoader    *rl; AliHMPIDParam *rp;
+AliRun     *a; AliRunLoader *al;   TGeoManager *g; //globals for easy manual manipulations
+AliHMPID   *h; AliLoader    *hl; AliHMPIDParam *hp;
 Bool_t isGeomType=kFALSE;
 
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 void GetParam()
 {
   isGeomType=!isGeomType;
-  if(g) delete g;  if(rp) delete rp; //delete current TGeoManager and AliHMPIDParam
+  if(g) delete g;  if(hp) delete hp; //delete current TGeoManager and AliHMPIDParam
   if(isGeomType) g=TGeoManager::Import("geometry.root");
   else           g=TGeoManager::Import("misaligned_geometry.root");
-  rp=AliHMPIDParam::Instance();
+  hp=AliHMPIDParam::Instance();
 }//GetParam()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 void Hmenu()
@@ -20,9 +20,9 @@ void Hmenu()
     al=AliRunLoader::Open();                                                //try to open galice.root from current dir 
     if(gAlice) delete gAlice;                                               //in case we execute this in aliroot delete default AliRun object 
     al->LoadgAlice(); a=al->GetAliRun();                                    //take new AliRun object from galice.root   
-    rl=al->GetDetectorLoader("HMPID");  r=(AliHMPID*)a->GetDetector("HMPID");  //get HMPID object from galice.root
+    hl=al->GetDetectorLoader("HMPID");  h=(AliHMPID*)a->GetDetector("HMPID");  //get HMPID object from galice.root
     
-    status+=Form(" with %i event(s)",al->GetNumberOfEvents()); status+=(r)? " with HMPID": " without HMPID";
+    status+=Form(" with %i event(s)",al->GetNumberOfEvents()); status+=(h)? " with HMPID": " without HMPID";
   }else  
     status+="No galice.root";
   
@@ -45,23 +45,24 @@ void Hmenu()
 void General()
 {         
   TControlBar *pMenu = new TControlBar("vertical","General purpose",100,50);  
-    pMenu->AddButton("       Debug ON      ","don();"                    ,"Switch debug on-off"                        );   
-    pMenu->AddButton("       Debug OFF     ","doff();"                   ,"Switch debug on-off"                        );   
-    pMenu->AddButton("        Geo GUI      ","geo();"                    ,"Shows geometry"                             ); 
-    pMenu->AddButton("        Browser      ","new TBrowser;"             ,"Start ROOT TBrowser"                        );
+    pMenu->AddButton("Debug ON","don();"                    ,"Switch debug on-off"                        );   
+    pMenu->AddButton("Debug OFF","doff();"                   ,"Switch debug on-off"                        );   
+    pMenu->AddButton("Geo GUI","geo();"                    ,"Shows geometry"                             ); 
+    pMenu->AddButton("Browser","new TBrowser;"             ,"Start ROOT TBrowser"                        );
   pMenu->Show();  
 }//General()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 void SimData()
 {
-  TControlBar *pMenu = new TControlBar("vertical","Sim data",340,50);  
-    pMenu->AddButton("Display ALL chambers ","ed();"      ,"Display Fast");
-    pMenu->AddButton("      HITS QA        ","hqa()"      ,"QA plots for hits: hqa()");
-    pMenu->AddButton("     Print Stack     ","stack();"   ,"To print hits:     hp(evt)");
-    pMenu->AddButton("     Print hits      ","hp();"      ,"To print hits:     hp(evt)");
-    pMenu->AddButton("    Print sdigits    ","sp();"      ,"To print sdigits:  sp(evt)");
-    pMenu->AddButton("    Print digits     ","dp();"      ,"To print digits:   dp(evt)");
-    pMenu->AddButton("   Print clusters    ","cp();"      ,"To print clusters: cp(evt)");
+  TControlBar *pMenu = new TControlBar("vertical","Sim data",310,50);  
+    pMenu->AddButton("Display from files","hed();"    ,"Display Fast");
+    pMenu->AddButton("Display simulated" ,"sed();"    ,"Display Fast");
+    pMenu->AddButton("HITS QA"           ,"hqa()"     ,"QA plots for hits: hqa()");
+    pMenu->AddButton("Print stack"       ,"stack();"  ,"To print hits:     hp(evt)");
+    pMenu->AddButton("Print hits"        ,"hp();"     ,"To print hits:     hp(evt)");
+    pMenu->AddButton("Print sdigits"     ,"sp();"     ,"To print sdigits:  sp(evt)");
+    pMenu->AddButton("Print digits"      ,"dp();"     ,"To print digits:   dp(evt)");
+    pMenu->AddButton("Print clusters"    ,"cp();"     ,"To print clusters: cp(evt)");
   pMenu->Show();         
 }//SimData()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -96,15 +97,15 @@ void don() {  Printf("DebugON");   AliLog::SetGlobalDebugLevel(AliLog::kDebug);}
 
 void geo (                       ) {gGeoManager->GetTopVolume()->Draw("ogl");}
   
-void du  (                       ) {r->Dump         (   );}                //utility display 
+void du  (                       ) {h->Dump         (   );}                //utility display 
 
-void hp  (Int_t evt=0            ) {r->HitPrint  (evt);}   //print hits for requested event
-void hq  (                       ) {r->HitQA     (   );}   //hits QA plots for all events 
-void sp  (Int_t evt=0            ) {r->SdiPrint  (evt);}   //print sdigits for requested event
-void sq  (Int_t evt=0            ) {r->SdiPrint  (evt);}   //print sdigits for requested event
-void dp  (Int_t evt=0            ) {r->DigPrint  (evt);}   //print digits for requested event
+void hp  (Int_t evt=0            ) {h->HitPrint  (evt);}   //print hits for requested event
+void hq  (                       ) {h->HitQA     (   );}   //hits QA plots for all events 
+void sp  (Int_t evt=0            ) {h->SdiPrint  (evt);}   //print sdigits for requested event
+void sq  (Int_t evt=0            ) {h->SdiPrint  (evt);}   //print sdigits for requested event
+void dp  (Int_t evt=0            ) {h->DigPrint  (evt);}   //print digits for requested event
 void dq  (                       ) {AliHMPIDReconstructor::DigQA     (al );}   //digits QA plots for all events
-void cp  (Int_t evt=0            ) {r->CluPrint  (evt);                   }   //print clusters for requested event
+void cp  (Int_t evt=0            ) {h->CluPrint  (evt);                   }   //print clusters for requested event
 void cq  (                       ) {AliHMPIDReconstructor::CluQA     (al );}   //clusters QA plots for all events
 
 void ep  (                       ) {AliHMPIDTracker::EsdQA(1);            } 
@@ -115,79 +116,6 @@ void mp  (Double_t probCut=0.7   ) {AliHMPIDTracker::MatrixPrint(probCut);}
 void stack(                     )   {AliHMPIDParam::Stack();}    
 void tid  (Int_t tid,Int_t evt=0)   {AliHMPIDParam::Stack(evt,tid);} 
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void tst()
-{
-  TCanvas *pC=new TCanvas("pads","View from electronics side, IP is behind the picture.",1000,900);  pC->ToggleEventStatus();
-  SimEvt();
-}//tst()
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void SimEvt()
-{
-//sample track         
-  Int_t ch=1;
-  Float_t th=8*TMath::DegToRad();                     
-  Float_t ph=gRandom->Rndm()*TMath::TwoPi(); 
-  Float_t radx=gRandom->Rndm()*AliHMPIDDigit::SizeAllX(); 
-  Float_t rady=gRandom->Rndm()*AliHMPIDDigit::SizeAllY();
-  Float_t simCkov=0.63,simErr=0; Int_t simN=8,recN=0;
-  
-  AliHMPIDRecon rec; rec.SetTrack(th,ph,radx,rady);                                                                
-//sample list of hits  
-  TClonesArray hl("AliHMPIDHit"); Int_t hc=0;  
-  TVector2 pc;//tmp position
-  
-  Int_t kCerenkov=50000050,kFeedback=50000051;
-  Int_t tn=0;
-                                                                                   new(hl[hc++]) AliHMPIDHit(ch,200e-9,kProton  ,tn++,radx  ,rady    );   //mip hit
-  for(int i=0;i<simN;i++){rec.TracePhot(simCkov,gRandom->Rndm()*TMath::TwoPi(),pc);new(hl[hc++]) AliHMPIDHit(ch,7.5e-9,kCerenkov,tn++,pc.X(),pc.Y());}  //photon hits
-  for(int i=0;i<10;i++)  { Float_t x=gRandom->Rndm()*130,y=gRandom->Rndm()*126;    new(hl[hc++]) AliHMPIDHit(ch,7.5e-9,kFeedback,tn++,   x  ,   y  );}  //bkg hits
-//do reconstruction
-                                                        TClonesArray sl("AliHMPIDDigit");              AliHMPIDv1::Hit2Sdi(&hl,&sl);                               
-  TObjArray dl(7);  for(Int_t i=0;i<7;i++) dl.AddAt(new TClonesArray("AliHMPIDDigit"),i);       AliHMPIDDigitizer::Sdi2Dig(&sl,&dl);     
-  TObjArray cl(7);  for(Int_t i=0;i<7;i++) cl.AddAt(new TClonesArray("AliHMPIDCluster"),i); AliHMPIDReconstructor::Dig2Clu(&dl,&cl);
-  
-  
-                             TClonesArray* pDigLst=(TClonesArray*)dl.At(ch);
-                             TClonesArray* pCluLst=(TClonesArray*)cl.At(ch);
-                    
-          Int_t recN=0;  Float_t recCkov=rec.CkovAngle(pCluLst,recN);    //reconstruct the ring
-  
-  AliESDtrack trk; trk.SetHMPIDtrk(radx,rady,th,ph);  trk.SetHMPIDmip(radx,rady,340,recN); trk.SetHMPIDsignal(recCkov); trk.SetHMPIDchi2(rec.CkovSigma2());
-  
-  Printf("Start of EVENT\n");                 
-  Printf("SDI------SDI---------SDI--------SDI------SDI------SDI");sl.Print();Printf("");
-  Printf("DIG------DIG---------DIG--------DIG------DIG------DIG");dl.Print();Printf("");                   
-  Printf("HIT------HIT---------HIT--------HIT------HIT------HIT");hl.Print();Printf("");
-  Printf("CLU------CLU---------CLU--------CLU------CLU------CLU");cl.Print();Printf("");                     
-  Printf("End of EVENT\n");                 
-  
-  
-  DrawEvt(&hl,pDigLst,pCluLst,&trk);  
-}//SimEvt()
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void DrawEvt(TClonesArray *pHitLst,TClonesArray *pDigLst,TClonesArray *pCluLst,AliESDtrack *pTrk)
-{
-//visualization  
-  TPolyLine    *pRing  =new TPolyLine;     pRing->SetLineColor(kMagenta);     
-  
-  Float_t recCkov=pTrk->GetHMPIDsignal();  Float_t recErr=TMath::Sqrt(pTrk->GetHMPIDchi2());
-  
-  Float_t trkTh,trkPh,trkX,trkY,mipX,mipY; pTrk->GetHMPIDtrk(trkX,trkY,trkTh,trkPh);  Int_t mipQ,recN; pTrk->GetHMPIDmip(mipX,mipY,mipQ,recN);
-  AliHMPIDRecon rec; rec.SetTrack(trkTh,trkPh,trkX,trkY);
-
-  TVector2 pos;    for(int j=0;j<100;j++){rec.TracePhot(recCkov,j*0.0628,pos); pRing->SetNextPoint(pos.X(),pos.Y());}  
-  
-  gPad->Clear(); TButton *pBtn=new TButton("Next","SimEvt()",0,0,0.07,0.05);   pBtn->Draw(); DrawPc(0);
-  
-  pRing->Draw();
-  pDigLst->Draw();
-  pCluLst->Draw();
-  pHitLst->Draw();  
-  TLatex txt; txt.SetTextSize(0.02);
-//  txt.DrawLatex(20,-5,Form("#theta=%.4f#pm%.5f with %2i #check{C}"          ,simCkov,simErr,simN));
-  txt.DrawLatex(25,-5,Form("#theta=%.4f#pm%.5f with %2i #check{C}"          ,recCkov,recErr,recN));
-//  txt.DrawLatex(0 ,127,Form("#theta=%.2f#circ   #phi=%.2f#circ @(%.2f,%.2f) ",th*TMath::RadToDeg(),ph*TMath::RadToDeg(),radx,rady));
-}//DrawEvt()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 void PrintMap()
 {
@@ -404,26 +332,6 @@ void HitQA(Double_t cut,Double_t cutele,Double_t cutR)
   gBenchmark->Show("HitsPlots");
 }//HitQA()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void RecWithStack()
-{
-  al->LoadHeader();al->LoadKinematics();
-  AliStack *pStk=al->Stack();  
-  
-  AliESD *pEsd=new AliESD;
-  for(Int_t iTrk=0;iTrk<pStk->GetNtrack();iTrk++){//stack loop
-    TParticle *pPart=pStk->Particle(iTrk);
-    if(pPart->GetPDG()->Charge()==0) continue; //neutral particles are not reconstructed
-    pEsd->AddTrack(new AliESDtrack(pPart));
-  }//stack loop
-  
-  pEsd->Print();
-  AliTracker::SetFieldMap(new AliMagF,1);
-  AliHMPIDTracker t; 
-  rl->LoadRecPoints(); 
-  t.LoadClusters(rl->TreeR()); 
-  t.PropagateBack(pEsd);
-  rl->UnloadRecPoints();
-}//RecWithStack()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 void CluQA(AliRunLoader *pAL)
 {
@@ -663,32 +571,120 @@ void DrawPc(Bool_t isFill)
   }//PC loop      
 }//DrawPc()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void ReadEvt() 
+void hed()
+{//event display from files
+  static TCanvas *pC=0;
+  static Int_t iEvt=0;
+  static Int_t iEvtTot=0;
+  static TFile *pEsdFl=0;
+  static TTree *pEsdTr=0;
+  static AliESD *pEsd=0;
+  
+  if(!pC){
+    if(hl==0) {Printf("hed: no HMPID loader");return;}
+    Printf("Opening session");
+    pEsdFl=TFile::Open("AliESDs.root");     if(!pEsdFl || !pEsdFl->IsOpen()) return;//open AliESDs.root
+    pEsdTr=(TTree*) pEsdFl->Get("esdTree"); if(!pEsdTr)                      return;//get ESD tree
+    pEsdTr->SetBranchAddress("ESD", &pEsd);
+    hl->LoadDigits(); hl->LoadRecPoints();
+    iEvtTot=pEsdTr->GetEntries();
+    pC=new TCanvas("hed","View from electronics side, IP is behind the picture.",1000,900);  pC->ToggleEventStatus(); pC->Divide(3,3);
+    pC->cd(7); TButton *pBtn=new TButton("Next","hed()",0,0,0.2,0.1);   pBtn->Draw(); 
+  }
+  if(iEvt<iEvtTot){
+    pEsdTr->GetEntry(iEvt); al->GetEvent(iEvt); hl->TreeD()->GetEntry(0); hl->TreeR()->GetEntry(0);
+    TLatex txt;   pC->cd(3);  txt.DrawLatex(0.2,0.2,Form("Event %i Total %i",iEvt,iEvtTot));
+    DrawEvt(pC,h->DigLst(),h->CluLst(),pEsd);
+    iEvt++;
+  }else{
+    Printf("Last event");
+    pC->Clear();
+  }
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void sed()
 {
-// Display digits, reconstructed tracks intersections and HMPID rings if available 
-// Arguments: none
-//   Returns: none    
-  TFile *pEsdFl=TFile::Open("AliESDs.root");     if(!pEsdFl || !pEsdFl->IsOpen()) return;//open AliESDs.root
-  TTree *pEsdTr=(TTree*) pEsdFl->Get("esdTree"); if(!pEsdTr)                      return;//get ESD tree
-                                                                 
-  AliESD *pEsd=0;  pEsdTr->SetBranchAddress("ESD", &pEsd);
+
+  static TCanvas *pC=0;
   
-  for(Int_t iEvt=0;iEvt<pEsdTr->GetEntries();iEvt++) {                //events loop
-    
-    
-    
-    pEsdTr->GetEntry(iEvt);                                       //get ESD for this event   
-    for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){//ESD tracks loop
-      AliESDtrack *pTrk = pEsd->GetTrack(iTrk);             //
-    }//ESD tracks loop
-    
-    al->GetEvent(iEvt);   rl->TreeD()->GetEntry(0); //get digits list
+  if(!pC){
+    pC=new TCanvas("hed","Simulated evets-View from electronics side, IP is behind the picture.",1000,900); pC->Divide(3,3);
+    pC->cd(7); TButton *pBtn=new TButton("Next","sed()",0,0,0.2,0.1);   pBtn->Draw(); 
+  }
+
+
+  
+  AliHMPIDRecon rec;  
+
+  TClonesArray lh("AliHMPIDHit"); 
+  TClonesArray ls("AliHMPIDDigit"); 
+  TObjArray    ld(7); for(Int_t i=0;i<7;i++) ld.AddAt(new TClonesArray("AliHMPIDDigit"),i);
+  TObjArray    lc(7); for(Int_t i=0;i<7;i++) lc.AddAt(new TClonesArray("AliHMPIDCluster"),i);
+  AliESD esd;
+  
+  
+  EsdFromStack(&esd);
+  HitsFromEsd(&esd,&lh);
+  
+
+
+
+             AliHMPIDv1::Hit2Sdi(&lh,&ls);                               
+      AliHMPIDDigitizer::Sdi2Dig(&ls,&ld);     
+  AliHMPIDReconstructor::Dig2Clu(&ld,&lc);
+//        AliHMPIDTracker::Recon(&esd,&cl);
+  
+  DrawEvt(pC,&ld,&lc,&esd);  
+}//SimEvt()
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void DrawEvt(TCanvas *pC,TObjArray *pDigLst,TObjArray *pCluLst,AliESD *pEsd)
+{//draws all the objects of current event
+
+  AliHMPIDRecon rec;  
+  TPolyMarker *pTxC[7];  TPolyMarker *pRin[7]; //intesections and rings
+  for(Int_t ch=0;ch<7;ch++){
+    pTxC[ch]=new TPolyMarker; pTxC[ch]->SetMarkerStyle(2); pTxC[ch]->SetMarkerColor(kRed); pTxC[ch]->SetMarkerSize(3);
+    pRin[ch]=new TPolyMarker; pRin[ch]->SetMarkerStyle(6); pRin[ch]->SetMarkerColor(kMagenta);
+  }
+  
+  
+  for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){//tracks loop to collect cerenkov rings and intersection points
+    AliESDtrack *pTrk=pEsd->GetTrack(iTrk);
+    Int_t ch=pTrk->GetHMPIDcluIdx();
+    if(ch<0) continue; //this track does not hit HMPID
+    ch/=1000000; 
+    Float_t th,ph,xPc,yPc; pTrk->GetHMPIDtrk(xPc,yPc,th,ph);  //get info on current track
+    pTxC[ch]->SetNextPoint(xPc,yPc);                          //add new intersection point
     
+    Float_t ckov=pTrk->GetHMPIDsignal();  Float_t err=TMath::Sqrt(pTrk->GetHMPIDchi2());
     
-  }//events loop
-  delete pEsd;  pEsdFl->Close();//close AliESDs.root
-//  rl->UnloadDigits();
-}//ReadEvt()
+    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());}      
+    }
+  }//tracks loop
+      
+  for(Int_t iCh=0;iCh<7;iCh++){//chambers loop    
+    switch(iCh){
+      case 6: pC->cd(1); break; case 5: pC->cd(2); break;
+      case 4: pC->cd(4); break; case 3: pC->cd(5); break; case 2: pC->cd(6); break;
+                                case 1: pC->cd(8); break; case 0: pC->cd(9); break;
+    }
+    gPad->SetEditable(kTRUE); gPad->Clear();
+    DrawPc(0);
+    ((TClonesArray*)pDigLst->At(iCh))->Draw();  //draw digits
+    ((TClonesArray*)pCluLst->At(iCh))->Draw();  //draw clusters
+                            pTxC[iCh]->Draw();  //draw intersections
+                            pRin[iCh]->Draw();  //draw rings
+    gPad->SetEditable(kFALSE);
+  }//chambers loop
+//  TLatex txt; txt.SetTextSize(0.02);
+//  txt.DrawLatex(20,-5,Form("#theta=%.4f#pm%.5f with %2i #check{C}"          ,simCkov,simErr,simN));
+//  txt.DrawLatex(25,-5,Form("#theta=%.4f#pm%.5f with %2i #check{C}"          ,recCkov,recErr,recN));
+//  txt.DrawLatex(0 ,127,Form("#theta=%.2f#circ   #phi=%.2f#circ @(%.2f,%.2f) ",th*TMath::RadToDeg(),ph*TMath::RadToDeg(),radx,rady));
+//  Printf("DIG------DIG---------DIG--------DIG------DIG------DIG");pDigLst->Print();Printf("");                   
+}//DrawEvt()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 void t1(Int_t case=1)
 {
@@ -722,7 +718,38 @@ void t1(Int_t case=1)
   cl->Print();  
 }//t1()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void EsdFromStack(AliESD *pEsd)
+{
+  al->LoadHeader();al->LoadKinematics();
+  AliStack *pStk=al->Stack();  
+  
+  for(Int_t iTrk=0;iTrk<pStk->GetNtrack();iTrk++){//stack loop
+    TParticle *pPart=pStk->Particle(iTrk);
+    if(pPart->GetPDG()->Charge()==0) continue; //neutral particles are not reconstructed
+    if(pPart->GetFirstMother()>0)    continue; //do not consider secondaries
+    AliESDtrack trk(pPart);
+    pEsd->AddTrack(&trk);
+  }//stack loop  
+}//EsdFromStack()
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void HitsFromEsd(AliESD *pEsd, TClonesArray *pHitLst)
+{
+  AliHMPIDRecon rec;
+  const Int_t kCerenkov=50000050,kFeedback=50000051;
+  Int_t hc=0; TVector2 pos;
+  for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){//tracks loop
+    AliESDtrack *pTrk=pEsd->GetTrack(iTrk);
+    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 th,ph,xPc,yPc,; pTrk->GetHMPIDtrk(xPc,yPc,th,ph);  
+                          new((*pHitLst)[hc++]) AliHMPIDHit(ch,200e-9,kProton  ,iTrk,xPc                ,yPc);                 //mip hit
+    for(int i=0;i<4;i++)  new((*pHitLst)[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++){
+      rec.TracePhot(ckov,gRandom->Rndm()*TMath::TwoPi(),pos);
+                          new((*pHitLst)[hc++]) AliHMPIDHit(ch,7.5e-9,kCerenkov,iTrk,pos.X(),pos.Y());}                      //photon hits  
+  }//tracks loop    
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++