Still optimizing HTA
authordibari <dibari@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 18 Sep 2008 18:06:33 +0000 (18:06 +0000)
committerdibari <dibari@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 18 Sep 2008 18:06:33 +0000 (18:06 +0000)
HMPID/AliHMPIDReconHTA.cxx
HMPID/AliHMPIDReconHTA.h
HMPID/AliHMPIDTracker.cxx
HMPID/AliHMPIDTracker.h
HMPID/HESDfromKin.C

index 4cd7457..8959f68 100644 (file)
@@ -129,10 +129,16 @@ Bool_t AliHMPIDReconHTA::CkovHiddenTrk(AliESDtrack *pTrk,TClonesArray *pCluLst,I
   
   pParam->SetRefIdx(nmean);
   
+  //
+  Float_t xra,yra,th,ph; pTrk->GetHMPIDtrk(xra,yra,th,ph);
+  Printf(" simulated phi %6.2f ",ph*TMath::RadToDeg());
+  //
+  
   if(!DoRecHiddenTrk()) {
     pTrk->SetHMPIDsignal(pParam->kNoPhotAccept);
     return kFALSE;
   }                                                                           //Do track and ring reconstruction,if problems returns 1
+  Printf("    fitted phi %6.2f ",fPhTrkFit*TMath::RadToDeg());
   
   pTrk->SetHMPIDtrk(fRadX,fRadY,fThTrkFit,fPhTrkFit);                                        //store track intersection info
   pTrk->SetHMPIDmip(fMipX,fMipY,(Int_t)fMipQ,fNClu);                                         //store mip info 
@@ -347,11 +353,14 @@ 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();
   
+  phiStart*=TMath::DegToRad();
   Double_t phitest = FindSimmPhi(); 
-
-  Printf("phiStart %f phiTest %f",phiStart*TMath::RadToDeg(),phitest*TMath::RadToDeg());
+    //Double_t phiStart = phitest;
+  phiStart = phitest;  
+  Printf("   started phi %6.2f ",phiStart*TMath::RadToDeg());
+  
+//  Printf("phiStart %f phiTest %f",phiStart*TMath::RadToDeg(),phitest*TMath::RadToDeg());
  
   return kTRUE;
 }
index a27db54..c78166e 100644 (file)
@@ -59,6 +59,8 @@ public :
   Double_t YClu         (Int_t i)const {return fYClu[i];}                                          //Getter of y clu
   Bool_t   ClCk         (Int_t i)const {return fClCk[i];}                                          //Getter of cluster flags
   Double_t CkovFit      ()const {return fCkovFit;}                                                 //Getter of ckov angle fitted
+  Double_t ThTrkIn      ()const {return fThTrkIn;}                                                 //Getter of theta started of the track
+  Double_t PhTrkIn      ()const {return fPhTrkIn;}                                                 //Getter of phi started of the track
   Double_t ThTrkFit     ()const {return fThTrkFit;}                                                //Getter of theta fitted of the track
   Double_t PhTrkFit     ()const {return fPhTrkFit;}                                                //Getter of phi fitted of the track
 //
@@ -74,6 +76,8 @@ protected:
   Double_t *fXClu;                             //container for x clus position
   Double_t *fYClu;                             //container for y clus position
   Bool_t   *fClCk;                             //flag if cluster is used in fitting
+  Double_t fThTrkIn;                           //theta started from ShapeModel
+  Double_t fPhTrkIn;                           //phi   started from ShapeModel
   Double_t fThTrkFit;                          //theta fitted of the track
   Double_t fPhTrkFit;                          //phi   fitted of the track
   Double_t fCkovFit;                           //estimated ring Cherenkov angle
@@ -86,7 +90,7 @@ private:
   AliHMPIDReconHTA(const AliHMPIDReconHTA& r);              //dummy copy constructor
   AliHMPIDReconHTA &operator=(const AliHMPIDReconHTA& r);   //dummy assignment operator
 //
-  ClassDef(AliHMPIDReconHTA,1)
+  ClassDef(AliHMPIDReconHTA,2)
 };
 
 #endif // #ifdef AliHMPIDReconHTA_cxx
