Optimization of the reconstruction code.
authordibari <dibari@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 31 Mar 2008 09:33:09 +0000 (09:33 +0000)
committerdibari <dibari@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 31 Mar 2008 09:33:09 +0000 (09:33 +0000)
HMPID/AliHMPIDRecon.cxx
HMPID/AliHMPIDRecon.h
HMPID/AliHMPIDTracker.cxx
HMPID/AliHMPIDTracker.h
HMPID/HESDfromKin.C
HMPID/Hconfig.C
HMPID/Hdisp.C

index a6c5892..5949ee6 100644 (file)
@@ -107,6 +107,10 @@ void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t n
 //              pCluLst  - list of clusters for this chamber   
 //   Returns:            - track ckov angle, [rad], 
     
+
+  const Int_t nMinPhotAcc = 3;                      // Minimum number of photons required to perform the pattern recognition
+  
+  
   Int_t nClusTot = pCluLst->GetEntries();
   if(nClusTot>fParam->MultCut()) fIsWEIGHT = kTRUE; // offset to take into account bkg in reconstruction
   else                           fIsWEIGHT = kFALSE;
@@ -137,14 +141,25 @@ void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t n
       }
     }
   }//clusters loop
+  
   fMipPos.Set(mipX,mipY);
-  if(fPhotCnt<=3) pTrk->SetHMPIDsignal(kNoPhotAccept);                                        //no reconstruction with <=3 photon candidates
+  
+  if(fPhotCnt<=nMinPhotAcc) {                                                                 //no reconstruction with <=3 photon candidates
+    pTrk->SetHMPIDsignal(kNoPhotAccept);                                                      //set the appropriate flag
+    pTrk->SetHMPIDmip(mipX,mipY,mipQ,fPhotCnt);                                               //store mip info 
+    return;
+  }
+  
+  if(mipId==-1)              {pTrk->SetHMPIDsignal(kMipQdcCut);  return;}                     //no clusters with QDC more the threshold at all
+  if(dMin>fParam->DistCut()) {pTrk->SetHMPIDsignal(kMipDistCut); return;}                     //closest cluster with enough charge is still too far from intersection
+  
+//PATTERN RECOGNITION STARTED: 
+  
   Int_t iNrec=FlagPhot(HoughResponse());                                                      //flag photons according to individual theta ckov with respect to most probable
   pTrk->SetHMPIDmip(mipX,mipY,mipQ,iNrec);                                                    //store mip info 
 
-  if(mipId==-1)              {pTrk->SetHMPIDsignal(kMipQdcCut);  return;}                     //no clusters with QDC more the threshold at all
-  if(dMin>fParam->DistCut()) {pTrk->SetHMPIDsignal(kMipDistCut); return;}                     //closest cluster with enough charge is still too far from intersection
   pTrk->SetHMPIDcluIdx(chId,mipId);                                                           //set index of cluster
+  
   if(iNrec<1){
     pTrk->SetHMPIDsignal(kNoPhotAccept);                                                      //no photon candidates are accepted
   }
index 6316e4d..ec1bc2a 100644 (file)
@@ -53,7 +53,7 @@ public :
                                 {fPc.Set(xPc,yPc);}                                                //set track impact to PC 
   void     SetMip       (Double_t xmip,Double_t ymip                                        )
                                 {fMipPos.Set(xmip,ymip);}                                          //set track impact to PC
-  enum ETrackingFlags {kMipDistCut=-9,kMipQdcCut=-5,kNoPhotAccept=-11};
+  enum ETrackingFlags {kNotPerformed=-20,kMipDistCut=-9,kMipQdcCut=-5,kNoPhotAccept=-11};
 //
 protected:
   Int_t     fPhotCnt;                           // counter of photons candidate
index 7b7b4fb..0907bfd 100644 (file)
@@ -46,13 +46,12 @@ Bool_t AliHMPIDTracker::GetTrackPoint(Int_t idx, AliTrackPoint& point) const
   return kTRUE;
 }//GetTrackPoint()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Int_t AliHMPIDTracker::IntTrkCha(AliESDtrack *pTrk,Float_t &xPc,Float_t &yPc)
