]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - HMPID/HESDfromKin.C
PID probability calculated using refractive index evaluated from OCDB information...
[u/mrichter/AliRoot.git] / HMPID / HESDfromKin.C
index ec1377a5419a5824549d47d1cbd260649aaf229a..042f4ca2c77ce4fb6e4bd56724c21882395a711a 100644 (file)
@@ -28,144 +28,102 @@ 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("-----------------------------------------------");
-  AliHMPIDTracker::SetFieldMap(gAL->GetAliRun()->Field(),kTRUE);
+  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("-----------------------------------------------");
+}
+
+  AliGRPManager *grpMan = new AliGRPManager();
+  
+  grpMan->ReadGRPEntry();
+  grpMan->SetMagField();
+
+  AliHMPIDTracker pTracker;
   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;
+//      Printf("track n. %i",i);
       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
         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
+      
+      if(phRa<0) phRa += TMath::TwoPi(); // to be verified
+      
+      trk.SetHMPIDtrk(xPc,yPc,thRa,phRa);                                                        //store initial infos
+      trk.SetLabel(i);
       pEsd->AddTrack(&trk);
-    }// track loop
     
-    AliHMPIDTracker::Recon(pEsd,pH->CluLst(),pNmean,pQthre);
-// Perform PID
+      Int_t status;
+      if(!htaCheck) status = pTracker.ReconFromKin  (pEsd,pH->CluLst(),pNmean,pQthre);
+      else          status = pTracker.ReconHiddenTrk(pEsd,pH->CluLst(),pNmean,pQthre);
+
+//      Printf("status %i",status);
+      if(status !=0) continue;
+
     
-    for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){                                       //ESD tracks loop
-      AliESDtrack *pTrk = pEsd->GetTrack(iTrk);                                                    //get reconstructed track    
-      AliHMPIDPid pID;
-      Double_t prob[5];
-      pID.FindPid(pTrk,5,prob);
-      trk.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);
-    }
+    }// track loop
     
-    gEsdTr->Fill();
-    pEsd->Reset();
+    if(!(iEvt%50)) Printf("Number of events processed: %i with tracks %i in HMPID",iEvt,nTrkHMPID);
+//      Printf("Number of events processed: %i with tracks %i in HMPID",iEvt,nTrkHMPID);
     
-  }// event loop
-  Printf("Events processed %i",iEvt);
-  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 mtid=-1;
-  Int_t iNevt=gAL->GetNumberOfEvents();
-  
-  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();
-    
-    for(Int_t i=0;i<pStack->GetNtrack();i++){
-      
-      TParticle *pTrack=pStack->Particle(i); 
-      mtid=pTrack->GetFirstMother();
-      
-      if(mtid>=0) continue; // only primaries
-      
-      //find the chamber that intersects HMPID
-      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;
+  if(htaCheck) {
+    fout->Write();
+    fout->Close();
+    delete fout;
+  }
   gAL->UnloadHeader();  gAL->UnloadKinematics();
-}//EsdHidden()
+}//SimEsd()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 Bool_t OpenCalib()
 {
   AliCDBManager* pCDB = AliCDBManager::Instance();
-  pCDB->SetDefaultStorage("local://$HOME");
+  pCDB->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+  pCDB->SetSpecificStorage("GRP/GRP/Data", Form("local://%s",gSystem->pwd()));
+  pCDB->SetRun(0);
   AliCDBEntry *pQthreEnt=pCDB->Get("HMPID/Calib/Qthre",0);
   AliCDBEntry *pNmeanEnt=pCDB->Get("HMPID/Calib/Nmean",0);
   
@@ -177,4 +135,3 @@ Bool_t OpenCalib()
   if(!pQthre || !pNmean) return kFALSE;  
   return kTRUE;
 }
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++