From 35e721bfde0f8d9d8c121a01f1c1b61154b14ce6 Mon Sep 17 00:00:00 2001 From: nick Date: Fri, 22 Sep 2006 11:18:53 +0000 Subject: [PATCH] 19-sep-2006 NvE Accuracy problem solved in Ali3Vector::GetOpeningAngle. 20-sep-2006 NvE AliSample::GetMedian corrected for empty or single value arrays. 22-sep-2006 NvE IceDwalk extended for specification of maximum number of hits per OM, to allow comparison with Sieglinde (which only uses the first OM hit). Also IceDwalk now introduces a device into the event named "IceDwalk" which contains the values of all the used reconstruction parameters. --- RALICE/Ali3Vector.cxx | 2 + RALICE/AliSample.cxx | 109 +++++++++++++++++++++++++++++--- RALICE/AliSample.h | 2 + RALICE/history.txt | 2 + RALICE/icepack/IceDwalk.cxx | 122 +++++++++++++++++++++++++++--------- RALICE/icepack/IceDwalk.h | 8 ++- RALICE/icepack/history.txt | 4 ++ 7 files changed, 208 insertions(+), 41 deletions(-) diff --git a/RALICE/Ali3Vector.cxx b/RALICE/Ali3Vector.cxx index df665121699..8645655e4c5 100644 --- a/RALICE/Ali3Vector.cxx +++ b/RALICE/Ali3Vector.cxx @@ -1045,6 +1045,8 @@ Double_t Ali3Vector::GetOpeningAngle(Ali3Vector& q,TString u) Double_t x=v1.Dot(v2); Double_t dx=fDresult; + if (x>1.) x=1; + if (x<-1.) x=-1; ang=acos(x); fDresult=0; if (fabs(x)<1.-dx) fDresult=dx/sqrt(1.-x*x); diff --git a/RALICE/AliSample.cxx b/RALICE/AliSample.cxx index 5a1abc89068..bd491681b6a 100644 --- a/RALICE/AliSample.cxx +++ b/RALICE/AliSample.cxx @@ -510,7 +510,13 @@ void AliSample::Data() cout << " " << fNames[i] << " : N = " << fN; cout << " Sum = " << fSum[i] << " Mean = " << fMean[i]; cout << " Var = " << fVar[i] << " Sigma = " << fSigma[i]; - if (fStore) cout << " Median = " << GetMedian(i+1); + if (fStore) + { + cout << endl; + cout << " Minimum = " << GetMinimum(i+1); + cout << " Maximum = " << GetMaximum(i+1); + cout << " Median = " << GetMedian(i+1); + } cout << endl; } } @@ -527,7 +533,13 @@ void AliSample::Data(Int_t i) cout << " " << fNames[i-1] << " : N = " << fN; cout << " Sum = " << fSum[i-1] << " Mean = " << fMean[i-1]; cout << " Var = " << fVar[i-1] << " Sigma = " << fSigma[i-1]; - if (fStore) cout << " Median = " << GetMedian(i); + if (fStore) + { + cout << endl; + cout << " Minimum = " << GetMinimum(i); + cout << " Maximum = " << GetMaximum(i); + cout << " Median = " << GetMedian(i); + } cout << endl; } } @@ -561,8 +573,8 @@ void AliSample::SetStoreMode(Int_t mode) // For normal statistics evaluation (e.g. mean, sigma, covariance etc...) // storage of entered data is not needed. This is the default mode // of operation and is the most efficient w.r.t. cpu time and memory. -// However, when calculation of a median is required, then the data -// storage mode has be activated. +// However, when calculation of a median, minimum or maximum is required, +// then the data storage mode has be activated. // // Note : Activation of storage mode can only be performed before the // first data item is entered. @@ -585,11 +597,12 @@ Int_t AliSample::GetStoreMode() const /////////////////////////////////////////////////////////////////////////// Float_t AliSample::GetMedian(Int_t i) { -// Provide the median of a certain variable +// Provide the median of a certain variable. +// For this functionality the storage mode has to be activated. if (fDim < i) { - cout << " *AliSample::mean* Error : Dimension less than " << i << endl; + cout << " *AliSample::GetMedian* Error : Dimension less than " << i << endl; return 0; } @@ -599,6 +612,18 @@ Float_t AliSample::GetMedian(Int_t i) return 0; } + if (fN<=0) return 0; + + Float_t median=0; + + if (fN==1) + { + if (i==1) median=fX->At(0); + if (i==2) median=fY->At(0); + if (i==3) median=fZ->At(0); + return median; + } + // Prepare temp. array to hold the ordered values if (!fArr) { @@ -646,7 +671,7 @@ Float_t AliSample::GetMedian(Int_t i) } } - Float_t median=0; + median=0; Int_t index=fN/2; if (fN%2) // Odd number of entries { @@ -659,3 +684,73 @@ Float_t AliSample::GetMedian(Int_t i) return median; } /////////////////////////////////////////////////////////////////////////// +Float_t AliSample::GetMinimum(Int_t i) const +{ +// Provide the minimum value of a certain variable. +// For this functionality the storage mode has to be activated. + + if (fDim < i) + { + cout << " *AliSample::GetMinimum* Error : Dimension less than " << i << endl; + return 0; + } + + if (!fStore) + { + cout << " *AliSample::GetMinimum* Error : Storage of data entries was not activated." << endl; + return 0; + } + + if (fN<=0) return 0; + + Float_t min=0; + + if (i==1) min=fX->At(0); + if (i==2) min=fY->At(0); + if (i==3) min=fZ->At(0); + + for (Int_t k=1; kAt(k)At(k); + if (i==2 && fY->At(k)At(k); + if (i==3 && fZ->At(k)At(k); + } + + return min; +} +/////////////////////////////////////////////////////////////////////////// +Float_t AliSample::GetMaximum(Int_t i) const +{ +// Provide the maxmum value of a certain variable. +// For this functionality the storage mode has to be activated. + + if (fDim < i) + { + cout << " *AliSample::GetMaximum* Error : Dimension less than " << i << endl; + return 0; + } + + if (!fStore) + { + cout << " *AliSample::GetMaximum* Error : Storage of data entries was not activated." << endl; + return 0; + } + + if (fN<=0) return 0; + + Float_t max=0; + + if (i==1) max=fX->At(0); + if (i==2) max=fY->At(0); + if (i==3) max=fZ->At(0); + + for (Int_t k=1; kAt(k)>max) max=fX->At(k); + if (i==2 && fY->At(k)>max) max=fY->At(k); + if (i==3 && fZ->At(k)>max) max=fZ->At(k); + } + + return max; +} +/////////////////////////////////////////////////////////////////////////// diff --git a/RALICE/AliSample.h b/RALICE/AliSample.h index 9f1ca69984b..b7197963547 100644 --- a/RALICE/AliSample.h +++ b/RALICE/AliSample.h @@ -31,6 +31,8 @@ class AliSample Float_t GetCov(Int_t i, Int_t j) const; // Covariance for i-th and j-th variable Float_t GetCor(Int_t i, Int_t j) const; // Correlation for i-th and j-th variable Float_t GetMedian(Int_t i); // Provide median for i-th variable + Float_t GetMinimum(Int_t i) const; // Provide the minimum value for i-th variable + Float_t GetMaximum(Int_t i) const; // Provide the maximum value for i-th variable void Data(); // Stat. info for the complete sample void Data(Int_t i); // Stat. info for the i-th variable void Data(Int_t i, Int_t j) const; // Stat. info for i-th and j-th variable diff --git a/RALICE/history.txt b/RALICE/history.txt index ddb2598d699..602b4352d77 100644 --- a/RALICE/history.txt +++ b/RALICE/history.txt @@ -705,3 +705,5 @@ The modified AliAttrib and AliSignal don't rely on a GetSize() anymore, but determine the number of stored values from the actual array contents. +19-sep-2006 NvE Accuracy problem solved in Ali3Vector::GetOpeningAngle. +20-sep-2006 NvE AliSample::GetMedian corrected for empty or single value arrays. diff --git a/RALICE/icepack/IceDwalk.cxx b/RALICE/icepack/IceDwalk.cxx index 9da005a72b0..e71bf6b4245 100644 --- a/RALICE/icepack/IceDwalk.cxx +++ b/RALICE/icepack/IceDwalk.cxx @@ -32,7 +32,12 @@ // This allows selection of events for processing with a certain maximum and/or // minimum number of good Amanda OMs firing. // By default the minimum and maximum are set to 0 and 999, respectively, -// in the constructor, which implies no multiplicity selection. +// in the constructor, which implies no multiplicity selection. +// The maximum number of hits per Amanda OM to be used for the reconstruction +// can be specified via the memberfunction SetMaxHitsA(). +// By default all the hits of an OM are used, but this may lead to large +// processing time in case many noise and/or afterpulse signals are not +// recognised by the hit cleaning procedure. // The various reconstruction steps are summarised as follows : // // 1) Construction of track elements (TE's). @@ -180,6 +185,7 @@ IceDwalk::IceDwalk(const char* name,const char* title) : TTask(name,title) fRjdmax=fDmin; fMaxmodA=999; fMinmodA=0; + fMaxhitsA=0; fTrackname="IceDwalk"; fCharge=0; } @@ -331,6 +337,22 @@ void IceDwalk::SetMinModA(Int_t nmin) fMinmodA=nmin; } /////////////////////////////////////////////////////////////////////////// +void IceDwalk::SetMaxHitsA(Int_t nmax) +{ +// Set the maximum number of hits per Amanda module to be processed. +// +// Special values : +// nmax = 0 : No maximum limit set; all available hits will be processed +// < 0 : No hits will be processed +// +// In case the user selects a maximum number of hits per module, all the +// hits of each module will be ordered w.r.t. increasing hit time (LE). +// This allows selection of processing e.g. only the first hits etc... +// By default the maximum number of hits per Amanda modules is set to 0 in the ctor, +// which implies just processing all available hits without any maximum limit. + fMaxhitsA=nmax; +} +/////////////////////////////////////////////////////////////////////////// void IceDwalk::SetTrackName(TString s) { // Set (alternative) name identifier for the produced first guess tracks. @@ -362,6 +384,37 @@ void IceDwalk::Exec(Option_t* opt) IceEvent* evt=(IceEvent*)parent->GetObject("IceEvent"); if (!evt) return; + // Enter the reco parameters as a device in the event + AliSignal params; + params.SetNameTitle("IceDwalk","IceDwalk reco parameters"); + params.SetSlotName("Dmin",1); + params.SetSlotName("Dtmarg",2); + params.SetSlotName("Tangmax",3); + params.SetSlotName("Rtangmax",4); + params.SetSlotName("Rtdmax",5); + params.SetSlotName("Jangmax",6); + params.SetSlotName("Rjangmax",7); + params.SetSlotName("Rjdmax",8); + params.SetSlotName("MaxmodA",9); + params.SetSlotName("MinmodA",10); + params.SetSlotName("MaxhitsA",11); + + params.SetSignal(fDmin,1); + params.SetSignal(fDtmarg,2); + params.SetSignal(fTangmax,3); + params.SetSignal(fRtangmax,4); + params.SetSignal(fRtdmax,5); + params.SetSignal(fJangmax,6); + params.SetSignal(fRjangmax,7); + params.SetSignal(fRjdmax,8); + params.SetSignal(fMaxmodA,9); + params.SetSignal(fMinmodA,10); + params.SetSignal(fMaxhitsA,11); + + evt->AddDevice(params); + + if (fMaxhitsA<0) return; + // Fetch all fired Amanda OMs for this event TObjArray* aoms=evt->GetDevices("IceAOM"); Int_t naoms=aoms->GetEntries(); @@ -378,20 +431,13 @@ void IceDwalk::Exec(Option_t* opt) } if (ngoodfMaxmodA) return; - const Float_t c=0.3; // Light speed in vacuum in meters per ns - const Float_t nice=1.33; // Refractive index of ice + const Float_t c=0.299792; // Light speed in vacuum in meters per ns + const Float_t nice=1.35634; // Refractive index of ice const Float_t thetac=acos(1./nice); // Cherenkov angle (in radians) - // Storage of track elements with various time difference margins. - // temap(i,j) holds the i-th track element (TE) with a time difference margin - // of less than j*3 nanoseconds. Currently we use a maximum margin of 30 ns. + // Storage of track elements. TObjArray tes; tes.SetOwner(); - AliObjMatrix temap; - - // Counter of TEs for each 3 ns margin slot - TArrayI ntes(fDtmarg/3); - if (ntes.GetSize()==0) ntes.Set(1); AliPosition r1; AliPosition r2; @@ -405,14 +451,14 @@ void IceDwalk::Exec(Option_t* opt) AliSignal* sx2=0; Float_t dist=0; Float_t t1,t2,dt,t0; - Float_t dtmax,dttest; + Float_t dtmax; TObjArray hits; + TObjArray* ordered; // Check the hits of Amanda OM pairs for posible track elements. // Also all the good hits are stored in the meantime (to save CPU time) // for hit association with the various track elements lateron. AliTrack* te=0; - Int_t ite=0; for (Int_t i1=0; i1At(i1); @@ -421,14 +467,27 @@ void IceDwalk::Exec(Option_t* opt) r1=omx1->GetPosition(); // Select all the good hits of this first OM hits1.Clear(); + // Determine the max. number of hits to be processed for this OM + ordered=0; + if (fMaxhitsA>0 && omx1->GetNhits()>fMaxhitsA) ordered=omx1->SortHits("LE",1,0,7); + nh1=0; for (Int_t j1=1; j1<=omx1->GetNhits(); j1++) { - sx1=omx1->GetHit(j1); + if (ordered) + { + if (nh1>=fMaxhitsA) break; + sx1=(AliSignal*)ordered->At(j1-1); + } + else + { + sx1=omx1->GetHit(j1); + } if (!sx1) continue; if (sx1->GetDeadValue("ADC") || sx1->GetDeadValue("LE") || sx1->GetDeadValue("TOT")) continue; hits1.Add(sx1); // Also store all good hits in the total hit array hits.Add(sx1); + nh1++; } // No further pair to be formed with the last OM in the list @@ -450,12 +509,25 @@ void IceDwalk::Exec(Option_t* opt) // Select all the good hits of this second OM hits2.Clear(); + // Determine the max. number of hits to be processed for this OM + ordered=0; + if (fMaxhitsA>0 && omx2->GetNhits()>fMaxhitsA) ordered=omx2->SortHits("LE",1,0,7); + nh2=0; for (Int_t j2=1; j2<=omx2->GetNhits(); j2++) { - sx2=omx2->GetHit(j2); + if (ordered) + { + if (nh2>=fMaxhitsA) break; + sx2=(AliSignal*)ordered->At(j2-1); + } + else + { + sx2=omx2->GetHit(j2); + } if (!sx2) continue; if (sx2->GetDeadValue("ADC") || sx2->GetDeadValue("LE") || sx2->GetDeadValue("TOT")) continue; hits2.Add(sx2); + nh2++; } nh2=hits2.GetEntries(); @@ -485,21 +557,12 @@ void IceDwalk::Exec(Option_t* opt) te=new AliTrack(); tes.Add(te); - ite++; if (dt<0) r12*=-1.; r0.SetTimestamp((AliTimestamp&)*evt); AliTimestamp* tsx=r0.GetTimestamp(); tsx->Add(0,0,(int)t0); te->SetReferencePoint(r0); te->Set3Momentum(r12); - dttest=dtmax; - for (Int_t jt=ntes.GetSize(); jt>0; jt--) - { - if (fabs(dt)>=dttest) break; - temap.EnterObject(ite,jt,te); - ntes.AddAt(ntes.At(jt-1)+1,jt-1); - dttest-=3.; - } } } } // end of loop over the second OM of the pair @@ -537,8 +600,7 @@ void IceDwalk::Exec(Option_t* opt) IceGOM* omx=(IceGOM*)sx1->GetDevice(); if (!omx) continue; r1=omx->GetPosition(); - d=tr0->GetDistance(r1); - d*=sin(thetac); + d=te->GetDistance(r1); r12=r1-(*tr0); dist=p.Dot(r12)+d*tan(thetac); tgeo=t0+dist/c; @@ -574,7 +636,6 @@ void IceDwalk::Exec(Option_t* opt) if (sx1) qtc=sx1->GetSignal("QTC"); if (qtc<0.7*qmax) { - temap.RemoveObjects(te); tes.RemoveAt(jtc); delete te; } @@ -585,7 +646,7 @@ void IceDwalk::Exec(Option_t* opt) if (!nte) return; // Order the track candidates w.r.t. decreasing number of associated hits - TObjArray* ordered=0; + ordered=0; ordered=evt->SortTracks(-1,&tes); TObjArray tcs(*ordered); @@ -604,7 +665,6 @@ void IceDwalk::Exec(Option_t* opt) AliSample pos; AliSample time; Float_t vec[3],err[3]; - Float_t edist=0; for (Int_t jtc1=0; jtc1GetTimestamp(); if (!ts2) continue; dist=x1->GetDistance(x2); - edist=x1->GetResultError(); dt=ts1->GetDifference(ts2,"ns"); if (dist>0) { r12=(*x2)-(*x1); if (dt<0) r12*=-1.; ang=te->GetOpeningAngle(r12,"deg"); - if (angGetPosition(vec,"car"); pos.Enter(vec[0],vec[1],vec[2]); @@ -698,6 +757,7 @@ void IceDwalk::Exec(Option_t* opt) // Merge jets within a certain opening to provide the final track(s). AliJet* jx1=0; AliJet* jx2=0; + Float_t edist=0; if (fJangmax>0) { for (Int_t jet1=0; jet1