+Int_t AliHMPIDTracker::IntTrkCha(AliESDtrack *pTrk,Float_t &xPc,Float_t &yPc,Float_t &xRa,Float_t &yRa,Float_t &theta,Float_t &phi)
 {
 // Static method to find intersection in between given track and HMPID chambers
 // Arguments: pTrk- ESD track; xPc,yPc- track intersection with PC in LORS [cm]
 //   Returns: intersected chamber ID or -1
   AliHMPIDParam *pParam=AliHMPIDParam::Instance();
-  Float_t xRa=0,yRa=0,theta=0,phi=0;                                                            //track intersection at PC and angles at RAD, LORS  
   for(Int_t i=AliHMPIDParam::kMinCh;i<=AliHMPIDParam::kMaxCh;i++){                              //chambers loop
     Double_t p1[3],n1[3]; pParam->Norm(i,n1); pParam->Point(i,p1,AliHMPIDParam::kRad);          //point & norm  for middle of radiator plane
     Double_t p2[3],n2[3]; pParam->Norm(i,n2); pParam->Point(i,p2,AliHMPIDParam::kPc);           //point & norm  for entrance to PC plane
@@ -61,12 +60,9 @@ Int_t AliHMPIDTracker::IntTrkCha(AliESDtrack *pTrk,Float_t &xPc,Float_t &yPc)
     pParam->Mars2LorsVec(i,n1,theta,phi);                                                       //track angles at RAD
     pParam->Mars2Lors   (i,p1,xRa,yRa);                                                         //TRKxRAD position
     pParam->Mars2Lors   (i,p2,xPc,yPc);                                                         //TRKxPC position
-    if(AliHMPIDParam::IsInside(xPc,yPc,pParam->DistCut())==kFALSE) continue;                    //not in active area  
-    pTrk->SetHMPIDtrk      (xRa,yRa,theta,phi);                                                 //store track intersection info
-    pTrk->SetHMPIDcluIdx   (i,0);
-    return i;
+    if(AliHMPIDParam::IsInside(xPc,yPc,pParam->DistCut())==kTRUE) return i;                     //return intersected chamber  
   }                                                                                             //chambers loop
-  return -1; //no intersection with HMPID chambers
+  return -1;                                                                                    //no intersection with HMPID chambers
 }//IntTrkCha()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 Int_t AliHMPIDTracker::LoadClusters(TTree *pCluTree)
@@ -98,11 +94,17 @@ Int_t AliHMPIDTracker::Recon(AliESDEvent *pEsd,TObjArray *pClus,TObjArray *pNmea
 // Arguments: pEsd- pointer ESD; pClu- pointer to clusters for all chambers; pNmean - pointer to all function Nmean=f(time)
 //   Returns: error code, 0 if no errors   
   AliHMPIDRecon recon;                                                                       //instance of reconstruction class, nothing important in ctor
-  Float_t xPc,yPc;
+  Float_t xPc,yPc,xRa,yRa,theta,phi;
   for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){                                       //ESD tracks loop
     AliESDtrack *pTrk = pEsd->GetTrack(iTrk);                                                    //get reconstructed track    
-    Int_t cham=IntTrkCha(pTrk,xPc,yPc);                                                          //get chamber intersected by this track 
-    if(cham<0) continue;                                                                         //no intersection at all, go after next track
+    Int_t cham=IntTrkCha(pTrk,xPc,yPc,xRa,yRa,theta,phi);                                        //get chamber intersected by this track 
+    if(cham<0) {                                                                                 //no intersection at all, go after next track
+      pTrk->SetHMPIDtrk(0,0,0,0);                                                                //no intersection found
+      pTrk->SetHMPIDcluIdx   (99,99999);                                                         //chamber not found, mip not yet considered
+      pTrk->SetHMPIDsignal(AliHMPIDRecon::kNotPerformed);                                        //ring reconstruction not yet performed
+      continue;                                                                         
+    }
+    pTrk->SetHMPIDtrk(xRa,yRa,theta,phi);                                                        //store initial infos
     Double_t nmean=((TF1*)pNmean->At(3*cham))->Eval(pEsd->GetTimeStamp());                       //C6F14 Nmean for this chamber
     Int_t hvsec = AliHMPIDParam::InHVSector(xPc,yPc);
     Double_t qthre=((TF1*)pQthre->At(6*cham+hvsec))->Eval(pEsd->GetTimeStamp());
index 705b36d..fd97e97 100644 (file)
@@ -25,7 +25,7 @@ public:
          Int_t       RefitInward    (AliESDEvent *                   )       {return 0;} //pure virtual from AliTracker 
          void        UnloadClusters (                           )       {         } //pure virtual from AliTracker 
 //private part  
-  static Int_t       IntTrkCha     (AliESDtrack *pTrk,Float_t &xPc,Float_t &yPc        );              //find track-PC intersection, retuns chamber ID
+  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       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    
   
index 6bb3999..ea63651 100644 (file)
@@ -56,10 +56,17 @@ void SimEsd(AliLoader *pHL,AliESDEvent *pEsd)
       if(mtid>=0) continue; // only primaries
       if(pTrack->GetPDG()->Charge()==0) continue;
       AliESDtrack trk(pTrack); 
-      Float_t xPc,yPc;
-      Int_t iCh=AliHMPIDTracker::IntTrkCha(&trk,xPc,yPc);                           //get chamber intersected by this track 
-      if(iCh<0) continue;                                                           //no intersection at all, go after next track
+      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) {
+        trk.SetHMPIDtrk(0,0,0,0);                                                                //no intersection found
+        trk.SetHMPIDcluIdx   (99,99999);                                                         //chamber not found, mip not yet considered
+        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);
+      trk.SetHMPIDcluIdx   (iCh,99999);                                                          //chamber not found, mip not yet considered
+      trk.SetHMPIDtrk(xRa,yRa,thRa,phRa);                                                        //store initial infos
       pEsd->AddTrack(&trk);
       AliHMPIDTracker::Recon(pEsd,pH->CluLst(),pNmean,pQthre);
     }// track loop
