HTA revisited algorithm.
authordibari <dibari@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 15 Sep 2008 16:18:17 +0000 (16:18 +0000)
committerdibari <dibari@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 15 Sep 2008 16:18:17 +0000 (16:18 +0000)
HMPID/AliHMPIDReconHTA.cxx
HMPID/AliHMPIDReconHTA.h
HMPID/AliHMPIDTracker.cxx
HMPID/AliHMPIDTracker.h
HMPID/HESDfromKin.C

index 68773e4..c291552 100644 (file)
 #include "AliHMPIDReconHTA.h"//class header
 #include "AliHMPIDCluster.h" //CkovHiddenTrk()
 #include "AliHMPIDRecon.h"   //FunMinPhot()
+#include <TFile.h>           //Database()
 #include <TMinuit.h>         //FitFree()
 #include <TClonesArray.h>    //CkovHiddenTrk()
 #include <AliESDtrack.h>     //CkovHiddenTrk()
-#include <TH2I.h>            //InitDatabase()
+#include <TH2F.h>            //InitDatabase()
 #include <TGraph.h>          //ShapeModel()
 #include <TSpline.h>         //ShapeModel()
 #include "TStopwatch.h"      //
 
-TH2I* AliHMPIDReconHTA::fgDatabase = 0x0;
+TH2F* AliHMPIDReconHTA::fgDatabase = new TH2F("deconv","database;d1;d2;thC+1000*thTrk",500,0,50,150,0,15);
+
 
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 AliHMPIDReconHTA::AliHMPIDReconHTA():
@@ -57,7 +59,7 @@ AliHMPIDReconHTA::AliHMPIDReconHTA():
 //hidden algorithm
 //..
   fParam->SetRefIdx(fParam->MeanIdxRad()); // initialization of ref index to a default one
