1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 //////////////////////////////////////////////////////////////////////////
18 // AliHMPIDReconHTA //
20 // HMPID class to perfom pattern recognition based on Hough transfrom //
21 // for single chamber //
22 //////////////////////////////////////////////////////////////////////////
24 #include "AliHMPIDReconHTA.h"//class header
25 #include "AliHMPIDCluster.h" //CkovHiddenTrk()
26 #include "AliHMPIDRecon.h" //FunMinPhot()
27 #include <TFile.h> //Database()
28 #include <TMinuit.h> //FitFree()
29 #include <TClonesArray.h> //CkovHiddenTrk()
30 #include <AliESDtrack.h> //CkovHiddenTrk()
31 #include <TH2F.h> //InitDatabase()
32 #include <TGraph.h> //ShapeModel()
33 #include <TSpline.h> //ShapeModel()
34 #include "TStopwatch.h" //
36 TH2F* AliHMPIDReconHTA::fgDatabase = new TH2F("deconv","database;d1;d2;thC+1000*thTrk",500,0,50,150,0,15);
39 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
40 AliHMPIDReconHTA::AliHMPIDReconHTA():
41 TTask("RichRec","RichPat"),
56 fParam(AliHMPIDParam::Instance())
61 fParam->SetRefIdx(fParam->MeanIdxRad()); // initialization of ref index to a default one
62 if(fgDatabase->GetEntries()<1) InitDatabase();
64 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
65 AliHMPIDReconHTA::~AliHMPIDReconHTA()
69 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
70 void AliHMPIDReconHTA::InitVars(Int_t n)
75 fXClu = new Double_t[n];
76 fYClu = new Double_t[n];
77 fClCk = new Bool_t[n];
78 for(Int_t i=0;i<n;i++) fClCk[i] = kTRUE;
81 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
82 void AliHMPIDReconHTA::DeleteVars()const
87 if(fXClu) delete fXClu;
88 if(fYClu) delete fYClu;
89 if(fClCk) delete fClCk;
91 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
92 Bool_t AliHMPIDReconHTA::CkovHiddenTrk(AliESDtrack *pTrk,TClonesArray *pCluLst,Int_t index, Double_t nmean)
94 // Pattern recognition method without any infos from tracking:HTA (Hidden Track Algorithm)...
95 // The method finds in the chmuber the cluster with the highest charge
96 // compatibile with a MIP, then the strategy is applied
97 // Arguments: pTrk - pointer to ESD track
98 // pCluLs - list of clusters for a given chamber
99 // pNmean - pointer to ref. index
100 // pQthre - pointer to qthre
101 // Returns: - 0=ok,1=not fitted
103 AliHMPIDParam *pParam = AliHMPIDParam::Instance();
105 if(!CluPreFilter(pCluLst)) return kFALSE;
110 fNClu = pCluLst->GetEntriesFast();
112 for (Int_t iClu=0;iClu<fNClu;iClu++){ //clusters loop
113 AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu); //get pointer to current cluster
114 fXClu[iClu] = pClu->X();fYClu[iClu] = pClu->Y(); //store x,y for fitting procedure
115 fClCk[iClu] = kTRUE; //all cluster are accepted at this stage to be reconstructed
122 sizeClu = pClu->Size();
124 fClCk[index] = kFALSE;
126 // Printf(" n. %d x %f y %f Q %f",iClu,pClu->X(),pClu->Y(),pClu->Q());
130 pParam->SetRefIdx(nmean);
133 Float_t xra,yra,th,ph; pTrk->GetHMPIDtrk(xra,yra,th,ph);
134 Printf(" simulated phi %6.2f ",ph*TMath::RadToDeg());
137 if(!DoRecHiddenTrk()) {
138 pTrk->SetHMPIDsignal(pParam->kNoPhotAccept);
140 } //Do track and ring reconstruction,if problems returns 1
141 Printf(" fitted phi %6.2f ",fPhTrkFit*TMath::RadToDeg());
143 pTrk->SetHMPIDtrk(fRadX,fRadY,fThTrkFit,fPhTrkFit); //store track intersection info
144 pTrk->SetHMPIDmip(fMipX,fMipY,(Int_t)fMipQ,fNClu); //store mip info
145 pTrk->SetHMPIDcluIdx(nCh,fIdxMip+1000*sizeClu); //set cham number, index of cluster + cluster size
146 pTrk->SetHMPIDsignal(fCkovFit); //find best Theta ckov for ring i.e. track
147 pTrk->SetHMPIDchi2(fCkovSig2); //errors squared
148 // Printf(" n clusters tot %i accepted %i",pCluLst->GetEntriesFast(),fNClu);
149 // Printf("CkovHiddenTrk: thetaC %f th %f ph %f",fCkovFit,fThTrkFit,fPhTrkFit);
154 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
155 Bool_t AliHMPIDReconHTA::DoRecHiddenTrk()
157 // Pattern recognition method without any infos from tracking...
158 // First a preclustering filter to avoid part of the noise
159 // Then only ellipsed-rings are fitted (no possibility,
160 // for the moment, to reconstruct very inclined tracks)
161 // Finally a fitting with (th,ph) free, starting by very close values
162 // previously evaluated.
165 Double_t thTrkRec,phiTrkRec,thetaCRec;
167 if(!FindShape(thTrkRec,phiTrkRec,thetaCRec)) return kFALSE;
169 if(!FitFree(thTrkRec,phiTrkRec)) {return kFALSE;}
173 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
174 Bool_t AliHMPIDReconHTA::CluPreFilter(TClonesArray *pCluLst)
176 // Pre-filter of bkg clusters
177 // Arguments: pSluLst - List of the clusters for a given chamber
178 // Returns: status - TRUE if filtering leaves enough photons, FALSE if not
180 Int_t nClusTot = pCluLst->GetEntriesFast();
181 if(nClusTot<4||nClusTot>100) {
188 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
189 Bool_t AliHMPIDReconHTA::FindShape(Double_t &thTrkRec,Double_t &phiTrkRec,Double_t &thetaCRec)
191 // Finds the estimates for phi and theta of the track and the ThetaCerenkov
192 // by using a database of the shapes of the rings
194 // Returns: thTrkRec - estimate of theta track
195 // phiTrkRec - estimate of phi track
196 // thetaCRec - estimate of ThetaCerenkov
197 // status - TRUE if a good solution is found, FALSE if not
199 Double_t *phiphot = new Double_t[fNClu];
200 Double_t *dist = new Double_t[fNClu];
201 Int_t *indphi = new Int_t[fNClu];
205 // Sort in phi angle...
206 // Printf(" mipX %f mipy %f",fMipX,fMipY);
207 for(Int_t i=0;i<fNClu;i++) {
213 phiphot[i] = (TMath::ATan2(fMipY-fYClu[i],fMipX-fXClu[i])+TMath::Pi())*TMath::RadToDeg();
214 dist[i]=TMath::Sqrt((fMipX-fXClu[i])*(fMipX-fXClu[i])+(fMipY-fYClu[i])*(fMipY-fYClu[i]));
215 // Printf(" n.%3i phiphot %f dist %f check %i",i,phiphot[i],dist[i],fClCk[i]);
218 TMath::Sort(fNClu,phiphot,indphi,kFALSE);
220 // Purify with a truncated mean;
224 for(Int_t i=0;i<fNClu;i++) {
225 if(!fClCk[indphi[i]]) continue; // Check if a good photon candidate or not
226 dMean +=dist[indphi[i]];
227 dMean2+=dist[indphi[i]]*dist[indphi[i]];
231 dMean /=(Double_t)np;
232 dMean2 /=(Double_t)np;
233 Double_t rms = TMath::Sqrt(dMean2 - dMean*dMean);
235 for(Int_t i=0;i<fNClu;i++) {
236 if(!fClCk[indphi[i]]) continue; // Check if a good photon candidate or not
237 if(TMath::Abs(dMean-dist[indphi[i]]) > 2*rms) {
238 fClCk[indphi[i]] = kFALSE;
244 // purify vectors for good photon candidates
247 Double_t *phiphotP = new Double_t[fNClu+1];
248 Double_t *distP = new Double_t[fNClu+1];
249 for(Int_t i=0;i<fNClu;i++) {
250 if(!fClCk[indphi[i]]) continue; // Check if a good photon candidate or not
251 phiphotP[npeff] = phiphot[indphi[i]];
252 distP[npeff] = dist[indphi[i]];
253 // Printf("n. %2i phi %f dist %f",npeff,phiphotP[npeff],distP[npeff]);
267 // for(Int_t i=0;i<npeff;i++) {Printf(" n. %d phiphot %f dist %f",i,phiphotP[i],distP[i]);}
271 if(ShapeModel(npeff,phiphotP,distP,xA,xB,phiTrkRec)) {
273 // Printf("FindShape: phi start %f xA %f yA %f",phiTrkRec*TMath::RadToDeg(),xA,xB);
274 if(xA < 50 && xB < 15) { // limits of the Database. See TH2F in InitDatabase()
276 Int_t bin = fgDatabase->FindBin(xA,xB);
278 Int_t compact = (Int_t)fgDatabase->GetBinContent(bin);
279 thetaCRec = (Double_t)(compact%1000);
280 thTrkRec = (Double_t)(compact/1000);
282 thTrkRec *= TMath::DegToRad();
283 thetaCRec *= TMath::DegToRad();
285 // Printf("FindShape: xA %f xB %f compact %d thTrk %f thC %f",xA,xB,compact,thTrkRec*TMath::RadToDeg(),thetaCRec*TMath::RadToDeg());
297 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
298 Bool_t AliHMPIDReconHTA::ShapeModel(Int_t np,Double_t *phiphot,Double_t *dist,Double_t &xA,Double_t &xB,Double_t &phiStart)
300 // Find a Spline curve to define dist. vs. phi angle
301 // in order to estimate the phi of the track
302 // Arguments: np - # points corresponding to # photon candidates
303 // dist - distance of each photon from MIP
304 // phiphot - phi of the photon in the DRS
305 // Returns: xA - min. distance from MIP
306 // xB - dist. from mip perpedicular to the major axis
307 // phiStart- estimate of the track phi
309 TGraph *phigr = new TGraph(np,phiphot,dist);
310 TSpline3 *sphi = new TSpline3("sphi",phigr);
311 if(!sphi) {Printf("Spline not created!Bye.");return kFALSE;}
313 Int_t locMin = TMath::LocMin(np,dist);
314 Int_t locMax = TMath::LocMax(np,dist);
316 Double_t minX = phiphot[locMin];
317 // Double_t minY = dist[locMin];
318 Double_t maxX = phiphot[locMax];
319 // Double_t maxY = dist[locMax];
321 Int_t ip[3] = {-1,0,1};
322 if(locMin==0 ) {ip[0]= 0;ip[1]= 1;ip[2]= 2;}
323 if(locMin==np-1) {ip[0]=-2;ip[1]=-1;ip[2]= 0;}
325 Double_t minXf = VertParab(phiphot[locMin+ip[0]],dist[locMin+ip[0]],
326 phiphot[locMin+ip[1]],dist[locMin+ip[1]],
327 phiphot[locMin+ip[2]],dist[locMin+ip[2]]);
328 if(minXf< phiphot[locMin+ip[0]] || minXf > phiphot[locMin+ip[2]]) minXf = minX;
330 ip[0]=-1;ip[1]=0;ip[2]=1;
331 if(locMax==0 ) {ip[0]= 0;ip[1]= 1;ip[2]= 2;}
332 if(locMax==np-1) {ip[0]=-2;ip[1]=-1;ip[2]= 0;}
334 Double_t maxXf = VertParab(phiphot[locMax+ip[0]],dist[locMax+ip[0]],
335 phiphot[locMax+ip[1]],dist[locMax+ip[1]],
336 phiphot[locMax+ip[2]],dist[locMax+ip[2]]);
337 if(maxXf< phiphot[locMax+ip[0]] || maxXf > phiphot[locMax+ip[2]]) maxXf = maxX;
340 if(TMath::Abs(maxXf-minXf)>30) {
341 xA = sphi->Eval(minXf);
342 if(xA < 0) xA = dist[sphi->FindX(xA)];
343 xB = sphi->Eval(minXf-90);
344 if(xB < 0) xB = dist[sphi->FindX(xB)];
345 phiStart = minXf-180; //open ring or acceptance effect...so believe to min phi angle!
347 phiStart = 0.5*(maxXf-180+minXf);
348 xA = sphi->Eval(phiStart+180);
349 if(xA < 0) xA = dist[sphi->FindX(xA)];
350 xB = sphi->Eval(phiStart+90);
351 if(xB < 0) xB = dist[sphi->FindX(xB)];
354 // Printf("ShapeModel: phiStart %f xA %f xB %f",phiStart,xA,xB);
357 phiStart*=TMath::DegToRad();
358 Double_t phitest = FindSimmPhi();
359 //Double_t phiStart = phitest;
361 Printf(" started phi %6.2f ",phiStart*TMath::RadToDeg());
363 // Printf("phiStart %f phiTest %f",phiStart*TMath::RadToDeg(),phitest*TMath::RadToDeg());
367 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
368 Double_t AliHMPIDReconHTA::VertParab(Double_t x1,Double_t y1,Double_t x2, Double_t y2, Double_t x3, Double_t y3)const
370 // It uses parabola from 3 points to evaluate the x-coord of the parab
371 // Arguments: xi,yi - points
372 // Returns: x-coord of the vertex
374 Double_t a = ((x1-x3)*(y1-y2)-(x1-x2)*(y1-y3))/((x1*x1-x2*x2)*(x1-x3)-(x1*x1-x3*x3)*(x1-x2));
375 Double_t b = (y1-y2 - a*(x1*x1-x2*x2))/(x1-x2);
376 // Double_t c = y1 - a*x1*x1-b*x1;
379 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
380 Bool_t AliHMPIDReconHTA::FitFree(Double_t thTrkRec,Double_t phiTrkRec)
382 // Fit performed by minimizing RMS/sqrt(n) of the
383 // photons reconstructed. First phi is fixed and theta
384 // is fouond, then (th,ph) of the track
385 // as free parameters
386 // Arguments: PhiRec phi of the track
389 TMinuit *pMinuit = new TMinuit(2);
391 gMinuit->SetObjectFit((TObject*)this); gMinuit->SetFCN(AliHMPIDReconHTA::FunMinPhot); //set fit function
392 Double_t aArg=-1,parStep,parLow,parHigh; Int_t iErrFlg; //tmp vars for TMinuit
397 pMinuit->mnexcm("SET PRI",&aArg,1,iErrFlg); //suspend all printout from TMinuit
398 pMinuit->mnexcm("SET NOW",&aArg,0,iErrFlg);
400 if(thTrkRec==0) thTrkRec = 3.*TMath::DegToRad(); // not to start from the edge...
402 pMinuit->mnparm(0," thTrk ",thTrkRec ,parStep=0.01,parLow=0,parHigh=TMath::PiOver4(),iErrFlg);
403 pMinuit->mnparm(1," phiTrk ",phiTrkRec,parStep=0.01,parLow=0,parHigh=TMath::TwoPi(),iErrFlg);
405 pMinuit->FixParameter(1);
406 pMinuit->mnexcm("SIMPLEX" ,&aArg,0,iErrFlg);
407 pMinuit->mnexcm("MIGRAD" ,&aArg,0,iErrFlg);
409 pMinuit->mnexcm("MIGRAD" ,&aArg,0,iErrFlg);
411 pMinuit->mnpout(0,sName,th,d1,d2,d3,iErrFlg);
412 pMinuit->mnpout(1,sName,ph,d1,d2,d3,iErrFlg);
416 par[0] = th;par[1] = ph;
417 pMinuit->Eval(2,grad,f,par,3);
419 // Printf("FitFree: theta %f phi %f",th,ph);
424 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
425 void AliHMPIDReconHTA::FunMinPhot(Int_t &/* */,Double_t* /* */,Double_t &f,Double_t *par,Int_t iflag)
427 // Minimization function to find best track and thetaC parameters
428 // Arguments: f = function value to minimize
429 // par = list of parameter to find
430 // iflag = flag status. See Minuit instructions
433 // Note: it is necessary to call an instance of AlihMPIDParam. Not possible to use fParam
434 // because of the static instantiation of the function in Minuit
436 AliHMPIDParam *pParam=AliHMPIDParam::Instance();
437 AliHMPIDReconHTA *pRecHTA=(AliHMPIDReconHTA*)gMinuit->GetObjectFit();
439 Double_t sizeCh = 0.5*pParam->RadThick()+pParam->WinThick()+pParam->GapThick();
440 Double_t thTrk = par[0];
441 Double_t phTrk = par[1];
442 Double_t xrad = pRecHTA->MipX() - sizeCh*TMath::Tan(thTrk)*TMath::Cos(phTrk);
443 Double_t yrad = pRecHTA->MipY() - sizeCh*TMath::Tan(thTrk)*TMath::Sin(phTrk);
444 pRecHTA->SetRadXY(xrad,yrad);
445 pRec.SetTrack(xrad,yrad,thTrk,phTrk);
447 Double_t meanCkov =0;
448 Double_t meanCkov2=0;
449 Double_t thetaCer,phiCer;
451 Int_t nClTot=pRecHTA->NClu();
453 for(Int_t i=0;i<nClTot;i++) {
454 if(!(pRecHTA->ClCk(i))) continue;
455 pRec.FindPhotCkov(pRecHTA->XClu(i),pRecHTA->YClu(i),thetaCer,phiCer);
456 meanCkov += thetaCer;
457 meanCkov2 += thetaCer*thetaCer;
460 if(nClAcc==0) {f=999;return;}
461 meanCkov /=(Double_t)nClAcc;
462 meanCkov2 /=(Double_t)nClAcc;
463 Double_t rms = TMath::Sqrt(TMath::Abs(meanCkov2 - meanCkov*meanCkov));
464 f = rms/TMath::Sqrt((Double_t)nClAcc);
468 Printf("FunMinPhot before: photons candidates %i used %i",nClTot,nClAcc);
470 Double_t meanCkov1=0;
471 Double_t meanCkov2=0;
472 for(Int_t i=0;i<nClTot;i++) {
473 if(!(pRec->ClCk(i))) continue;
474 pRec->FindPhotCkov(pRec->XClu(i),pRec->YClu(i),thetaCer,phiCer);
475 if(TMath::Abs(thetaCer-meanCkov)<2*rms) {
476 meanCkov1 += thetaCer;
477 meanCkov2 += thetaCer*thetaCer;
479 } else pRec->SetClCk(i,kFALSE);
482 Double_t rms2 = (meanCkov2 - meanCkov*meanCkov*nClAcc)/nClAcc;
483 Printf("FunMinPhot after: photons candidates %i used %i thetaCer %f",nClTot,nClAcc,meanCkov1);
484 pRec->SetCkovFit(meanCkov1);
485 pRec->SetCkovSig2(rms2);
486 pRec->SetNClu(nClAcc);
488 // Printf("FunMinPhot: reconstructed theta Cerenkov %f with %d photons",meanCkov,nClAcc);
489 pRecHTA->SetCkovFit(meanCkov);
490 pRecHTA->SetCkovSig2(rms*rms);
491 pRecHTA->SetNClu(nClAcc);
494 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
495 void AliHMPIDReconHTA::InitDatabase()
497 // Construction a database of ring shapes on fly
500 // N.B. fgDatabase points to a TH2I with x-min dist from MIP
501 // y-dist from the ring of the MIP perpendicular to major axis
502 // The content is the packed info of track theta and thetaC in degrees
505 // TFile *pout = new TFile("./database.root","recreate");
510 AliInfo(Form("database HTA is being built.Please, wait..."));
512 Double_t x[3]={0,0,0},y[3];
516 if(!fParam) fParam=AliHMPIDParam::Instance();
517 Double_t thetaMax = TMath::ACos(1./fParam->MeanIdxRad());
518 Double_t thTrkMax = 1./TMath::ASin(fParam->MeanIdxRad());
523 // TH2F *fgDatabase = new TH2F("deconv","database;d1;d2;thC+1000*thTrk",500,0,50,150,0,15);
529 for(Int_t i=0;i<nstepx;i++) { //loop on thetaC
530 for(Int_t j=0;j<nstepy;j++) { //loop on theta particle
531 Double_t thetaC = thetaMax/nstepx*((Double_t)i+0.5);
532 Double_t thTrk = thTrkMax/nstepy*((Double_t)j+0.5);
536 Double_t sizeCh = 0.5*fParam->RadThick()+fParam->WinThick()+fParam->GapThick();
537 Double_t xmip = xrad + sizeCh*TMath::Tan(thTrk)*TMath::Cos(phTrk);
538 Double_t ymip = yrad + sizeCh*TMath::Tan(thTrk)*TMath::Sin(phTrk);
540 Double_t dist1,dist2;
542 //first point at phi=0
544 rec.SetTrack(xrad,yrad,thTrk,phTrk);
546 pos=rec.TracePhot(thetaC,0);
549 dist1 = 0; //open ring...anly the distance btw mip and point at 180 will be considered
551 x[0] = pos.X(); y[0] = pos.Y();
552 dist1 = TMath::Sqrt((x[0]-xmip)*(x[0]-xmip)+(y[0]-ymip)*(y[0]-ymip));
555 //second point at phi=180
557 rec.SetTrack(xrad,yrad,thTrk,phTrk);
558 pos=rec.TracePhot(thetaC,TMath::Pi());
560 if(pos.X()==-999) {Printf("it should not happens!Bye");return;}
561 x[1] = pos.X(); y[1] = pos.Y();
562 if((x[1]-xmip)*(x[0]-xmip)>0) continue; // to avoid circles out mips (for very low ThetaC)
563 dist2 = TMath::Sqrt((x[1]-xmip)*(x[1]-xmip)+(y[1]-ymip)*(y[1]-ymip));
565 // Double_t distA = dist1+dist2;
566 Double_t distA = dist2; // only the minimum: problem of acceptance
568 //second point at phi=90
570 rec.SetTrack(xrad,yrad,thTrk,phTrk);
571 pos=rec.TracePhot(thetaC,TMath::PiOver2());
573 if(pos.X()==-999) continue;
574 x[2] = pos.X(); y[2] = pos.Y();
575 Double_t distB = TMath::Sqrt((x[2]-xmip)*(x[2]-xmip)+(y[2]-ymip)*(y[2]-ymip));
576 // compact the infos...
577 Int_t compact = (Int_t)(thetaC*TMath::RadToDeg())+1000*(Int_t)(thTrk*TMath::RadToDeg());
578 Int_t bin = fgDatabase->FindBin(distA,distB);
579 if(fgDatabase->GetBinContent(bin)==0) fgDatabase->Fill(distA,distB,compact);
584 // fgDatabase = deconv;
587 Double_t nSecs = timer.CpuTime();
588 AliInfo(Form("database HTA successfully open in %3.1f sec.(CPU). Reconstruction is started.",nSecs));
594 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
595 void AliHMPIDReconHTA::FillZeroChan()const
597 //If fills eventually channel without entries
598 //inthe histo "database" jyst interpolating the neighboring cells
599 // Arguments: histogram pointer of the database
602 Int_t nbinx = fgDatabase->GetNbinsX();
603 Int_t nbiny = fgDatabase->GetNbinsY();
604 for(Int_t i = 0;i<nbinx;i++) {
605 for(Int_t j = 0;j<nbiny;j++) {
606 if(fgDatabase->GetBinContent(i,j) == 0) {
607 Int_t nXmin = i-1; Int_t nXmax=i+1;
608 Int_t nYmin = j-1; Int_t nYmax=j+1;
612 for(Int_t ix=nXmin;ix<=nXmax;ix++) {
613 if(ix<0||ix>nbinx) continue;
614 for(Int_t iy=nYmin;iy<=nYmax;iy++) {
615 if(iy<0||iy>nbiny) continue;
616 meanC += (Int_t)fgDatabase->GetBinContent(ix,iy)%1000;
617 meanTrk+= (Int_t)fgDatabase->GetBinContent(ix,iy)/1000;
620 meanC/=nc; meanTrk/=nc;
621 Int_t compact = (Int_t)meanC+1000*(Int_t)meanTrk;
622 if(compact>0)fgDatabase->SetCellContent(i,j,compact);
628 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
630 // stima gli angoli con il metodo dei minimi quadrati che sfrutta le distanze...
632 Double_t AliHMPIDReconHTA::FindSimmPhi()
634 //It finds the phi of the ring
635 //by using the min. dist. algorithm
640 Float_t xrotsumm =0; Float_t yrotsumm =0; Float_t xx =0;
641 Float_t yy =0; Float_t xy =0;
645 for(Int_t i=0;i<fNClu;i++) {
646 if(!fClCk[i]) continue;
648 xrotsumm+=fXClu[i]; // summ xi
649 yrotsumm+=fYClu[i]; // summ yi
650 xx+=fXClu[i]*fXClu[i]; // summ xixi
651 yy+=fYClu[i]*fYClu[i]; // summ yiyi
652 xy+=fXClu[i]*fYClu[i]; // summ yixi
655 //_____calc. met min quadr using effective distance _________________________________________________
658 coeff[0]= xy-xrotsumm*yrotsumm/np;
659 coeff[1]= yrotsumm*yrotsumm/np - xrotsumm*xrotsumm/np - yy + xx;
660 coeff[2]= xrotsumm*yrotsumm/np- xy;
662 //____________________________________________________________________________________________________
669 n1=(yrotsumm-m1*xrotsumm)/np;
670 n2=(yrotsumm-m2*xrotsumm)/np;
672 // le due soluzioni.................
674 Double_t d1 =(1/(m1*m1+1))*(yy+m1*m1*xx+np*n1*n1-2*m1*xy-2*n1*yrotsumm+2*m1*n1*xrotsumm);
675 Double_t d2 =(1/(m2*m2+1))*(yy+m2*m2*xx+np*n2*n2-2*m2*xy-2*n2*yrotsumm+2*m2*n2*xrotsumm);
677 if(d1 > d2) mMin = m2; else mMin = m1;
679 Double_t PhiTrk= TMath::ATan(mMin);
681 // positive angles...
682 if(PhiTrk<0) PhiTrk+=180*TMath::DegToRad();
687 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
688 Int_t AliHMPIDReconHTA::r2(Double_t *coef, Double_t &x1, Double_t &x2)
692 // Arguments: coef 2 1 0: ax^2+bx+c=0
693 // Returns: n. of solutions
700 Double_t delta = b*b-4*a*c;
701 if(delta<0) {return 0;}
707 x1 = (-b+TMath::Sqrt(delta))/(2*a);
708 x2 = (-b-TMath::Sqrt(delta))/(2*a);