@@ -105,15 +112,16 @@ void SimEsdHidden(AliLoader *pHL,AliESDEvent *pEsd)
       
       //find the chamber that intersects HMPID
       AliESDtrack trk(pTrack);
-      Float_t xPc,yPc;
-      Int_t iCh=AliHMPIDTracker::IntTrkCha(&trk,xPc,yPc);                           //get chamber intersected by this track 
+      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,&trk,(TClonesArray *)pClus->At(iCh),pNmean,pQthre)!=0) continue;
+      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()));
index 95ce276..e5c7442 100644 (file)
@@ -639,10 +639,10 @@ void HmpConfig::WriteBatch()
   if(fDetBG->GetButton(kTRD  )->GetState())  det+="TRD ";
   if(fDetBG->GetButton(kTOF  )->GetState())  det+="TOF ";
   if(!fVerBG->GetButton(kNo)->GetState())    det+="HMPID ";
-  char *sBatchName="Hbatch";
-  FILE *fp=fopen(Form("%s.C",sBatchName),"w"); if(!fp){Info("CreateBatch","Cannot open output file: %s.C",sBatchName);return;}
+  char *sBatchName="sim";
+  FILE *fp=fopen(Form("%s.C",sBatchName),"w"); if(!fp){Info("CreateSim","Cannot open output file: %s.C",sBatchName);return;}
   
-                                                    fprintf(fp,"void %s(Int_t iNevt,Bool_t isDbg,char *sCfg)\n{\n",sBatchName);
+                                                    fprintf(fp,"void %s(Int_t iNevt=1,Bool_t isDbg=kFALSE,char *sCfg=\"Config.C\")\n{\n",sBatchName);
                                                     fprintf(fp,"  gSystem->Exec(\"rm -rf hlt hough gphysi* fort* ZZZ* raw*\");  //remove garbage\n"); 
                                                     fprintf(fp,"  gBenchmark->Start(\"ALICE\"); TDatime time;      //start benchmarking\n\n");
                                                        
@@ -668,11 +668,25 @@ void HmpConfig::WriteBatch()
     else if(fRawBG->GetButton(kRoo)->GetState())    fprintf(fp,"  pSim->SetWriteRawData(\"%s\",\"raw.root\");     //raw data as ROOT\n",det.Data());
 
                                                     fprintf(fp,"  pSim->SetRunHLT(\"\");                           //no HLT stuff\n");   
-                                                    fprintf(fp,"  pSim->Run(iNevt);                              //run iNevt events\n  delete pSim;\n\n");
+                                                    fprintf(fp,"  pSim->SetQA(kFALSE);                             //no QA\n");
+                                                    fprintf(fp,"  pSim->Run(iNevt);                                //run iNevt events\n  delete pSim;\n\n");
   }//sim section
+                                                    fprintf(fp,"  cout<<\"!!!!!!!!!!!!Info in <sim.C>: Start time: \";time.Print();\n");
+                                                    fprintf(fp,"  cout<<\"!!!!!!!!!!!!Info in <sim.C>: Stop  time: \";time.Set();  time.Print();\n");
+                                                    fprintf(fp,"  gBenchmark->Show(\"ALICE\");\n");
+  
+                                                    fprintf(fp,"  gSystem->Exec(\"aliroot rec.C\");\n");
+                                                    fprintf(fp,"  gSystem->Exec(\"touch ZZZ______finished_______SSS\");\n}\n");
+  fclose(fp);  
+  char *sBatchName="rec";
+  FILE *fp=fopen(Form("%s.C",sBatchName),"w"); if(!fp){Info("CreateRec","Cannot open output file: %s.C",sBatchName);return;}
+  
+                                                    fprintf(fp,"void %s()\n{\n",sBatchName);
+                                                    fprintf(fp,"  gSystem->Exec(\"rm -rf RRR* \");  //remove garbage\n"); 
 
   if(fRecB->GetState()){
                                                     fprintf(fp,"  AliReconstruction *pRec=new AliReconstruction;\n");
+                                                    fprintf(fp,"  gBenchmark->Start(\"ALICE\"); TDatime time;      //start benchmarking\n\n");
                                                     
     //---------------------------------------------
     if     (fTrkBG->GetButton(kRecoPar)->GetState())
@@ -710,11 +724,11 @@ void HmpConfig::WriteBatch()
                                                     fprintf(fp,"  pRec->Run();delete pRec;\n\n");         
   }//rec part                                                       
 //benchmarks  