-  if(!fgDatabase) InitDatabase();
+  if(fgDatabase->GetEntries()<1) InitDatabase();
 }
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 AliHMPIDReconHTA::~AliHMPIDReconHTA()
@@ -87,60 +89,61 @@ void AliHMPIDReconHTA::DeleteVars()const
   if(fClCk) delete fClCk;
 }
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Bool_t AliHMPIDReconHTA::CkovHiddenTrk(AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t nmean, Double_t qthre)
+Bool_t AliHMPIDReconHTA::CkovHiddenTrk(AliESDtrack *pTrk,TClonesArray *pCluLst,Int_t index, Double_t nmean)
 {
 // Pattern recognition method without any infos from tracking:HTA (Hidden Track Algorithm)...
 // The method finds in the chmuber the cluster with the highest charge
 // compatibile with a MIP, then the strategy is applied
 // Arguments:  pTrk     - pointer to ESD track
 //             pCluLs   - list of clusters for a given chamber 
-//             nmean    - mean freon ref. index
+//             pNmean   - pointer to ref. index
+//             pQthre   - pointer to qthre
 //   Returns:           - 0=ok,1=not fitted 
   
   AliHMPIDParam *pParam = AliHMPIDParam::Instance(); 
-  pParam->SetRefIdx(nmean);
 
-  if(!CluPreFilter(pCluLst)) {return kFALSE;}
-  
-  Float_t mipX=-1,mipY=-1;Int_t mipId=-1,mipQ=-1;                                                                           
-  Double_t qRef = 0;
+  if(!CluPreFilter(pCluLst)) return kFALSE;
+
   Int_t nCh=0;
   Int_t sizeClu=0;
-  for (Int_t iClu=0;iClu<pCluLst->GetEntriesFast();iClu++){                                   //clusters loop
+  
+  fNClu = pCluLst->GetEntriesFast();
+    
+  for (Int_t iClu=0;iClu<fNClu;iClu++){                                                       //clusters loop
     AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu);                       //get pointer to current cluster    
-    nCh = pClu->Ch();
     fXClu[iClu] = pClu->X();fYClu[iClu] = pClu->Y();                                          //store x,y for fitting procedure
     fClCk[iClu] = kTRUE;                                                                      //all cluster are accepted at this stage to be reconstructed
-    if(pClu->Q()>qthre) fClCk[iClu] = kFALSE;                                                 // not a good photon in any case (multiple MIPs)
-    if(pClu->Q()>qRef){                                                                       //searching the highest charge to select a MIP      
-      qRef = pClu->Q();
-      mipId=iClu; mipX=pClu->X();mipY=pClu->Y();mipQ=(Int_t)pClu->Q();
-      sizeClu=pClu->Size();
-    }                                                                                    
-//    Printf(" n. %d x %f y %f Q %f",iClu,pClu->X(),pClu->Y(),pClu->Q());
+    
+    if(iClu == index) {
+
+      fMipX = pClu->X();
+      fMipY = pClu->Y();
+      fMipQ = pClu->Q();
+      sizeClu = pClu->Size();
+      nCh = pClu->Ch();
+      fClCk[index] = kFALSE;
+      fIdxMip = index;
+ //    Printf(" n. %d x %f y %f Q %f",iClu,pClu->X(),pClu->Y(),pClu->Q());
+    }
   }//clusters loop
-
-  fNClu = pCluLst->GetEntriesFast();
-  if(qRef>qthre){                                                                     //charge compartible with MIP clusters
-    fIdxMip = mipId;
-    fClCk[mipId] = kFALSE;
-    fMipX = mipX; fMipY=mipY; fMipQ = qRef;
-//    Printf(" mipId %d x %f y %f Q %f",fIdxMip,fMipX,fMipY,fMipQ);
-    if(!DoRecHiddenTrk()) {
-      pTrk->SetHMPIDsignal(pParam->kNoPhotAccept);
-      return kFALSE;
-    }                                                                           //Do track and ring reconstruction,if problems returns 1
-    pTrk->SetHMPIDtrk(fRadX,fRadY,fThTrkFit,fPhTrkFit);                                        //store track intersection info
-    pTrk->SetHMPIDmip(fMipX,fMipY,(Int_t)fMipQ,fNClu);                                         //store mip info 
-    pTrk->SetHMPIDcluIdx(nCh,fIdxMip+1000*sizeClu);                                            //set cham number, index of cluster + cluster size
-    pTrk->SetHMPIDsignal(fCkovFit);                                                            //find best Theta ckov for ring i.e. track
-    pTrk->SetHMPIDchi2(fCkovSig2);                                                             //errors squared
+  
+  pParam->SetRefIdx(nmean);
+  
+  if(!DoRecHiddenTrk()) {
+    pTrk->SetHMPIDsignal(pParam->kNoPhotAccept);
+    return kFALSE;
+  }                                                                           //Do track and ring reconstruction,if problems returns 1
+  
+  pTrk->SetHMPIDtrk(fRadX,fRadY,fThTrkFit,fPhTrkFit);                                        //store track intersection info
+  pTrk->SetHMPIDmip(fMipX,fMipY,(Int_t)fMipQ,fNClu);                                         //store mip info 
+  pTrk->SetHMPIDcluIdx(nCh,fIdxMip+1000*sizeClu);                                            //set cham number, index of cluster + cluster size
+  pTrk->SetHMPIDsignal(fCkovFit);                                                            //find best Theta ckov for ring i.e. track
+  pTrk->SetHMPIDchi2(fCkovSig2);                                                             //errors squared
 //    Printf(" n clusters tot %i accepted %i",pCluLst->GetEntriesFast(),fNClu);
 //    Printf("CkovHiddenTrk: thetaC %f th %f ph %f",fCkovFit,fThTrkFit,fPhTrkFit);
-    return kTRUE;
-  }
   
-  return kFALSE;
+  return kTRUE;
+  
 }//CkovHiddenTrk()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 Bool_t AliHMPIDReconHTA::DoRecHiddenTrk()
