X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=RALICE%2Ficepack%2FIceDwalk.cxx;h=296ce99f682404abbd1ef0bb3d93dc70b8a7691e;hb=caa58e1ad7e2710be3a39cacc822c19108b1cf53;hp=683601553420687535ed8e3c34e6999c1c661890;hpb=8a394508e6b1df25049071497914c3ba97843c85;p=u%2Fmrichter%2FAliRoot.git diff --git a/RALICE/icepack/IceDwalk.cxx b/RALICE/icepack/IceDwalk.cxx index 68360155342..296ce99f682 100644 --- a/RALICE/icepack/IceDwalk.cxx +++ b/RALICE/icepack/IceDwalk.cxx @@ -18,6 +18,10 @@ /////////////////////////////////////////////////////////////////////////// // Class IceDwalk // TTask derived class to perform direct walk track reconstruction. +// +// In case an event has been rejected by an AliEventSelector (based) processor, +// this task (and its sub-tasks) is not executed. +// // The procedure is based on the method described in the Amanda publication // in Nuclear Instruments and Methods A524 (2004) 179-180. // However, the Amanda method has been extended with the intention to @@ -79,10 +83,10 @@ // tres : time residual // Difference between the observed hit time and the time expected // for a direct photon hit. -// dhit : Distance between the hit and the TE +// dhit : Distance traveled by the cherenkov photon from the track to the hit position // lambda : Photon scattering length in ice // -// By default F is set to 2, but this can be modified via the memberfunction +// By default F is set to 3.07126, but this can be modified via the memberfunction // SetMaxDhit(). // // 3) Construction of track candidates (TC's). @@ -257,9 +261,10 @@ IceDwalk::IceDwalk(const char* name,const char* title) : TTask(name,title) // The newly introduced angular separation parameter for jet merging // is initialised as half the value of the angular separation parameter // for track candidate clustering. + fEvt=0; fDmin=75; fDtmarg=0; - fMaxdhit=2; + fMaxdhit=3.07126; fTangmax=15; fTdistmax=20; fTinvol=1; @@ -271,7 +276,7 @@ IceDwalk::IceDwalk(const char* name,const char* title) : TTask(name,title) fMinmodA=0; fMaxhitsA=1; fVgroup=1; - fTrackname="IceDwalk"; + fTrackname=""; fCharge=0; } /////////////////////////////////////////////////////////////////////////// @@ -439,8 +444,8 @@ void IceDwalk::SetTrackName(TString s) // Set (alternative) name identifier for the produced first guess tracks. // This allows unique identification of (newly) produced direct walk tracks // in case of re-processing of existing data with different criteria. -// By default the produced first guess tracks have the name "IceDwalk" -// which is set in the constructor of this class. +// By default the produced first guess tracks have the name of the class +// by which they were produced. fTrackname=s; } /////////////////////////////////////////////////////////////////////////// @@ -462,12 +467,19 @@ void IceDwalk::Exec(Option_t* opt) if (!parent) return; - IceEvent* evt=(IceEvent*)parent->GetObject("IceEvent"); - if (!evt) return; + fEvt=(IceEvent*)parent->GetObject("IceEvent"); + if (!fEvt) return; + + // Only process accepted events + AliDevice* seldev=(AliDevice*)fEvt->GetDevice("AliEventSelector"); + if (seldev) + { + if (seldev->GetSignal("Select") < 0.1) return; + } // Enter the reco parameters as a device in the event AliSignal params; - params.SetNameTitle("IceDwalk","IceDwalk reco parameters"); + params.SetNameTitle(ClassName(),"Reco parameters"); params.SetSlotName("Dmin",1); params.SetSlotName("Dtmarg",2); params.SetSlotName("Maxdhit",3); @@ -498,12 +510,12 @@ void IceDwalk::Exec(Option_t* opt) params.SetSignal(fMaxhitsA,13); params.SetSignal(fVgroup,14); - evt->AddDevice(params); + fEvt->AddDevice(params); if (fMaxhitsA<0) return; // Fetch all fired Amanda OMs for this event - TObjArray* aoms=evt->GetDevices("IceAOM"); + TObjArray* aoms=fEvt->GetDevices("IceAOM"); if (!aoms) return; Int_t naoms=aoms->GetEntries(); if (!naoms) return; @@ -519,16 +531,7 @@ void IceDwalk::Exec(Option_t* opt) } if (ngoodfMaxmodA) return; - const Float_t pi=acos(-1.); - const Float_t c=0.299792458; // Light speed in vacuum in meters per ns - const Float_t npice=1.31768387; // Phase refractive index (c/v_phase) of ice - const Float_t ngice=1.35075806; // Group refractive index (c/v_group) of ice - const Float_t lambda=33.3; // Light scattering length in ice - const Float_t thetac=acos(1./npice); // Cherenkov angle (in radians) - - // Angular reduction of complement of thetac due to v_phase and v_group difference - Float_t alphac=0; - if (fVgroup) alphac=atan((1.-npice/ngice)/sqrt(npice*npice-1.)); + const Float_t c=0.299792458; // Light speed in vacuum in meters per ns // Storage of track elements. TObjArray tes; @@ -544,13 +547,13 @@ void IceDwalk::Exec(Option_t* opt) Int_t nh1,nh2; AliSignal* sx1=0; AliSignal* sx2=0; - Float_t dist=0,dist2=0; + Float_t dist=0; Float_t t1,t2,dt,t0; Float_t dtmax; TObjArray hits; TObjArray* ordered; - // Check the hits of Amanda OM pairs for posible track elements. + // Check the hits of Amanda OM pairs for possible 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; @@ -653,7 +656,7 @@ void IceDwalk::Exec(Option_t* opt) te=new AliTrack(); tes.Add(te); if (dt<0) r12*=-1.; - r0.SetTimestamp((AliTimestamp&)*evt); + r0.SetTimestamp((AliTimestamp&)*fEvt); AliTimestamp* tsx=r0.GetTimestamp(); tsx->Add(0,0,(int)t0); te->SetReferencePoint(r0); @@ -664,11 +667,59 @@ void IceDwalk::Exec(Option_t* opt) } // end of loop over first OM of the pair // Association of hits to the various track elements. + Float_t qmax=0; + Int_t nahmax=0; + AssociateHits(tes,hits,qmax,nahmax); + + // Selection on quality (Q value) in case of multiple track candidates + SelectQvalue(tes,qmax); + + Int_t nte=tes.GetEntries(); + if (!nte) return; + + // Clustering of track candidates into jets + TObjArray jets; + jets.SetOwner(); + ClusterTracks(tes,jets,qmax); + + Int_t njets=jets.GetEntries(); + if (!njets) return; + + // Order the jets w.r.t. decreasing quality value + ordered=fEvt->SortJets(-2,&jets); + TObjArray jets2(*ordered); + + // Merging f jets + MergeJets(jets2); + + // Production and storage of the final tracks + StoreTracks(jets2); +} +/////////////////////////////////////////////////////////////////////////// +void IceDwalk::AssociateHits(TObjArray& tes,TObjArray& hits,Float_t& qmax,Int_t& nahmax) +{ + // Association of hits to the various track elements. + + const Float_t pi=acos(-1.); + const Float_t c=0.299792458; // Light speed in vacuum in meters per ns + const Float_t npice=1.31768387; // Phase refractive index (c/v_phase) of ice + const Float_t ngice=1.35075806; // Group refractive index (c/v_group) of ice + const Float_t lambda=33.3; // Light scattering length in ice + const Float_t thetac=acos(1./npice); // Cherenkov angle (in radians) + + // Angular reduction of complement of thetac due to v_phase and v_group difference + Float_t alphac=0; + if (fVgroup) alphac=atan((1.-npice/ngice)/sqrt(npice*npice-1.)); + Int_t nte=tes.GetEntries(); Int_t nh=hits.GetEntries(); Float_t d=0; Ali3Vector p; - Float_t tgeo,tres; + AliPosition r1; + AliPosition r2; + Ali3Vector r12; + Float_t t1; + Float_t dist,t0,tgeo,tres; AliSample levers; // Statistics of the assoc. hit lever arms levers.SetStoreMode(1);// Enable median calculation AliSample hprojs; // Statistics of the assoc. hit position projections on the track w.r.t. r0 @@ -698,26 +749,28 @@ void IceDwalk::Exec(Option_t* opt) fit.AddNamedSlot("term3"); fit.AddNamedSlot("term4"); fit.AddNamedSlot("term5"); - Float_t qtc=0,qmax=0; - Int_t nah,nahmax=0; // Number of associated hits for a certain TE + Float_t qtc=0; + Int_t nah; // Number of associated hits for a certain TE Float_t lmin,lmax,spanl,medianl,meanl,sigmal,spreadl,expspreadl; Float_t hproj,hprojmin,hprojmax,span,median,mean,sigma,spread,expspread; Float_t mediant,meant,sigmat,spreadt; Float_t term1,term2,term3,term4,term5; + qmax=0; + nahmax=0; for (Int_t jte=0; jteGetReferencePoint(); AliTimestamp* tt0=tr0->GetTimestamp(); - t0=evt->GetDifference(tt0,"ns"); + t0=fEvt->GetDifference(tt0,"ns"); p=te->Get3Momentum(); levers.Reset(); hprojs.Reset(); times.Reset(); for (Int_t jh=0; jhGetDevice(); if (!omx) continue; @@ -730,6 +783,8 @@ void IceDwalk::Exec(Option_t* opt) t1=sx1->GetSignal("LE",7); tres=t1-tgeo; + d=d/sin(thetac); // The distance traveled by a cherenkov photon + if (tres<-30 || tres>300 || d>fMaxdhit*lambda) continue; // Associate this hit to the TE @@ -811,14 +866,22 @@ void IceDwalk::Exec(Option_t* opt) te->SetFitDetails(fit); if (qtc>qmax) qmax=qtc; } - +} +/////////////////////////////////////////////////////////////////////////// +void IceDwalk::SelectQvalue(TObjArray& tes,Float_t qmax) +{ // Perform selection on Q value in case of multiple track candidates + + Int_t nte=tes.GetEntries(); + Int_t nah; + Float_t qtc; + Ali3Vector p; for (Int_t jtc=0; jtcGetNsignals(); - sx1=(AliSignal*)te->GetFitDetails(); + AliSignal* sx1=(AliSignal*)te->GetFitDetails(); qtc=-1; if (sx1) qtc=sx1->GetSignal("QTC"); if (!nah || qtc<0.8*qmax) @@ -837,10 +900,10 @@ void IceDwalk::Exec(Option_t* opt) } } tes.Compress(); - nte=tes.GetEntries(); - - if (!nte) return; - +} +/////////////////////////////////////////////////////////////////////////// +void IceDwalk::ClusterTracks(TObjArray& tes,TObjArray& jets,Float_t qmax) +{ // Cluster track candidates within a certain opening angle into jets. // Also the track should be within a certain maximum distance of the // starting track in order to get clustered. @@ -848,17 +911,19 @@ void IceDwalk::Exec(Option_t* opt) // crossing the detector a very different locations (e.g. muon bundles). // The average r0 and t0 of the constituent tracks will be taken as the // jet reference point. - TObjArray jets; - jets.SetOwner(); - AliTrack* te2=0; + + Int_t nte=tes.GetEntries(); Float_t ang=0; AliSample pos; AliSample time; Float_t vec[3],err[3]; - nahmax=0; // Determine the max. number of associated hits for the jets + AliPosition r0; + Float_t t0,dist,dist2; + Int_t nah=0,nahmax=0; // Determine the max. number of associated hits for the jets + Float_t qtc; for (Int_t jtc1=0; jtc1GetReferencePoint(); if (!x1) continue; @@ -870,12 +935,12 @@ void IceDwalk::Exec(Option_t* opt) time.Reset(); x1->GetPosition(vec,"car"); pos.Enter(vec[0],vec[1],vec[2]); - t0=evt->GetDifference(ts1,"ns"); + t0=fEvt->GetDifference(ts1,"ns"); time.Enter(t0); for (Int_t jtc2=0; jtc2GetOpeningAngle(*te2,"deg"); if (ang<=fTangmax) @@ -898,7 +963,7 @@ void IceDwalk::Exec(Option_t* opt) { x2->GetPosition(vec,"car"); pos.Enter(vec[0],vec[1],vec[2]); - t0=evt->GetDifference(ts2,"ns"); + t0=fEvt->GetDifference(ts2,"ns"); time.Enter(t0); jx->AddTrack(te2); } @@ -913,7 +978,7 @@ void IceDwalk::Exec(Option_t* opt) } r0.SetPosition(vec,"car"); r0.SetPositionErrors(err,"car"); - r0.SetTimestamp((AliTimestamp&)*evt); + r0.SetTimestamp((AliTimestamp&)*fEvt); AliTimestamp* jt0=r0.GetTimestamp(); t0=time.GetMean(1); jt0->Add(0,0,(int)t0); @@ -928,7 +993,7 @@ void IceDwalk::Exec(Option_t* opt) } else // Only keep single-track jets which have qtc=qmax { - sx1=(AliSignal*)te->GetFitDetails(); + AliSignal* sx1=(AliSignal*)te->GetFitDetails(); qtc=-1; if (sx1) qtc=sx1->GetSignal("QTC"); if (qtc>=(qmax-1.e-10)) @@ -943,11 +1008,11 @@ void IceDwalk::Exec(Option_t* opt) } } } - Int_t njets=jets.GetEntries(); + Int_t njets=jets.GetEntries(); if (!njets) return; - // The sum of 0.15*(nah-nahmax) and average qtc value per track for this jet + // The sum of 0.15*(nah-nahmax) and average qtc value per track for each jet // will be stored as the jet energy to enable sorting on this value lateron Float_t sortval=0; Int_t ntk=0; @@ -961,11 +1026,10 @@ void IceDwalk::Exec(Option_t* opt) if (ntk) sortval+=jx->GetMomentum()/float(ntk); jx->SetScalar(sortval); } - - // Order the jets w.r.t. decreasing value of the above sortval - ordered=evt->SortJets(-2,&jets); - TObjArray jets2(*ordered); - +} +/////////////////////////////////////////////////////////////////////////// +void IceDwalk::MergeJets(TObjArray& jets2) +{ // Merge jets within a certain opening to provide the final track(s). // Also the jet should be within a certain maximum distance of the // starting jet in order to get merged. @@ -973,9 +1037,18 @@ void IceDwalk::Exec(Option_t* opt) // crossing the detector a very different locations (e.g. muon bundles). // The average r0 and t0 of the constituent jets will be taken as the // final reference point. + + Int_t njets=jets2.GetEntries(); AliJet* jx1=0; AliJet* jx2=0; Int_t merged=1; + Int_t ntk,nah,nahmax; + Float_t ang,dist,dist2,t0; + AliSample pos; + AliSample time; + AliPosition r0; + Float_t vec[3],err[3]; + Float_t sortval; if (fJangmax>=0) { while (merged) @@ -994,7 +1067,7 @@ void IceDwalk::Exec(Option_t* opt) time.Reset(); x1->GetPosition(vec,"car"); pos.Enter(vec[0],vec[1],vec[2]); - t0=evt->GetDifference(ts1,"ns"); + t0=fEvt->GetDifference(ts1,"ns"); time.Enter(t0); for (Int_t jet2=0; jet2GetPosition(vec,"car"); pos.Enter(vec[0],vec[1],vec[2]); - t0=evt->GetDifference(ts2,"ns"); + t0=fEvt->GetDifference(ts2,"ns"); time.Enter(t0); for (Int_t jtk=1; jtk<=jx2->GetNtracks(); jtk++) { - te=jx2->GetTrack(jtk); + AliTrack* te=jx2->GetTrack(jtk); if (!te) continue; jx1->AddTrack(te); } @@ -1043,7 +1116,7 @@ void IceDwalk::Exec(Option_t* opt) } r0.SetPosition(vec,"car"); r0.SetPositionErrors(err,"car"); - r0.SetTimestamp((AliTimestamp&)*evt); + r0.SetTimestamp((AliTimestamp&)*fEvt); AliTimestamp* jt0=r0.GetTimestamp(); t0=time.GetMean(1); jt0->Add(0,0,(int)t0); @@ -1055,8 +1128,8 @@ void IceDwalk::Exec(Option_t* opt) jets2.Compress(); - // The sum of 0.15*(nah-nahmax) and average qtc value per track for this jet - // will be stored as the jet energy to enable sorting on this value lateron + // The sum of 0.15*(nah-nahmax) and average qtc value per track for each jet + // will be stored as the jet energy to enable sorting on this value for (Int_t jjet=0; jjetSetScalar(sortval); } - // Order the jets w.r.t. decreasing value of the above sortval - ordered=evt->SortJets(-2,&jets2); + // Order the jets w.r.t. decreasing quality value + TObjArray* ordered=fEvt->SortJets(-2,&jets2); njets=ordered->GetEntries(); jets2.Clear(); for (Int_t icopy=0; icopyGetReferencePoint(); if (!ref) continue; - evt->AddTrack(t); - AliTrack* trk=evt->GetTrack(evt->GetNtracks()); + fEvt->AddTrack(t); + AliTrack* trk=fEvt->GetTrack(fEvt->GetNtracks()); if (!trk) continue; - trk->SetId(evt->GetNtracks(1)+1); + trk->SetId(fEvt->GetNtracks(1)+1); p=jx->Get3Momentum(); p/=p.GetNorm(); trk->Set3Momentum(p); @@ -1114,7 +1197,7 @@ void IceDwalk::Exec(Option_t* opt) if (!tx) continue; for (Int_t is=1; is<=tx->GetNsignals(); is++) { - sx1=tx->GetSignal(is); + AliSignal* sx1=tx->GetSignal(is); if (sx1) sx1->AddTrack(*trk); } }