index 4777e4b..351716f 100644 (file)
@@ -1,7 +1,8 @@
 #include "AliHMPIDTracker.h"     //class header
-#include "AliHMPIDtrack.h"     //class header
+#include "AliHMPIDtrack.h"       //class header
 #include "AliHMPIDCluster.h"     //GetTrackPoint(),PropagateBack() 
 #include "AliHMPIDParam.h"       //GetTrackPoint(),PropagateBack()
+#include "AliHMPIDPid.h"         //Recon(),reconHTA()
 #include "AliHMPIDRecon.h"       //Recon()
 #include "AliHMPIDReconHTA.h"    //ReconHTA()
 #include <AliESDEvent.h>         //PropagateBack(),Recon()  
@@ -255,6 +256,13 @@ Int_t AliHMPIDTracker::Recon(AliESDEvent *pEsd,TObjArray *pClus,TObjArray *pNmea
     //
     recon.SetImpPC(xPc,yPc);                                                                     //store track impact to PC
     recon.CkovAngle(pTrk,(TClonesArray *)pClus->At(cham),index,nmean);                           //search for Cerenkov angle of this track
+    
+    AliHMPIDPid pID;
+    Double_t prob[5];
+    pID.FindPid(pTrk,5,prob);
+    pTrk->SetHMPIDpid(prob);
+//      Printf(" Prob e- %6.2f mu %6.2f pi %6.2f k %6.2f p %6.2f",prob[0]*100,prob[1]*100,prob[2]*100,prob[3]*100,prob[4]*100);
+
   }//iTrk
 
   return 0; // error code: 0=no error;
@@ -345,6 +353,12 @@ Int_t AliHMPIDTracker::ReconHiddenTrk(AliESDEvent *pEsd,TObjArray *pClus,TObjArr
     //
 //    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
+    
+    AliHMPIDPid pID;
+    Double_t prob[5];
+    pID.FindPid(pTrk,5,prob);
+    pTrk->SetHMPIDpid(prob);
+//      Printf(" Prob e- %6.2f mu %6.2f pi %6.2f k %6.2f p %6.2f",prob[0]*100,prob[1]*100,prob[2]*100,prob[3]*100,prob[4]*100);
   }//iTrk
 
   return 0; // error code: 0=no error;
index ecb6130..0374017 100644 (file)
@@ -31,8 +31,8 @@ public:
   static Int_t       IntTrkCha     (AliESDtrack *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       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(AliESDEvent *pEsd,TObjArray *pClus,TObjArray *pNmean, TObjArray *pQthre);//do actual job with Hidden Track Algorithm    
+  static Int_t       Recon         (AliESDEvent *pEsd,TObjArray *pClus,TObjArray *pNmean=0,TObjArray *pQthre=0);//do actual job, returns status code  
+  static Int_t       ReconHiddenTrk(AliESDEvent *pEsd,TObjArray *pClus,TObjArray *pNmean=0,TObjArray *pQthre=0);//do actual job with Hidden Track Algorithm    
   
   
 protected:
index e4f079f..3efbca4 100644 (file)
@@ -55,6 +55,7 @@ void SimEsd(AliLoader *pHL,AliESDEvent *pEsd,Bool_t htaCheck)
     Printf("-----------------------------------------------");
 }
   AliHMPIDTracker::SetFieldMap(gAL->GetAliRun()->Field(),kTRUE);
+  AliHMPIDTracker pTracker;
   AliHMPID *pH=(AliHMPID*)gAL->GetAliRun()->GetDetector("HMPID");
   Int_t iNevt=gAL->GetNumberOfEvents();
   pEsd->SetMagneticField(AliHMPIDTracker::GetBz());
@@ -69,7 +70,7 @@ void SimEsd(AliLoader *pHL,AliESDEvent *pEsd,Bool_t htaCheck)
       if(pTrack->GetPDG()->Charge()==0) continue;
       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 