@@ -224,10 +227,17 @@ Bool_t AliHMPIDReconHTA::FindShape(Double_t &thTrkRec,Double_t &phiTrkRec,Double
 
   Bool_t status;
     
-// Sort in phi angle...  
+// Sort in phi angle...
+//  Printf(" mipX %f mipy %f",fMipX,fMipY);
   for(Int_t i=0;i<fNClu;i++) {
+    if(!fClCk[i]) {
+      phiphot[i] = 999.;
+      dist[i] = 999.;
+      continue;
+    }
     phiphot[i] = (TMath::ATan2(fMipY-fYClu[i],fMipX-fXClu[i])+TMath::Pi())*TMath::RadToDeg();
     dist[i]=TMath::Sqrt((fMipX-fXClu[i])*(fMipX-fXClu[i])+(fMipY-fYClu[i])*(fMipY-fYClu[i]));
+//    Printf(" n.%3i  phiphot %f dist %f check %i",i,phiphot[i],dist[i],fClCk[i]);
   }
   
   TMath::Sort(fNClu,phiphot,indphi,kFALSE);
@@ -265,6 +275,7 @@ Bool_t AliHMPIDReconHTA::FindShape(Double_t &thTrkRec,Double_t &phiTrkRec,Double
     if(!fClCk[indphi[i]]) continue;                                                  // Check if a good photon candidate or not
     phiphotP[npeff] = phiphot[indphi[i]];
     distP[npeff]    = dist[indphi[i]];
+//    Printf("n. %2i phi %f dist %f",npeff,phiphotP[npeff],distP[npeff]);
     npeff++;
   }
   
@@ -281,26 +292,26 @@ Bool_t AliHMPIDReconHTA::FindShape(Double_t &thTrkRec,Double_t &phiTrkRec,Double
 //  for(Int_t i=0;i<npeff;i++) {Printf(" n. %d phiphot %f dist %f",i,phiphotP[i],distP[i]);}
   
   Double_t xA,xB;
+  status = kFALSE;
   if(ShapeModel(npeff,phiphotP,distP,xA,xB,phiTrkRec)) {
     
-//    Printf("FindShape: phi start %f",phiTrkRec*TMath::RadToDeg());
+//    Printf("FindShape: phi start %f xA %f yA %f",phiTrkRec*TMath::RadToDeg(),xA,xB);
+    if(xA < 50 && xB < 15) {                     // limits of the Database. See TH2F in InitDatabase()
 
-    Int_t bin = fgDatabase->FindBin(xA,xB);
-    Int_t compact = (Int_t)fgDatabase->GetBinContent(bin);
-    thetaCRec =  (Double_t)(compact%1000);
-    thTrkRec  =  (Double_t)(compact/1000);
+      Int_t bin = fgDatabase->FindBin(xA,xB);
+      if(bin>0) {
+        Int_t compact = (Int_t)fgDatabase->GetBinContent(bin);
+        thetaCRec =  (Double_t)(compact%1000);
+        thTrkRec  =  (Double_t)(compact/1000);
 
-    thTrkRec *= TMath::DegToRad(); 
-    thetaCRec *= TMath::DegToRad();
+        thTrkRec *= TMath::DegToRad(); 
+        thetaCRec *= TMath::DegToRad();
 
-//    Printf("FindShape: xA %f xB %f compact %d thTrk %f thC %f",xA,xB,compact,thTrkRec*TMath::RadToDeg(),thetaCRec*TMath::RadToDeg());
-        
-    status = kTRUE;
-    
-  } else {
-    
-    status = kFALSE;
-    
+    //    Printf("FindShape: xA %f xB %f compact %d thTrk %f thC %f",xA,xB,compact,thTrkRec*TMath::RadToDeg(),thetaCRec*TMath::RadToDeg());
+
+        status = kTRUE;
+      }
+    }
   }
 
   delete [] phiphotP;
@@ -350,8 +361,6 @@ Bool_t AliHMPIDReconHTA::ShapeModel(Int_t np,Double_t *phiphot,Double_t *dist,Do
                              phiphot[locMax+ip[2]],dist[locMax+ip[2]]);
   if(maxXf< phiphot[locMax+ip[0]] || maxXf > phiphot[locMax+ip[2]]) maxXf = maxX;
   
-//  Printf(" phi at mindist %f and found %f",minX,minXf);
-//  Printf(" phi at maxdist %f and found %f",maxX,maxXf);
 //  
   if(TMath::Abs(maxXf-minXf)>30) {
     xA = sphi->Eval(minXf);
@@ -370,6 +379,7 @@ Bool_t AliHMPIDReconHTA::ShapeModel(Int_t np,Double_t *phiphot,Double_t *dist,Do
 //  Printf("ShapeModel: phiStart %f xA %f xB %f",phiStart,xA,xB);
   
   phiStart*=TMath::DegToRad();
   return kTRUE;
 }
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -510,6 +520,8 @@ void AliHMPIDReconHTA::InitDatabase()
 //        The content is the packed info of track theta and thetaC in degrees
 //                        thetaC+1000*thTrk
 //
+//  TFile *pout = new TFile("./database.root","recreate");
+  
   TStopwatch timer;
   timer.Start();
   
@@ -526,7 +538,7 @@ void AliHMPIDReconHTA::InitDatabase()
   Int_t nstepx = 1000;
   Int_t nstepy = 1000;
 
-  TH2I *deconv = new TH2I("deconv","database;d1;d2;thC+1000*thTrk",500,0,50,150,0,15);
+//  TH2F *fgDatabase = new TH2F("deconv","database;d1;d2;thC+1000*thTrk",500,0,50,150,0,15);
   //
   Double_t xrad = 0;
   Double_t yrad = 0;
@@ -581,31 +593,35 @@ void AliHMPIDReconHTA::InitDatabase()
       Double_t distB   = TMath::Sqrt((x[2]-xmip)*(x[2]-xmip)+(y[2]-ymip)*(y[2]-ymip));
 // compact the infos...      
       Int_t compact = (Int_t)(thetaC*TMath::RadToDeg())+1000*(Int_t)(thTrk*TMath::RadToDeg());
-      Int_t bin = deconv->FindBin(distA,distB);
-      if(deconv->GetBinContent(bin)==0) deconv->Fill(distA,distB,compact);
+      Int_t bin = fgDatabase->FindBin(distA,distB);
+      if(fgDatabase->GetBinContent(bin)==0) fgDatabase->Fill(distA,distB,compact);
     }
   }
 
-  FillZeroChan(deconv);
-  fgDatabase = deconv;
+  FillZeroChan();
+//  fgDatabase = deconv;
 
   timer.Stop();
   Double_t nSecs = timer.CpuTime();  
   AliInfo(Form("database HTA successfully open in %3.1f sec.(CPU). Reconstruction is started.",nSecs));
+  
+//  pout->Write();
+//  pout->Close();
+  
 }
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPIDReconHTA::FillZeroChan(TH2I *deconv)const
+void AliHMPIDReconHTA::FillZeroChan()const
 {
   //If fills eventually channel without entries
   //inthe histo "database" jyst interpolating the neighboring cells
   // Arguments: histogram pointer of the database
   //   Returns: none
   //
-  Int_t nbinx = deconv->GetNbinsX();
-  Int_t nbiny = deconv->GetNbinsY();
+  Int_t nbinx = fgDatabase->GetNbinsX();
+  Int_t nbiny = fgDatabase->GetNbinsY();
   for(Int_t i = 0;i<nbinx;i++) {
     for(Int_t j = 0;j<nbiny;j++) {
-      if(deconv->GetBinContent(i,j) == 0) {
+      if(fgDatabase->GetBinContent(i,j) == 0) {
         Int_t nXmin = i-1; Int_t nXmax=i+1;
         Int_t nYmin = j-1; Int_t nYmax=j+1;
         Int_t nc = 0;
@@ -615,13 +631,13 @@ void AliHMPIDReconHTA::FillZeroChan(TH2I *deconv)const
           if(ix<0||ix>nbinx) continue;
           for(Int_t iy=nYmin;iy<=nYmax;iy++) {
             if(iy<0||iy>nbiny) continue;
-            meanC  +=  (Int_t)deconv->GetBinContent(ix,iy)%1000;
-            meanTrk+=  (Int_t)deconv->GetBinContent(ix,iy)/1000;
+            meanC  +=  (Int_t)fgDatabase->GetBinContent(ix,iy)%1000;
+            meanTrk+=  (Int_t)fgDatabase->GetBinContent(ix,iy)/1000;
             nc++;
           }
           meanC/=nc; meanTrk/=nc;
           Int_t compact = (Int_t)meanC+1000*(Int_t)meanTrk;
-          if(compact>0)deconv->SetCellContent(i,j,compact);
+          if(compact>0)fgDatabase->SetCellContent(i,j,compact);
         }
       }
     }
index 6b673b6..f7d86b5 100644 (file)
 
 
 #include <TTask.h>        //base class
+#include <TH2F.h>         //InitDatabase()
 
 class TClonesArray; //CkovAngle()
 class AliESDtrack;  //CkovAngle()
-class TH2I;         //InitDatabase()
 class AliHMPIDParam;//General pourpose
 
 class AliHMPIDReconHTA : public TTask 
@@ -30,9 +30,9 @@ public :
   void     InitVars         (Int_t n);                                                             //init space for variables
   void     DeleteVars       ()const;                                                               //delete variables
   void     InitDatabase     ();                                                                    //initialization of database
-  TH2I*    DBHTA            ()     {return fgDatabase;}                                            //pointer for HTA database of rings
-  void     FillZeroChan     (TH2I *deconv)const;                                                   //complete the DB
-  Bool_t   CkovHiddenTrk    (AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t nmean,Double_t qthre);//Pattern recognition without trackinf information
+  TH2F*    DBHTA            ()     {return fgDatabase;}                                            //pointer for HTA database of rings
+  void     FillZeroChan     ()const;                                                               //complete the DB
+  Bool_t   CkovHiddenTrk    (AliESDtrack *pTrk,TClonesArray *pClu,Int_t index, Double_t nmean);    //Pattern recognition without trackinf information
   Bool_t   CluPreFilter     (TClonesArray *pClu               );                                   //Pre clustering filter to cut bkg clusters
   Bool_t   DoRecHiddenTrk   (                                 );                                   //Calling to the fitted procedures
   Bool_t   FindShape        (Double_t &thTrkRec,Double_t &phiTrkRec,Double_t &thetaCRec);          //Find shape of the ring
@@ -78,13 +78,13 @@ protected:
   Double_t fCkovSig2;                          //estimated error^2 on ring Cherenkov angle
   
   AliHMPIDParam *fParam;                       //Pointer to AliHMPIDParam
-  static TH2I* fgDatabase;                     //database for ring shapes
+  static TH2F* fgDatabase;                     //database for ring shapes
 //
 private:
   AliHMPIDReconHTA(const AliHMPIDReconHTA& r);              //dummy copy constructor
   AliHMPIDReconHTA &operator=(const AliHMPIDReconHTA& r);   //dummy assignment operator
 //
-  ClassDef(AliHMPIDReconHTA,0)
+  ClassDef(AliHMPIDReconHTA,1)
 };
 
 #endif // #ifdef AliHMPIDReconHTA_cxx
index e388ef5..4777e4b 100644 (file)
@@ -260,20 +260,95 @@ Int_t AliHMPIDTracker::Recon(AliESDEvent *pEsd,TObjArray *pClus,TObjArray *pNmea
   return 0; // error code: 0=no error;
 }//Recon()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Int_t AliHMPIDTracker::ReconHiddenTrk(Int_t iCh,Int_t iHVsec,AliESDtrack *pTrk,TClonesArray *pCluLst,TObjArray *pNmean,TObjArray *pQthre)
+Int_t AliHMPIDTracker::ReconHiddenTrk(AliESDEvent *pEsd,TObjArray *pClus,TObjArray *pNmean, TObjArray *pQthre)
 {
 // Static method to reconstruct Theta Ckov for all valid tracks of a given event.
 // Arguments: pEsd- pointer ESD; pClu- pointer to clusters for all chambers; pNmean - pointer to all function Nmean=f(time), pQthre - pointer to all function Qthre=f(time)
 //   Returns: error code, 0 if no errors
-  AliHMPIDReconHTA reconHTA;                                                                          //instance of reconstruction class, nothing important in ctor
-  Double_t nmean=((TF1*)pNmean->At(3*iCh))->Eval(0);                                            //C6F14 Nmean for this chamber
-  Double_t qthre = 0;
-  if(pQthre->GetEntriesFast()==AliHMPIDParam::kMaxCh+1)                                         //
-    qthre=((TF1*)pQthre->At(iCh))->Eval(0);                                                     //just for backward compatibi
-  else  qthre=((TF1*)pQthre->At(6*iCh+iHVsec))->Eval(0);                                        //
-  if(pCluLst->GetEntriesFast()<4) return 1;                                                     //min 4 clusters (3 + 1 mip) to find a ring! 
-  if(reconHTA.CkovHiddenTrk(pTrk,pCluLst,nmean,qthre)) return 0;                                   //search for track parameters and Cerenkov angle of this track
-  else return 1;                                                                                // error code: 0=no error,1=fit not performed;
+  
+  AliHMPIDReconHTA reconHTA;                                                                     //instance of reconstruction class, nothing important in ctor
+  
+  AliHMPIDParam *pParam = AliHMPIDParam::Instance();                                             //Instance of AliHMPIDParam
+  
+  for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){                                        //loop on the ESD tracks in the event
+    
+    AliESDtrack *pTrk = pEsd->GetTrack(iTrk);                                                     //here it is simulated or just empty track
+    Int_t ipCh;
+    ipCh = pTrk->GetHMPIDcluIdx();ipCh/=1000000;
+    if(ipCh<0) continue;
+
+    TClonesArray *pMipCluLst=(TClonesArray *)pClus->At(ipCh);                                   //get the list of clusters
+    Int_t nMipClusTot = pMipCluLst->GetEntries();                                               //total number of clusters in the given chamber
+    
+    Double_t qMip=-1;
+    Int_t chMip=-1;    
+    Double_t xMip,yMip;
+    Int_t indMip;
+    Int_t cluMipSiz;
+
+    for (Int_t iClu=0; iClu<nMipClusTot;iClu++) {                                               //clusters loop
+      
+      AliHMPIDCluster *pClu=(AliHMPIDCluster*)pMipCluLst->UncheckedAt(iClu);                    //get the cluster
+      Double_t qClus = pClu->Q();
+      if(qClus>qMip) {
+        qMip  = qClus;
+        chMip = pClu->Ch();
+        xMip = pClu->X();
+        yMip = pClu->Y();
+        indMip = iClu;
+        cluMipSiz = pClu->Size();
+      }
+    }//clus loop
+
+    if(chMip<0) return 1;    
+    
+    Int_t hvsec;
+    Double_t qthre;
+// evaluate qThre
+    if(pQthre->GetEntriesFast()==pParam->kMaxCh+1) {                                              // just for backward compatibility
+      qthre=((TF1*)pQthre->At(chMip))->Eval(pEsd->GetTimeStamp());                                //
+    } else {                                                                                      // in the past just 1 qthre
+      hvsec = pParam->InHVSector(yMip);                                                           //  per chamber
+      if(hvsec>=0) qthre=((TF1*)pQthre->At(6*chMip+hvsec))->Eval(pEsd->GetTimeStamp());           //
+    }
+//
+    if(qMip<qthre) {
+      pTrk->SetHMPIDmip(xMip,yMip,(Int_t)qMip,0);                                                        //store mip info in any case 
+      pTrk->SetHMPIDcluIdx(chMip,indMip+1000*cluMipSiz);                                                          
+      pTrk->SetHMPIDsignal(pParam->kMipQdcCut);
+      return 1;                                                                                   //charge compatible with MIP clusters      
+    }
+      
+    pTrk->SetHMPIDmip(xMip,yMip,(Int_t)qMip,0);                                                   //store mip info in any case 
+    pTrk->SetHMPIDcluIdx(chMip,indMip+1000*cluMipSiz);                                                          
+
+    Double_t yRa = yMip;                                                                        //just an approx...
+    Double_t nmean;
+    //evaluate nMean
+    if(pNmean->GetEntries()==21) {                                                              //for backward compatibility
+      nmean=((TF1*)pNmean->At(3*chMip))->Eval(pEsd->GetTimeStamp());                            //C6F14 Nmean for this chamber
+    } else {
+      Int_t iRad     = pParam->Radiator(yRa);                                                   //evaluate the radiator involved
+      if(iRad < 0) {
+       nmean = -1;
+      } else {
+      Double_t tLow  = ((TF1*)pNmean->At(6*chMip+2*iRad  ))->Eval(pEsd->GetTimeStamp());        //C6F14 low  temp for this chamber
+      Double_t tHigh = ((TF1*)pNmean->At(6*chMip+2*iRad+1))->Eval(pEsd->GetTimeStamp());        //C6F14 high temp for this chamber
+      Double_t tExp  = pParam->FindTemp(tLow,tHigh,yRa);                                        //estimated temp for that chamber at that y
+      nmean = pParam->NIdxRad(AliHMPIDParam::Instance()->GetEPhotMean(),tExp);                  //mean ref idx @ a given temp
+      }
+      if(nmean < 0){                                                                            //track didn' t pass through the radiator
+         pTrk->SetHMPIDsignal(AliHMPIDRecon::kNoRad);                                           //set the appropriate flag
+         return 1;
+      }
+    }
+    //
+//    Printf(" qthre %f nmean %f index %i cham %i",qthre,nmean,indMip,chMip);
+    reconHTA.CkovHiddenTrk(pTrk,(TClonesArray *)pClus->At(ipCh),indMip,nmean);                  //search for track parameters and Cerenkov angle of this track
+  }//iTrk
+
+  return 0; // error code: 0=no error;
+
 }//Recon()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 void AliHMPIDTracker::FillClusterArray(TObjArray* array) const {