-                                                    fprintf(fp,"  cout<<\"!!!!!!!!!!!!Info in <my/Batch.C>: Start time: \";time.Print();\n");
-                                                    fprintf(fp,"  cout<<\"!!!!!!!!!!!!Info in <my/Batch.C>: Stop  time: \";time.Set();  time.Print();\n");
+                                                    fprintf(fp,"  cout<<\"!!!!!!!!!!!!Info in <rec.C>: Start time: \";time.Print();\n");
+                                                    fprintf(fp,"  cout<<\"!!!!!!!!!!!!Info in <rec.C>: Stop  time: \";time.Set();  time.Print();\n");
                                                     fprintf(fp,"  gBenchmark->Show(\"ALICE\");\n");
   
-                                                    fprintf(fp,"  gSystem->Exec(\"touch ZZZ______finished_______ZZZ\");\n}\n");
+                                                    fprintf(fp,"  gSystem->Exec(\"touch ZZZ______finished_______RRR\");\n}\n");
   fclose(fp);  
 }//WriteBatch()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
index d6bd9ca..2504d44 100644 (file)
@@ -320,10 +320,11 @@ void RenderEsd(AliESDEvent *pEsd)
   AliHMPIDRecon rec;
   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(); //get track and chamber intersected by it
-    if(ch<0) continue;                                                          //this track does not intersect any chamber
-    Float_t thRa,phRa,xRa,yRa; pTrk->GetHMPIDtrk(xRa,yRa,thRa,phRa);            //get info on current track
-    ch/=1000000;                            
-    Float_t xPc=0,yPc=0; AliHMPIDTracker::IntTrkCha(pTrk,xPc,yPc);              //find again intersection of track with PC--> it is not stored in ESD!
+    ch/=1000000;
+    if(ch<AliHMPIDParam::AliHMPIDParam::kMinCh||ch>AliHMPIDParam::kMaxCh) continue;//this track does not intersect any chamber
+    Float_t xPc,yPc,xRa,yRa,thRa,phRa; 
+    Int_t chamb = AliHMPIDTracker::IntTrkCha(pTrk,xPc,yPc,xRa,yRa,thRa,phRa);   //find again intersection of track with PC--> it is not stored in ESD!
+    if(ch!=chamb){(" CHAMBER MISMATCH: in ESDTrack chamber %i - in IntTrkCha chamber %i",ch,chamb);End();}
     Int_t npTrk = fRenTxC[ch]->SetNextPoint(xPc,yPc);                           //add this intersection point
     Float_t ckov=pTrk->GetHMPIDsignal();                                        //get ckov angle stored for this track  
     if(ckov>0){
@@ -568,9 +569,9 @@ void DisplayInfo(Int_t evt,Int_t px, Int_t py)
   else if (symbol==kTrack || symbol==kRing) {
     if(symbol==kRing) index = b->GetUniqueID();
     AliESDtrack *pTrk=fEsd->GetTrack(index);
-    Float_t thRa,phRa,xRa,yRa; pTrk->GetHMPIDtrk(xRa,yRa,thRa,phRa);
-    Float_t xPc,yPc; AliHMPIDTracker::IntTrkCha(pTrk,xPc,yPc);
-    text0="";text0.Append(Form("TRACK n.%d: x %6.2f y %6.2f at PC plane",index,xPc,yPc));
+    Float_t xPc,yPc,thRa,phRa,xRa,yRa;
+    Int_t ch = AliHMPIDTracker::IntTrkCha(pTrk,xPc,yPc,xRa,yRa,thRa,phRa);
+    text0="";text0.Append(Form("TRACK n.%d: x %6.2f y %6.2f at PC plane in chamber %i",index,xPc,yPc,ch));
     text2="";text2.Append(Form("p = %7.2f GeV/c",pTrk->GetP()));
     Float_t ckov=pTrk->GetHMPIDsignal();                             
     Double_t prob[5];