for(Int_t iCh=0;iCh<7;iCh++){
// evaluate High Voltage
-// TObjArray *pHV=(TObjArray*)pMap->GetValue(Form("HMP_DET/HMP_MP%i/HMP_MP%i_PW/HMP_MP%i_SEC0/HMP_MP%i_SEC0_HV.actual.vMon",iCh,iCh,iCh,iCh)); //HV
+ TObjArray *pHV=(TObjArray*)pMap->GetValue(Form("HMP_DET/HMP_MP%i/HMP_MP%i_PW/HMP_MP%i_SEC0/HMP_MP%i_SEC0_HV.actual.vMon",iCh,iCh,iCh,iCh)); TIter nextHV(pHV);
+ TGraph *pGrHV=new TGraph; cnt=0;
+ while((pVal=(AliDCSValue*)nextHV())) pGrHV->SetPoint(cnt++,pVal->GetTimeStamp(),pVal->GetFloat()); //P
+ if( cnt!=0) pGrHV->Fit(new TF1(Form("HV%i",iCh),"2000+x*[0]",fStartTime,fEndTime),"Q"); //clm: if no DCS map entry don't fit
+ delete pGrHV;
+
// evaluate Pressure
TObjArray *pP =(TObjArray*)pMap->GetValue(Form("HMP_DET/HMP_MP%i/HMP_MP%i_GAS/HMP_MP%i_GAS_PMWC.actual.value" ,iCh,iCh,iCh)); TIter nextP(pP);
TGraph *pGrP=new TGraph; cnt=0;
while((pVal=(AliDCSValue*)nextP())) pGrP->SetPoint(cnt++,pVal->GetTimeStamp(),pVal->GetFloat()); //P
if( cnt!=0) pGrP->Fit(new TF1(Form("P%i",iCh),"1005+x*[0]",fStartTime,fEndTime),"Q"); //clm: if no DCS map entry don't fit
- delete pGrP;
+ delete pGrP;
+
// evaluate Qthre
- arQthre.AddAt(new TF1(Form("HMP_Qthre%i",iCh),"100",fStartTime,fEndTime),iCh);
+ arQthre.AddAt(new TF1(Form("HMP_Qthre%i",iCh),Form("3*10^(3.01e-3*HV%i - 4.72)+170745848*exp(-P%i*0.0162012)",iCh),fStartTime,fEndTime),iCh);
+
+
// evaluate UserCut
Int_t nSigmaUserCut = 3;
TObject *pUserCut = new TObject();pUserCut->SetUniqueID(nSigmaUserCut);
arUserCut.AddAt(pUserCut,iCh);
+
// evaluate Temperatures
for(Int_t iRad=0;iRad<3;iRad++){
TObjArray *pT1=(TObjArray*)pMap->GetValue(Form("HMP_DET/HMP_MP%i/HMP_MP%i_LIQ_LOOP.actual.sensors.Rad%iIn_Temp",iCh,iCh,iRad)); TIter nextT1(pT1);//Tin
}
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t nmean)
+void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t nmean,Double_t qthre)
{
// Pattern recognition method based on Hough transform
// Arguments: pTrk - track for which Ckov angle is to be found
for (Int_t iClu=0; iClu<pCluLst->GetEntriesFast();iClu++){//clusters loop
AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu); //get pointer to current cluster
chId=pClu->Ch();
- if(pClu->Q()>pParam->QCut()){ //charge compartible with MIP clusters
+ if(pClu->Q()>qthre){ //charge compartible with MIP clusters
Float_t dX=fPc.X()-pClu->X(),dY=fPc.Y()-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 compatible with photon cluster
thetaCer= dirCkovTRS.Theta(); //actual value of thetaCerenkov of the photon
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Double_t AliHMPIDRecon::FindRingArea(Double_t ckovAng)const
+Double_t AliHMPIDRecon::FindRingArea(Double_t ckovAngMin,Double_t ckovAngMax)const
{
-// Find area inside the cerenkov ring which lays inside PCs
-// Arguments: ckovAng - cerenkov angle
+// Find area between 2 cerenkov angles in the PC acceptance
+// Arguments: ckovAngMin - cerenkov angle Min
+// Arguments: ckovAngMax - cerenkov angle Max
// Returns: area of the ring in cm^2 for given theta ckov
const Int_t kN=100;
- Double_t area=0;
+ Int_t np=0;
+ Double_t xP[2*kN],yP[2*kN];
- TVector2 pos1=TracePhot(ckovAng,0);//trace the first photon
- for(Int_t i=1;i<kN;i++){
- TVector2 pos2=TracePhot(ckovAng,Double_t(TMath::TwoPi()*(i+1)/kN));//trace the next photon
- if(pos1.X()==-999||pos2.X()==-999) return 0; //no area: open ring
- area+=TMath::Abs((pos1-fTrkPos).X()*(pos2-fTrkPos).Y()-(pos1-fTrkPos).Y()*(pos2-fTrkPos).X()); //add area of the triangle...
- pos1=pos2; // actual photon becomes the first one
+ Double_t area=0;
+//--- find points from first ring
+ for(Int_t i=0;i<kN;i++){
+ TVector2 pos=TracePhot(ckovAngMin,Double_t(TMath::TwoPi()*(i+1)/kN));//trace the next photon
+ if(pos.X()==-999) continue; //no area: open ring
+ if(AliHMPIDParam::IsInside(pos.X(),pos.Y(),0)) continue;
+ xP[np] = pos.X();
+ yP[np] = pos.Y();
+ np++;
+ }
+//--- find points from last ring
+ for(Int_t i=kN-1;i>=0;i--){
+ TVector2 pos=TracePhot(ckovAngMax,Double_t(TMath::TwoPi()*(i+1)/kN));//trace the next photon
+ if(pos.X()==-999) continue;
+ if(AliHMPIDParam::IsInside(pos.X(),pos.Y(),0)) continue;
+ xP[np] = pos.X();
+ yP[np] = pos.Y();
+ np++;
}
+//--calculate delta area from array of points...
+// for(Int_t i=0;i<np;i++) {
+ area = 1;
area*=0.5;
return area;
}//FindRingArea()
Int_t iInsideCnt = 0; //count photons which Theta ckov inside the window
for(Int_t i=0;i<fPhotCnt;i++){//photon candidates loop
+ fPhotFlag[i] = 0;
if(fPhotCkov[i] >= tmin && fPhotCkov[i] <= tmax) {
fPhotFlag[i]=2;
iInsideCnt++;
Int_t bin = (Int_t)(0.5+angle/(fDTheta));
Double_t weight=1.;
if(fIsWEIGHT){
- Double_t lowerlimit = ((Double_t)bin)*fDTheta - 0.5*fDTheta; Double_t upperlimit = ((Double_t)bin)*fDTheta + 0.5*fDTheta;
- Double_t diffArea = FindRingArea(upperlimit)-FindRingArea(lowerlimit);
+ Double_t lowerlimit = ((Double_t)bin)*fDTheta - 0.5*fDTheta; Double_t upperlimit = ((Double_t)bin)*fDTheta + 0.5*fDTheta;
+ Double_t diffArea = FindRingArea(lowerlimit,upperlimit);
if(diffArea>0) weight = 1./diffArea;
}
photsw->Fill(angle,weight);
// From here HTA....
//
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Bool_t AliHMPIDRecon::CkovHiddenTrk(AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t nmean)
+Bool_t AliHMPIDRecon::CkovHiddenTrk(AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t nmean, Double_t qthre)
{
// Pattern recognition method without any infos from tracking:HTA (Hidden Track Algorithm)...
// The method finds in the chmber the cluster with the highest charge
// nmean - mean freon ref. index
// Returns: - 0=ok,1=not fitted
- AliHMPIDParam *pParam=AliHMPIDParam::Instance();
-
fRadNmean=nmean;
if(pCluLst->GetEntriesFast()>100) return kFALSE; //boundary check for CluX,CluY...
}//clusters loop
fNClu = pCluLst->GetEntriesFast();
- if(qRef>pParam->QCut()){ //charge compartible with MIP clusters
+ if(qRef>qthre){ //charge compartible with MIP clusters
fIdxMip = mipId;
fClCk[mipId] = kFALSE;
fMipX = mipX; fMipY=mipY; fMipQ = qRef;
virtual ~AliHMPIDRecon() {}
- void CkovAngle (AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t nmean ); //reconstructed Theta Cerenkov
+ void CkovAngle (AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t nmean, Double_t qthre); //reconstructed Theta Cerenkov
Bool_t FindPhotCkov (Double_t cluX,Double_t cluY,Double_t &thetaCer,Double_t &phiCer ); //find ckov angle for single photon candidate
Double_t FindRingCkov (Int_t iNclus ); //best ckov for ring formed by found photon candidates
- Double_t FindRingArea (Double_t ckov )const;//estimated area of ring in cm^2
+ Double_t FindRingArea (Double_t ckovAngMin,Double_t ckovAngMax )const;//estimated area of delta ring in cm^2 to weight Hough Transform
Int_t FlagPhot (Double_t ckov ); //is photon ckov near most probable track ckov
Double_t HoughResponse( ); //most probable track ckov angle
void Propagate (const TVector3 dir, TVector3 &pos,Double_t z )const;//propagate photon alogn the line
Double_t Sigma2 (Double_t ckovTh,Double_t ckovPh )const;//photon candidate sigma^2
enum ETrackingFlags {kMipDistCut=-9,kMipQdcCut=-5,kNoPhotAccept=-11};
// HTA hidden track algorithm
- Bool_t CkovHiddenTrk (AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t nmean); //Pattern recognition without trackinf information
+ Bool_t CkovHiddenTrk (AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t nmean,Double_t qthre);//Pattern recognition without trackinf information
Bool_t CluPreFilter (TClonesArray *pClu ); //Pre clustering filter to cut bkg clusters
Bool_t DoRecHiddenTrk (TClonesArray *pClu ); //Calling to the fitted procedures
Bool_t FitEllipse (Double_t &phiRec ); //Fit clusters with a conical section (kTRUE only for ellipses)
if(!pNmeanEnt) AliFatal("No Nmean C6F14 ");
if(!pQthreEnt) AliFatal("No Qthre");
- return Recon(pEsd,fClu,(TObjArray*)pNmeanEnt->GetObject());
+ return Recon(pEsd,fClu,(TObjArray*)pNmeanEnt->GetObject(),(TObjArray*)pQthreEnt->GetObject());
}//PropagateBack()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Int_t AliHMPIDTracker::Recon(AliESDEvent *pEsd,TObjArray *pClus,TObjArray *pNmean)
+Int_t AliHMPIDTracker::Recon(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)
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
Double_t nmean=((TF1*)pNmean->At(3*cham))->Eval(pEsd->GetTimeStamp()); //C6F14 Nmean for this chamber
+ Double_t qthre=((TF1*)pQthre->At(cham)) ->Eval(pEsd->GetTimeStamp()); //Qthre for this chamber
recon.SetImpPC(xPc,yPc); //store track impact to PC
- recon.CkovAngle(pTrk,(TClonesArray *)pClus->At(cham),nmean); //search for Cerenkov angle of this track
+ recon.CkovAngle(pTrk,(TClonesArray *)pClus->At(cham),nmean,qthre); //search for Cerenkov angle of this track
+ Printf("AliHMPIDTracker::Recon: nmean %f, qthre %f",nmean,qthre);
} //ESD tracks loop
return 0; // error code: 0=no error;
}//Recon()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Int_t AliHMPIDTracker::ReconHiddenTrk(Int_t iCh,AliESDtrack *pTrk,TClonesArray *pCluLst,TObjArray *pNmean)
+Int_t AliHMPIDTracker::ReconHiddenTrk(Int_t iCh,AliESDtrack *pTrk,TClonesArray *pCluLst,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)
+// 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
AliHMPIDRecon recon; //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=((TF1*)pQthre->At(iCh)) ->Eval(0); //C6F14 Nmean for this chamber
if(pCluLst->GetEntriesFast()<4) return 1; //min 4 clusters (3 + 1 mip) to find a ring!
- if(recon.CkovHiddenTrk(pTrk,pCluLst,nmean)) return 0; //search for track parameters and Cerenkov angle of this track
+ if(recon.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;
}//Recon()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
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 Recon (AliESDEvent *pEsd,TObjArray *pCluAll,TObjArray *pNmean=0); //do actual job, returns status code
- static Int_t ReconHiddenTrk(Int_t iCh,AliESDtrack *pTrk,TClonesArray *pClus,TObjArray *pNmean);//do actual job with Hidden Track Algorithm
+ 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,AliESDtrack *pTrk,TClonesArray *pClus,TObjArray *pNmean, TObjArray *pQthre);//do actual job with Hidden Track Algorithm
protected:
TObjArray *fClu; //! each chamber holds it's one list of clusters
ClassDef(AliHMPIDTracker,0)