+      Int_t iCh=pTracker.IntTrkCha(&trk,xPc,yPc,xRa,yRa,thRa,phRa);         //get chamber intersected by this track 
       if(iCh<0) {
         trk.SetHMPIDtrk(0,0,0,0);                                                                //no intersection found
         trk.SetHMPIDcluIdx   (99,99999);                                                         //chamber not found, mip not yet considered
@@ -78,52 +79,20 @@ void SimEsd(AliLoader *pHL,AliESDEvent *pEsd,Bool_t htaCheck)
       }
       nTrkHMPID++;
       trk.SetHMPIDcluIdx   (iCh,99999);                                                          //chamber not found, mip not yet considered
+      
+      if(phRa<0) phRa += TMath::TwoPi(); // to be verified
+      
       trk.SetHMPIDtrk(xRa,yRa,thRa,phRa);                                                        //store initial infos
       pEsd->AddTrack(&trk);
     }// track loop
     
     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);
+    Int_t status;
+    if(!htaCheck) status = pTracker.Recon         (pEsd,pH->CluLst(),pNmean,pQthre);
+    else          status = pTracker.ReconHiddenTrk(pEsd,pH->CluLst(),pNmean,pQthre);
 
-        // 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);
-      pTrk->SetHMPIDpid(prob);
-//      Printf(" Prob e- %6.2f mu %6.2f pi %6.2f k %6.2f p %6.2f",prob[0]*100,prob[1]*100,prob[2]*100,prob[3]*100,prob[4]*100);
-    }
+    if(status !=0) continue;
     
     gEsdTr->Fill();
     pEsd->Reset();
@@ -136,79 +105,7 @@ void SimEsd(AliLoader *pHL,AliESDEvent *pEsd,Bool_t htaCheck)
     delete fout;
   }
   gAL->UnloadHeader();  gAL->UnloadKinematics();
-}//Esd()
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void SimEsdHidden(AliLoader *pHL,AliESDEvent *pEsd)
-{
-  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("----------------------------------------------");
-  AliHMPIDTracker::SetFieldMap(gAL->GetAliRun()->Field(),kTRUE);
-  AliHMPID *pH=(AliHMPID*)gAL->GetAliRun()->GetDetector("HMPID");
-  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%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); 
-      
-      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
-      if(iCh<0) continue;                                                           //no intersection at all, go after next track
-      trk.SetHMPIDtrk(xRa,yRa,thRa,phRa);                                                        //store initial infos
-      Float_t radX,radY,thetaTrk,phiTrk;
-      trk.GetHMPIDtrk(radX,radY,thetaTrk,phiTrk);
-//      Printf("simulated track theta %f phi %f",thetaTrk*rd,phiTrk*rd);
-      TObjArray *pClus = pH->CluLst();
-      
-      if(AliHMPIDTracker::ReconHiddenTrk(iCh,0,&trk,(TClonesArray *)pClus->At(iCh),pNmean,pQthre)!=0) continue;
-      
-      pEsd->AddTrack(&trk);
-      Double_t thetaCerSim = TMath::ACos(pTrack->Energy()/(1.292*pTrack->P()));
-//      Printf(" theta Cerenkov simulated     %f",thetaCerSim);
-//      Printf(" theta Cerenkov reconstructed %f",trk.GetHMPIDsignal());
-      Float_t thetaFit,phiFit;
-      trk.GetHMPIDtrk(radX,radY,thetaFit,phiFit);
-//      Printf("reconstr. track theta %f phi %f",thetaFit*rd,phiFit*rd);
-      
-      // fill useful histos
-      hdC->Fill(trk.GetHMPIDsignal()-thetaCerSim);
-      hCer->Fill(trk.GetHMPIDsignal());
-      htvsp->Fill(pTrack->P(),trk.GetHMPIDsignal());
-      hdth->Fill((thetaFit-thetaTrk)*1000);
-      hdph->Fill((phiFit-phiTrk)*1000);
-    }// track loop
-    pEsd->SetMagneticField(AliHMPIDTracker::GetBz());
-    gEsdTr->Fill();
-    pEsd->Reset();
-  }// event loop
-  Printf("Events processed %i",iEvt);
-  fout->Write();
-  fout->Close();
-  delete fout;
-  gAL->UnloadHeader();  gAL->UnloadKinematics();
-}//EsdHidden()
+}//SimEsd()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 Bool_t OpenCalib()
 {