index b95f4a8..ecb6130 100644 (file)
@@ -32,7 +32,7 @@ public:
   static Int_t       IntTrkCha     (Int_t ch,AliHMPIDtrack *pTrk,Float_t &xPc,Float_t &yPc,Float_t &xRa,Float_t &yRa,Float_t &theta,Float_t &phi);//find track-PC intersection, retuns chamber ID
 
   static Int_t       Recon         (AliESDEvent *pEsd,TObjArray *pCluAll,TObjArray *pNmean=0,TObjArray *pQthre=0);//do actual job, returns status code  
-  static Int_t       ReconHiddenTrk(Int_t iCh,Int_t iHVsec,AliESDtrack *pTrk,TClonesArray *pClus,TObjArray *pNmean, TObjArray *pQthre);//do actual job with Hidden Track Algorithm    
+  static Int_t       ReconHiddenTrk(AliESDEvent *pEsd,TObjArray *pClus,TObjArray *pNmean, TObjArray *pQthre);//do actual job with Hidden Track Algorithm    
   
   
 protected:
index 14c9b6c..e057ccf 100644 (file)
@@ -28,32 +28,44 @@ void HESDfromKin(const char *name="default")
     
   TString ttl=name;
   Bool_t htaCheck=ttl.Contains("HTA");
-  if(!htaCheck) SimEsd(pHL,pEsd); else SimEsdHidden(pHL,pEsd);
+//  if(!htaCheck) SimEsd(pHL,pEsd); else SimEsdHidden(pHL,pEsd);
+  SimEsd(pHL,pEsd,htaCheck);
   
   pEsdFl->cd();
   pEsdFl->Write();pEsdFl->Close();        
 }
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void SimEsd(AliLoader *pHL,AliESDEvent *pEsd)
+void SimEsd(AliLoader *pHL,AliESDEvent *pEsd,Bool_t htaCheck)
 {
-  Printf("-----------------------------------------------");
-  Printf("| SimESD: Utility to embed ESD from kinematics|");
-  Printf("-----------------------------------------------");
+  if(htaCheck) {
+    TFile *fout = new TFile("HTA.root","recreate");
+    TH1F *hdC   = new TH1F("dC"  ,";delta Cerenkov (rad)",100,-0.2,0.2);
+    TH1F *hCer  = new TH1F("Cer" ,"Theta Cerenkov (rad)",250,0.,0.75);
+    TH2F *htvsp = new TH2F("tvsp",";momentum (GeV/c);theta Cerenkov (rad)",100,0.,5.,1000,0.,0.75);
+    TH1F *hdth  = new TH1F("dth" ,";Delta theta Trk (mrad)",100,-250,250);
+    TH1F *hdph  = new TH1F("dph" ,";Delta phi Trk (mrad)",100,-500,500);
+    Double_t rd=TMath::RadToDeg();
+    Printf("----------------------------------------------");
+    Printf("| SimHTA:Utility to embed ESD from kinematics|");
+    Printf("|     with  Hidden Track Algorithm (HTA)     |");
+    Printf("----------------------------------------------");
+  } else {
+    Printf("-----------------------------------------------");
+    Printf("| SimESD: Utility to embed ESD from kinematics|");
+    Printf("-----------------------------------------------");
+}
   AliHMPIDTracker::SetFieldMap(gAL->GetAliRun()->Field(),kTRUE);
   AliHMPID *pH=(AliHMPID*)gAL->GetAliRun()->GetDetector("HMPID");
-  Int_t mtid=-1;
   Int_t iNevt=gAL->GetNumberOfEvents();
   pEsd->SetMagneticField(AliHMPIDTracker::GetBz());
-  Printf("Number of events to process: %i",iNevt);
   for(Int_t iEvt=0;iEvt<iNevt;iEvt++){//events loop
-    if(!(iEvt%50)) Printf("Events processed %i",iEvt);
     gAL->GetEvent(iEvt);    
     pHL->TreeR()->GetEntry(0);
     AliStack *pStack=gAL->Stack();
+    Int_t nTrkHMPID=0;
     for(Int_t i=0;i<pStack->GetNtrack();i++){
+      if(!pStack->IsPhysicalPrimary(i)) continue;
       TParticle *pTrack=pStack->Particle(i); 
-      mtid=pTrack->GetFirstMother();
-      if(mtid>=0) continue; // only primaries
       if(pTrack->GetPDG()->Charge()==0) continue;
       AliESDtrack trk(pTrack); 
       Float_t xPc,yPc,xRa,yRa,thRa,phRa;
@@ -64,17 +76,48 @@ void SimEsd(AliLoader *pHL,AliESDEvent *pEsd)
         trk.SetHMPIDsignal(AliHMPIDRecon::kNotPerformed);                                        //ring reconstruction not yet performed
         continue;                                                           //no intersection at all, go after next track
       }
-//      Printf(" particle analyzed %d with mtid %d is %s hitting chamber %d",i,mtid,pTrack->GetName(),iCh);
+      nTrkHMPID++;
       trk.SetHMPIDcluIdx   (iCh,99999);                                                          //chamber not found, mip not yet considered
       trk.SetHMPIDtrk(xRa,yRa,thRa,phRa);                                                        //store initial infos
       pEsd->AddTrack(&trk);
     }// track loop
     
-    AliHMPIDTracker::Recon(pEsd,pH->CluLst(),pNmean,pQthre);
+    if(!(iEvt%50)) Printf("Number of events processed: %i with tracks %i in HMPID",iEvt,nTrkHMPID);
+    if(!htaCheck) AliHMPIDTracker::Recon(pEsd,pH->CluLst(),pNmean,pQthre);
 // Perform PID
     
+    
     for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){                                       //ESD tracks loop
+      
       AliESDtrack *pTrk = pEsd->GetTrack(iTrk);                                                    //get reconstructed track    
+      
+      Float_t radX,radY,thetaTrk,phiTrk;
+      pTrk->GetHMPIDtrk(radX,radY,thetaTrk,phiTrk);
+      
+      if(htaCheck) {
+          Int_t iCh = pTrk->GetHMPIDcluIdx();                                                                  //chamber not found, mip not yet considered
+          iCh/=1000000;
+  //      Printf("simulated track theta %f phi %f",thetaTrk*rd,phiTrk*rd);
+        TObjArray *pClus = pH->CluLst();
+
+        if(AliHMPIDTracker::ReconHiddenTrk(pEsd,pClus,pNmean,pQthre)!=0) continue;
+
+        Double_t thetaCerSim = TMath::ACos(pTrack->Energy()/(1.292*pTrack->P()));
+        Printf("-----------------------------------------------------------");
+        Printf(" theta Cerenkov simulated     %f",thetaCerSim);
+        Printf(" theta Cerenkov reconstructed %f",pTrk->GetHMPIDsignal());
+        Float_t thetaFit,phiFit;
+        pTrk->GetHMPIDtrk(radX,radY,thetaFit,phiFit);
+//        Printf("reconstr. track theta %f phi %f",thetaFit*rd,phiFit*rd);
+
+        // fill useful histos
+        hdC->Fill(pTrk->GetHMPIDsignal()-thetaCerSim);
+        hCer->Fill(pTrk->GetHMPIDsignal());
+        htvsp->Fill(pTrk->P(),pTrk->GetHMPIDsignal());
+        hdth->Fill((thetaFit-thetaTrk)*1000);
+        hdph->Fill((phiFit-phiTrk)*1000);
+      }
+      
       AliHMPIDPid pID;
       Double_t prob[5];
       pID.FindPid(pTrk,5,prob);
@@ -105,25 +148,25 @@ void SimEsdHidden(AliLoader *pHL,AliESDEvent *pEsd)
   Printf("----------------------------------------------");
   AliHMPIDTracker::SetFieldMap(gAL->GetAliRun()->Field(),kTRUE);
   AliHMPID *pH=(AliHMPID*)gAL->GetAliRun()->GetDetector("HMPID");
-  Int_t mtid=-1;
   Int_t iNevt=gAL->GetNumberOfEvents();
+  pEsd->SetMagneticField(AliHMPIDTracker::GetBz());
   
   Printf("Number of events to process: %i",iNevt);
   
   for(Int_t iEvt=0;iEvt<iNevt;iEvt++){//events loop
-    if(!(iEvt%50)) Printf("Events processed %i",iEvt);
+    if(!(iEvt%1)) Printf("Events processed %i",iEvt);
     gAL->GetEvent(iEvt);    
     pHL->TreeR()->GetEntry(0);
     AliStack *pStack=gAL->Stack();
     
     for(Int_t i=0;i<pStack->GetNtrack();i++){
       
+      (!pStack->IsPhysicalPrimary(i)) continue;
       TParticle *pTrack=pStack->Particle(i); 
-      mtid=pTrack->GetFirstMother();
-      
-      if(mtid>=0) continue; // only primaries
       
+      if(pTrack->GetPDG()->Charge()==0) continue;
       //find the chamber that intersects HMPID
+      pTrack->Print();
       AliESDtrack trk(pTrack);
       Float_t xPc,yPc,xRa,yRa,thRa,phRa;
       Int_t iCh=AliHMPIDTracker::IntTrkCha(&trk,xPc,yPc,xRa,yRa,thRa,phRa);         //get chamber intersected by this track