From c630aafd1da7fec82908f28924c025969497227d Mon Sep 17 00:00:00 2001 From: hristov Date: Tue, 22 Jul 2003 15:56:14 +0000 Subject: [PATCH] Implementing ESD functionality in the NewIO (Yu.Belikov) --- HBTAN/AliHBTReaderPPprod.cxx | 8 +- HBTAN/AliHBTReaderTPC.cxx | 7 +- ITS/AliITSComparisonV2.C | 96 +-- ITS/AliITSFindClustersV2.C | 266 ++----- ITS/AliITSFindTracksV2.C | 149 ++-- ITS/AliITSTrackerV1.cxx | 11 +- ITS/AliITSclustererV2.cxx | 112 ++- ITS/AliITSclustererV2.h | 9 +- ITS/AliITSpidESD.cxx | 2 +- ITS/AliITStestV2.C | 55 +- ITS/AliITStrackerV2.cxx | 449 +++++++---- ITS/AliITStrackerV2.h | 27 +- ITS/ITSLinkDef.h | 2 +- ITS/libITS.pkg | 1 + STEER/AliESDanalysis.C | 168 ++++- STEER/AliESDpid.cxx | 12 + STEER/AliESDtest.C | 233 ++++-- STEER/AliESDtrack.cxx | 33 + STEER/AliESDtrack.h | 18 +- STEER/AliTracker.h | 10 +- TOF/AliTOF.cxx | 8 +- TOF/AliTOFSDigitizer.cxx | 1 + TOF/AliTOFSDigits2Digits.C | 207 +++--- TOF/AliTOFpidESD.cxx | 451 +++++++++++ TOF/AliTOFpidESD.h | 80 ++ TOF/TOFLinkDef.h | 2 + TOF/libTOF.pkg | 2 +- TPC/AliTPC.cxx | 13 +- TPC/AliTPC.h | 2 +- TPC/AliTPCComparison.C | 63 +- TPC/AliTPCFindClusters.C | 127 ++-- TPC/AliTPCFindTracks.C | 101 +-- TPC/AliTPCHits2Digits.C | 50 +- TPC/AliTPCclusterer.cxx | 85 +-- TPC/AliTPCclusterer.h | 45 +- TPC/AliTPCpidESD.cxx | 2 +- TPC/AliTPCtest.C | 11 +- TPC/AliTPCtrack.cxx | 9 + TPC/AliTPCtrack.h | 1 + TPC/AliTPCtracker.cxx | 1358 ++++++++++++++++++++++++---------- TPC/AliTPCtracker.h | 89 ++- TPC/TPCLinkDef.h | 2 +- TRD/AliTRDclusterizer.cxx | 14 +- TRD/AliTRDdigits2cluster.C | 4 +- TRD/AliTRDtrack.cxx | 67 +- TRD/AliTRDtracker.cxx | 540 +++++++++++--- TRD/AliTRDtracker.h | 101 ++- 47 files changed, 3430 insertions(+), 1673 deletions(-) create mode 100644 TOF/AliTOFpidESD.cxx create mode 100644 TOF/AliTOFpidESD.h diff --git a/HBTAN/AliHBTReaderPPprod.cxx b/HBTAN/AliHBTReaderPPprod.cxx index 33ce438b836..6ab976238a4 100644 --- a/HBTAN/AliHBTReaderPPprod.cxx +++ b/HBTAN/AliHBTReaderPPprod.cxx @@ -163,14 +163,16 @@ Int_t AliHBTReaderPPprod::Read(AliHBTRun* particles, AliHBTRun *tracks) AliTPCtrack *iotrack=0; fClustersFile->cd(); - AliTPCtracker *tracker = new AliTPCtracker(TPCParam,currentEvent,""); +printf("This method is not converted to the NewIO !\n"); //I.B. +return 1; //I.B. + AliTPCtracker *tracker = new AliTPCtracker(TPCParam); //I.B. + //AliTPCtracker *tracker = new AliTPCtracker(TPCParam,currentEvent,""); if (!tracker) { Error("AliHBTReaderPPprod::Read","Can't get a tracker !\n"); return 3; } - tracker->LoadInnerSectors(); - tracker->LoadOuterSectors(); +tracker->LoadClusters(0);//I.Belikov, "0" should be a pointer to a tree for (i=0; iLoadInnerSectors(); - tracker->LoadOuterSectors(); + tracker->LoadClusters(0);//I.Belikov, "0" must be a pointer to a tree for (Int_t i=0; i #include #include "AliRun.h" + #include "AliHeader.h" + #include "AliStack.h" #include "AliRunLoader.h" #include "AliLoader.h" #include "AliITSLoader.h" @@ -32,41 +38,43 @@ struct GoodTrackITS { Float_t x,y,z; }; -Int_t AliITSComparisonV2(Int_t event=0) { +extern AliRun *gAlice; + +Int_t AliITSComparisonV2() { cerr<<"Doing comparison...\n"; - if (gAlice) - { + if (gAlice) { delete gAlice->GetRunLoader(); delete gAlice; gAlice=0; - } + } AliRunLoader *rl = AliRunLoader::Open("galice.root"); - if (!rl) - { + if (!rl) { cerr<<"AliITSComparisonV2.C :Can't start sesion !\n"; return 1; - } + } rl->LoadgAlice(); if (rl->GetAliRun()) - AliKalmanTrack::SetConvConst(1000/0.299792458/rl->GetAliRun()->Field()->SolenoidField()); - else - { - cerr<<"AliITSComparisonV2.C :Can't get AliRun !\n"; - return 1; - } + AliKalmanTrack:: + SetConvConst(1000/0.299792458/rl->GetAliRun()->Field()->SolenoidField()); + else { + cerr<<"AliITSComparisonV2.C :Can't get AliRun !\n"; + return 1; + } rl->UnloadgAlice(); AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader"); - if (itsl == 0x0) - { + if (itsl == 0x0) { cerr<<"AliITSComparisonV2.C : Can not find TPCLoader\n"; delete rl; return 3; - } + } const Int_t MAX=15000; - Int_t good_tracks_its(GoodTrackITS *gt, const Int_t max, const char* evfoldname = AliConfig::fgkDefaultEventFolderName);//declaration only + Int_t good_tracks_its( + GoodTrackITS *gt, const Int_t max, + const char* evfoldname = AliConfig::fgkDefaultEventFolderName + );//declaration only Int_t nentr=0; TObjArray tarray(2000); {/* Load tracks */ @@ -299,50 +307,40 @@ Int_t AliITSComparisonV2(Int_t event=0) { text->SetTextSize(0.05); text->Draw(); -// PH taken from v3-09-09, but not used -// TCanvas *c2=new TCanvas("c2","",320,32,530,590); - -// TPad *p6=new TPad("p6","",0.,0.,1.,.5); p6->Draw(); -// p6->cd(); p6->SetFillColor(42); p6->SetFrameFillColor(10); -// he->SetFillColor(2); he->SetFillStyle(3005); -// he->SetXTitle("Arbitrary Units"); -// he->Fit("gaus"); c2->cd(); - -// TPad *p7=new TPad("p7","",0.,0.5,1.,1.); p7->Draw(); -// p7->cd(); p7->SetFillColor(42); p7->SetFrameFillColor(10); -// hep->SetXTitle("p (Gev/c)"); hep->SetYTitle("dE/dX (Arb. Units)"); -// hep->Draw(); c1->cd(); + TCanvas *c2=new TCanvas("c2","",320,32,530,590); + TPad *p6=new TPad("p6","",0.,0.,1.,.5); p6->Draw(); + p6->cd(); p6->SetFillColor(42); p6->SetFrameFillColor(10); + he->SetFillColor(2); he->SetFillStyle(3005); + he->SetXTitle("Arbitrary Units"); + he->Fit("gaus"); c2->cd(); + + TPad *p7=new TPad("p7","",0.,0.5,1.,1.); p7->Draw(); + p7->cd(); p7->SetFillColor(42); p7->SetFrameFillColor(10); + hep->SetXTitle("p (Gev/c)"); hep->SetYTitle("dE/dX (Arb. Units)"); + hep->Draw(); c1->cd(); - return 0; } Int_t good_tracks_its(GoodTrackITS *gt, const Int_t max, const char* evfoldname) { - - Int_t nt=0; - AliRunLoader* rl = AliRunLoader::GetRunLoader(evfoldname); - if (rl == 0x0) - { + if (rl == 0x0) { ::Fatal("AliTPCComparison.C::good_tracks_its", "Can not find Run Loader in Folder Named %s", evfoldname); - } + } AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader"); - if (itsl == 0x0) - { + if (itsl == 0x0) { cerr<<"AliITSComparisonV2.C : Can not find TPCLoader\n"; delete rl; return 3; - } + } rl->LoadgAlice(); rl->LoadHeader(); Int_t np = rl->GetHeader()->GetNtrack(); - - Int_t *good=new Int_t[np]; Int_t k; for (k=0; kLoadRawClusters(); - TTree *cTree=itsl->TreeC(); + itsl->LoadRecPoints(); + TTree *cTree=itsl->TreeR(); if (!cTree) { cerr<<"Can't get cTree !\n"; exit(7); } @@ -379,8 +377,11 @@ Int_t good_tracks_its(GoodTrackITS *gt, const Int_t max, const char* evfoldname) while (ncl--) { AliITSclusterV2 *pnt=(AliITSclusterV2*)clusters->UncheckedAt(ncl); Int_t l0=pnt->GetLabel(0); + if (l0>=np) {cerr<<"Wrong label: "<GetLabel(1); + if (l1>=np) {cerr<<"Wrong label: "<GetLabel(2); + if (l2>=np) {cerr<<"Wrong label: "<=0) good[l0]|=mask; if (l1>=0) good[l1]|=mask; @@ -403,13 +404,12 @@ Int_t good_tracks_its(GoodTrackITS *gt, const Int_t max, const char* evfoldname) while (in>>lab>>code>>px>>py>>pz>>x>>y>>z) { if (good[lab] != 0x3F) continue; TParticle *p = (TParticle*)stack->Particle(lab); - if (p == 0x0) - { + if (p == 0x0) { cerr<<"Can not get particle "<GetPdgCode(); //**** px py pz - in global coordinate system diff --git a/ITS/AliITSFindClustersV2.C b/ITS/AliITSFindClustersV2.C index aba52f2e7a1..82808a0539e 100644 --- a/ITS/AliITSFindClustersV2.C +++ b/ITS/AliITSFindClustersV2.C @@ -1,223 +1,101 @@ -/**************************************************************** -* This macro converts AliITSRecPoint(s) to AliITSclusterV2(s) * -* Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch * -*****************************************************************/ +/**************************************************************************** + * Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch * + ****************************************************************************/ -#ifndef __CINT__ +#if !defined(__CINT__) || defined(__MAKECINT__) #include #include "AliRun.h" #include "AliRunLoader.h" - #include "AliLoader.h" + #include "AliITSLoader.h" #include "AliITS.h" #include "AliITSgeom.h" #include "AliITSclustererV2.h" + #include "TTree.h" #include "TStopwatch.h" #endif -Int_t AliITSFindClustersV2(Char_t SlowOrFast='f') -{ +extern AliRun *gAlice; - cerr<<"AliITSRecPoint(s) -> AliITSclusterV2(s)...\n"; - - if (gAlice) - { +Int_t AliITSFindClustersV2(Int_t nev=5, Char_t SlowOrFast='s') { + + cerr<<"Looking for clusters...\n"; + + if (gAlice) { delete gAlice->GetRunLoader(); delete gAlice; gAlice=0; - } + } - AliRunLoader* rl = AliRunLoader::Open("galice.root"); - if (rl == 0x0) - { - cerr<<"AliITSHits2DigitsDefault.C : Can not open session RL=NULL" - << endl; - return 3; - } + AliRunLoader *rl = AliRunLoader::Open("galice.root"); + if (rl == 0x0) { + cerr<<"AliITSFindClustersV2.C : Can not open session RL=NULL"<< endl; + return 1; + } - Int_t retval = rl->LoadgAlice(); - if (retval) - { - cerr<<"AliITSHits2DigitsDefault.C : LoadgAlice returned error" - << endl; - delete rl; - return 3; - } - gAlice=rl->GetAliRun(); - rl->LoadHeader(); - retval = rl->LoadKinematics(); - if (retval) - { - cerr<<"AliITSHits2DigitsDefault.C : LoadKinematics returned error" - << endl; - delete rl; - return 3; - } - - AliITSLoader* gime = (AliITSLoader*)rl->GetLoader("ITSLoader"); - if (gime == 0x0) - { - cerr<<"AliITSHits2DigitsDefault.C : can not get ITS loader" - << endl; - } - - rl->GetEvent(0); + AliITSLoader *itsl = (AliITSLoader*)rl->GetLoader("ITSLoader"); + if (itsl == 0x0) { + cerr<<"AliITSFindClustersV2.C : can not get ITS loader"<< endl; + return 2; + } + + rl->LoadKinematics(); + + Int_t retval = rl->LoadgAlice(); + if (retval) { + cerr<<"AliITSFindClustersV2.C : LoadgAlice returned error"<< endl; + delete rl; + return 3; + } + gAlice=rl->GetAliRun(); AliITS *ITS = (AliITS*)gAlice->GetModule("ITS"); if (!ITS) { cerr<<"Can't find the ITS !\n"; delete rl; return 3; } AliITSgeom *geom=ITS->GetITSgeom(); - - TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000); - - gime->LoadRawClusters("recreate"); - - if (SlowOrFast=='f') - { - gime->SetRecPointsFileName("ITS.FastRecPoints.root"); - } - if (gime->LoadRecPoints()) - { - cerr<<"Load Rec Pints returned error !\n"; - delete rl; - return 4; - } - - TClonesArray *points = new TClonesArray("AliITSRecPoint",10000); - - Float_t lp[5]; - Int_t lab[6]; - - Int_t iEvent; - for (iEvent = 0; iEvent< rl->GetNumberOfEvents() ; iEvent++) - { - - rl->GetEvent(iEvent); - - TTree *cTree = gime->TreeC(); - if (cTree == 0x0) - { - gime->MakeTree("C"); - cTree = gime->TreeC(); - } - - cTree->Branch("Clusters",&clusters); - - TTree *pTree=gime->TreeR(); - if (pTree == 0x0) - { - cerr<<"Can not get TreeR !\n"; delete rl;return 5; - } - TBranch *branch = 0; - if (SlowOrFast=='f') { - branch = pTree->GetBranch("ITSRecPointsF"); - } - else { - branch = pTree->GetBranch("ITSRecPoints"); - } - if (!branch) - { - cerr<<"Can't get ITSRecPoints branch !\n"; - delete rl; - return 6; - } - - branch->SetAddress(&points); - - AliStack* stack = rl->Stack(); - if (stack == 0x0) - { - cerr<<"AliITSFindClustersV2.C : Can not get stack" - << endl; - delete rl; - return 3; - } - - TClonesArray &cl=*clusters; - Int_t nclusters=0; - Int_t nentr=(Int_t)branch->GetEntries(); - - cerr<<"Number of entries: "<Clear(); - clusters->Clear(); - branch->GetEvent(i); - Int_t ncl=points->GetEntriesFast(); if (ncl==0){cTree->Fill();continue;} - Int_t lay,lad,det; geom->GetModuleId(i,lay,lad,det); - if ( (lay<0) || (lad<0) || (det<0)) - { - ::Error("AliITSFindClustersV2.C","No such a module %d",i); - continue; - } - Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift); - Double_t rot[9]; geom->GetRotMatrix(lay,lad,det,rot); - Double_t yshift = x*rot[0] + y*rot[1]; - Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1); - nclusters+=ncl; - - Float_t kmip=1; // ADC->mip normalization factor for the SDD and SSD - if(lay==4 || lay==3){kmip=280.;}; - if(lay==6 || lay==5){kmip=38.;}; - - for (Int_t j=0; jUncheckedAt(j); - //Float_t lp[5]; - lp[0]=-p->GetX()-yshift; if (lay==1) lp[0]=-lp[0]; - lp[1]=p->GetZ()+zshift; - lp[2]=p->GetSigmaX2(); - lp[3]=p->GetSigmaZ2(); - lp[4]=p->GetQ(); lp[4]/=kmip; - //Int_t lab[6]; - lab[0]=p->GetLabel(0);lab[1]=p->GetLabel(1);lab[2]=p->GetLabel(2); - lab[3]=ndet; - - Int_t label=lab[0]; - if (label>=0) { - TParticle *part=(TParticle*)stack->Particle(label); - if (part == 0x0) - cerr<<"Can not get particle with label "<P() < 0.005) { - Int_t m=part->GetFirstMother(); - if (m<0) {cerr<<"Primary momentum: "<P()<Particle(label); - } - if (lab[1]<0) lab[1]=label; - else if (lab[2]<0) lab[2]=label; - else cerr<<"No empty labels !\n"; - } - - new(cl[j]) AliITSclusterV2(lab,lp); - } - cTree->Fill(); -// clusters->Delete(); points->Delete(); - } - gime->WriteRawClusters("OVERWRITE"); - cerr<<"Number of clusters: "<LoadRecPoints("recreate"); + if (SlowOrFast=='s') itsl->LoadDigits("read"); + else itsl->LoadHits("read"); + + AliITSclustererV2 clusterer(geom); + + TStopwatch timer; + if (nev>rl->GetNumberOfEvents()) nev=rl->GetNumberOfEvents(); + for (Int_t i=0; iGetEvent(i); + + TTree *out=itsl->TreeR(); + if (!out) { + itsl->MakeTree("R"); + out=itsl->TreeR(); + } + + if (SlowOrFast=='s') { + TTree *in=itsl->TreeD(); + if (!in) { + cerr<<"Can't get digits tree !\n"; + return 4; + } + clusterer.Digits2Clusters(in,out); + } else { + TTree *in=itsl->TreeH(); + if (!in) { + cerr<<"Can't get hits tree !\n"; + return 5; + } + clusterer.Hits2Clusters(in,out); + } + + itsl->WriteRecPoints("OVERWRITE"); } - - delete clusters; - delete points; + timer.Stop(); timer.Print(); + delete rl; - return 0; + return 0; } - - - - - - - - - - - diff --git a/ITS/AliITSFindTracksV2.C b/ITS/AliITSFindTracksV2.C index 25723ccbb61..6fc08504ef5 100644 --- a/ITS/AliITSFindTracksV2.C +++ b/ITS/AliITSFindTracksV2.C @@ -1,98 +1,109 @@ -#ifndef __CINT__ +/**************************************************************************** + * Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch * + ****************************************************************************/ + +#if !defined(__CINT__) || defined(__MAKECINT__) #include "Riostream.h" - #include "TFile.h" + #include "TStopwatch.h" + #include "AliRun.h" + #include "AliRunLoader.h" + #include "AliTPCLoader.h" + #include "AliITSLoader.h" + #include "AliITS.h" #include "AliITSgeom.h" #include "AliITStrackerV2.h" #endif -Int_t AliITSFindTracksV2(Int_t nev=1) { //number of events to process +extern AliRun *gAlice; + +Int_t AliITSFindTracksV2(Int_t nev=5) { //number of events to process cerr<<"Looking for tracks...\n"; - if (gAlice) - { + if (gAlice) { delete gAlice->GetRunLoader(); delete gAlice; gAlice=0; - } + } - AliRunLoader* rl = AliRunLoader::Open("galice.root"); - if (rl == 0x0) - { - cerr<<"AliITSHits2DigitsDefault.C : Can not open session RL=NULL" - << endl; - return 3; - } + AliRunLoader* rl = AliRunLoader::Open("galice.root"); + if (rl == 0x0) { + cerr<<"AliITSFindTracks.C : Can not open session RL=NULL"<< endl; + return 3; + } - Int_t retval = rl->LoadgAlice(); - if (retval) - { - ::Error("AliITSHits2DigitsDefault.C","LoadgAlice returned error"); - delete rl; + Int_t retval = rl->LoadgAlice(); + if (retval) { + cerr<<"AliITSFindTracksV2.C : LoadgAlice returned error"<LoadHeader(); - if (retval) - { - ::Error("AliITSHits2DigitsDefault.C","LoadHeader returned error"); + } + retval = rl->LoadHeader(); + if (retval) { + cerr<<"AliITSFindTracksV2.C : LoadHeader returned error"<GetAliRun(); - - - AliITSLoader* itsloader = (AliITSLoader*)rl->GetLoader("ITSLoader"); - if (itsloader == 0x0) - { - ::Error("AliITSHits2DigitsDefault.C","can not get ITS loader"); + } + gAlice=rl->GetAliRun(); + + AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader"); + if (itsl == 0x0) { + cerr<<"AliITSFindTracksV2.C : Can not get ITS loader"<GetLoader("TPCLoader"); - if (tpcloader == 0x0) - { - cerr<<"AliITSHits2DigitsDefault.C : can not get TPC loader" - << endl; - } + } - rl->GetEvent(0); + AliTPCLoader* tpcl = (AliTPCLoader*)rl->GetLoader("TPCLoader"); + if (tpcl == 0x0) { + cerr<<"AliITSFindTracksV2.C : can not get TPC loader"<LoadTracks("recreate"); - tpcloader->LoadTracks("read"); - itsloader->LoadRawClusters("read"); - - AliITS* dITS = (AliITS*)gAlice->GetDetector("ITS"); - if(!dITS) - { - ::Error("AliITSHits2DigitsDefault.C","Can not find ITS detector."); + AliITS *dITS = (AliITS*)gAlice->GetDetector("ITS"); + if (!dITS) { + cerr<<"AliITSFindClusters.C : Can not find the ITS detector !"<GetITSgeom(); - if(geom == 0x0) - { - ::Error("AliITSHits2DigitsDefault.C","Can not get geometry from ITS detector."); - return 6; - } // end if !GetITSgeom() + AliITStrackerV2 tracker(geom); - TStopwatch timer; - - for (Int_t i = 0;i < rl->GetNumberOfEvents(); i++) - { - AliITStrackerV2* tracker = new AliITStrackerV2(geom,i); - Int_t rc=tracker->Clusters2Tracks(); - if (rc) - { - ::Error("AliITSHits2DigitsDefault.C", - "AliITStrackerV2::Clusters2Tracks returned errror for event %d",i); - delete tracker; - break; + tpcl->LoadTracks("read"); + itsl->LoadTracks("recreate"); + itsl->LoadRecPoints("read"); + + TStopwatch timer; + if (nev>rl->GetNumberOfEvents()) nev=rl->GetNumberOfEvents(); + Int_t rc=0; + for (Int_t i=0; iGetEvent(i); + + TTree *cTree=itsl->TreeR(); + if (!cTree) { + cerr<<"AliITSFindTracksV2.C : Can't get the clusters tree !"<TreeT(); + if (!tpcTree) { + cerr<<"AliITSFindTracksV2.C : Can't get the TPC track tree !"<TreeT(); + if (!itsTree) { + itsl->MakeTree("T"); + itsTree=itsl->TreeT(); + } + + tracker.LoadClusters(cTree); + rc=tracker.Clusters2Tracks(tpcTree,itsTree); + tracker.UnloadClusters(); + + itsl->WriteTracks("OVERWRITE"); + } timer.Stop(); timer.Print(); - delete tracker; + delete rl; + return rc; } diff --git a/ITS/AliITSTrackerV1.cxx b/ITS/AliITSTrackerV1.cxx index 98c8e7e320b..b13f01093f7 100644 --- a/ITS/AliITSTrackerV1.cxx +++ b/ITS/AliITSTrackerV1.cxx @@ -510,14 +510,15 @@ void AliITSTrackerV1::DoTracking(Int_t evNumber,Int_t minTr,Int_t maxTr, cf->cd(); TString foldname(fITS->GetLoader()->GetEventFolder()->GetName()); - AliTPCtracker *tracker = new AliTPCtracker(digp,evNumber,foldname); - //PH AliTPCtracker *tracker = new AliTPCtracker(digp); //I.B. + printf("This method is not converted to the NewIO !\n"); //I.B. + return; //I.B. + AliTPCtracker *tracker = new AliTPCtracker(digp); //I.B. //PH tracker->SetEventNumber(evNumber); //I.B. // Load clusters - //tracker->LoadInnerSectors(); //I.B. - //tracker->LoadOuterSectors(); //I.B. - tracker->LoadClusters(); //I.B. + printf("This method is not converted to the NewIO !\n"); //I.B. + return; //I.B. + tracker->LoadClusters(0); //I.B. // Load tracks TFile *tf=TFile::Open("AliTPCtracksSorted.root"); diff --git a/ITS/AliITSclustererV2.cxx b/ITS/AliITSclustererV2.cxx index 75c2d546aab..452c53bd6b1 100644 --- a/ITS/AliITSclustererV2.cxx +++ b/ITS/AliITSclustererV2.cxx @@ -6,10 +6,6 @@ //uncomment the line below for running with V1 cluster finder classes //#define V1 -//PH 19/05/2003 This class hast to be adapted to NewIO - -#include - #include "AliRun.h" #include "AliITSclustererV2.h" @@ -18,7 +14,6 @@ #include "AliITSRawStreamSDD.h" #include "AliITSRawStreamSSD.h" -#include #include #include #include @@ -39,7 +34,10 @@ AliITSclustererV2::AliITSclustererV2(const AliITSgeom *geom) { fI=0; Int_t mmax=geom->GetIndexMax(); - if (mmax>2200) {cerr<<"Too many ITS subdetectors !\n"; exit(1);} + if (mmax>2200) { + Error("AliITSclustererV2","Too many ITS subdetectors !"); + exit(1); + } Int_t m; for (m=0; mGetModuleId(m,lay,lad,det); @@ -89,26 +87,15 @@ AliITSclustererV2::AliITSclustererV2(const AliITSgeom *geom) { fTanN=0.0075; } -void AliITSclustererV2::Digits2Clusters(const TFile *in, TFile *out) { +Int_t AliITSclustererV2::Digits2Clusters(TTree *dTree, TTree *cTree) { //------------------------------------------------------------ // This function creates ITS clusters //------------------------------------------------------------ Int_t ncl=0; - TDirectory *savedir=gDirectory; - if (!out->IsOpen()) { - cerr<<"AliITSclustererV2::Digits2Clusters(): output file not open !\n"; - return; - } - - Char_t name[100]; - sprintf(name,"TreeD%d",fEvent); - - //TTree *dTree=(TTree*)((TFile*)in)->Get(name); - TTree *dTree=gAlice->TreeD(); if (!dTree) { - cerr<<"Input tree "<SetBranchAddress("ITSDigitsSSD",&digitsSSD); - Int_t mmax=(Int_t)dTree->GetEntries(); - - out->cd(); - - sprintf(name,"TreeC_ITS_%d",fEvent); - TTree cTree(name,"ITS clusters V2"); TClonesArray *clusters=new TClonesArray("AliITSclusterV2",1000); - cTree.Branch("Clusters",&clusters); + TBranch *branch=cTree->GetBranch("Clusters"); + if (!branch) cTree->Branch("Clusters",&clusters); + else branch->SetAddress(&clusters); + + Int_t mmax=(Int_t)dTree->GetEntries(); for (fI=0; fIGetEvent(fI); @@ -139,14 +124,15 @@ void AliITSclustererV2::Digits2Clusters(const TFile *in, TFile *out) { ncl+=clusters->GetEntriesFast(); - cTree.Fill(); + cTree->Fill(); digitsSPD->Clear(); digitsSDD->Clear(); digitsSSD->Clear(); clusters->Clear(); } - cTree.Write(); + + //cTree->Write(); delete clusters; @@ -154,11 +140,9 @@ void AliITSclustererV2::Digits2Clusters(const TFile *in, TFile *out) { delete digitsSDD; delete digitsSSD; - //delete dTree; + Info("Digits2Clusters","Number of found clusters : %d",ncl); - cerr<<"Number of found clusters : "<cd(); + return 0; } void AliITSclustererV2::Digits2Clusters(TFile *out) { @@ -167,7 +151,7 @@ void AliITSclustererV2::Digits2Clusters(TFile *out) { //------------------------------------------------------------ TDirectory *savedir=gDirectory; if (!out->IsOpen()) { - cerr<<"AliITSclustererV2::Digits2Clusters(): output file not open !\n"; + Error("Digits2Clusters","Output file not open !"); return; } out->cd(); @@ -229,9 +213,13 @@ static void CheckLabels(Int_t lab[3]) { label=-3; while (part->P() < 0.005) { Int_t m=part->GetFirstMother(); - if (m<0) {cerr<<"Primary momentum: "<P()<P()); + break; + } if (part->GetStatusCode()>0) { - cerr<<"Primary momentum: "<P()<P()); + break; } label=m; part=(TParticle*)gAlice->Particle(label); @@ -266,40 +254,30 @@ void AliITSclustererV2::RecPoints2Clusters } } -void AliITSclustererV2::Hits2Clusters(const TFile *in, TFile *out) { +Int_t AliITSclustererV2::Hits2Clusters(TTree *hTree, TTree *cTree) { //------------------------------------------------------------ // This function creates ITS clusters //------------------------------------------------------------ - TDirectory *savedir=gDirectory; - - if (!out->IsOpen()) { - cerr<<"AliITSclustererV2::Hits2Clusters: output file not open !\n"; - return; - } - if (!gAlice) { - cerr<<"AliITSclustererV2::Hits2Clusters : gAlice==0 !\n"; - return; + Error("Hits2Clusters","gAlice==0 !"); + return 1; } AliITS *its = (AliITS*)gAlice->GetModule("ITS"); if (!its) { - cerr<<"AliITSclustererV2::Hits2Clusters : Can't find the ITS !\n"; - return; + Error("Hits2Clusters","Can't find the ITS !"); + return 2; } AliITSgeom *geom=its->GetITSgeom(); Int_t mmax=geom->GetIndexMax(); its->InitModules(-1,mmax); - its->FillModules(gAlice->TreeH(),0); + its->FillModules(hTree,0); - out->cd(); - - Char_t name[100]; - sprintf(name,"TreeC_ITS_%d",fEvent); - TTree cTree(name,"ITS clusters V2"); TClonesArray *clusters=new TClonesArray("AliITSclusterV2",1000); - cTree.Branch("Clusters",&clusters); + TBranch *branch=cTree->GetBranch("Clusters"); + if (!branch) cTree->Branch("Clusters",&clusters); + else branch->SetAddress(&clusters); static TClonesArray *points=its->RecPoints(); AliITSsimulationFastPoints sim; @@ -312,16 +290,17 @@ void AliITSclustererV2::Hits2Clusters(const TFile *in, TFile *out) { its->ResetRecPoints(); ncl+=clusters->GetEntriesFast(); - cTree.Fill(); + cTree->Fill(); clusters->Clear(); } - cTree.Write(); - cerr<<"Number of found fast clusters : "<Write(); delete clusters; - savedir->cd(); + return 0; } //*********************************** @@ -375,8 +354,10 @@ FindClustersSPD(const TClonesArray *digits, TClonesArray *clusters) { if (!bins[k].IsNotUsed()) continue; Int_t ni=0, idx[200]; FindCluster(k,kNzBins,bins,ni,idx); - if (ni==200) {cerr<<"SPD: Too big cluster !\n"; continue;} - + if (ni==200) { + Info("FindClustersSPD","Too big cluster !"); + continue; + } Int_t lab[4]; lab[0]=-2; lab[1]=-2; @@ -957,8 +938,7 @@ FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) { c[*n].SetQ(0.5*q); (*n)++; if (*n==MAX) { - cerr<< - "AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n"; + Error("FindClustersSSD","Too many 1D clusters !"); return; } c[*n].SetY(y/q+0.5*nd); @@ -968,7 +948,7 @@ FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) { } (*n)++; if (*n==MAX) { - cerr<<"AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n"; + Error("FindClustersSSD","Too many 1D clusters !"); return; } y=q=qmax=0.; @@ -996,7 +976,7 @@ FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) { c[*n].SetQ(0.5*q); (*n)++; if (*n==MAX) { - cerr<<"AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n"; + Error("FindClustersSSD","Too many 1D clusters !"); return; } c[*n].SetY(y/q+0.5*nd); @@ -1006,7 +986,7 @@ FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) { } (*n)++; if (*n==MAX) { - cerr<<"AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n"; + Error("FindClustersSSD","Too many 1D clusters !"); return; } diff --git a/ITS/AliITSclustererV2.h b/ITS/AliITSclustererV2.h index 737cd383a99..177a79c0c95 100644 --- a/ITS/AliITSclustererV2.h +++ b/ITS/AliITSclustererV2.h @@ -9,22 +9,23 @@ // // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch //-------------------------------------------------------------- -#include +#include class TFile; +class TTree; class TClonesArray; class AliITSgeom; class AliITSclusterV2; class AliITSRawStream; -class AliITSclustererV2 { +class AliITSclustererV2 : public TObject { public: AliITSclustererV2(){ fEvent=0; fI=0; } AliITSclustererV2(const AliITSgeom *geom); void SetEvent(Int_t event) { fEvent=event; } - void Digits2Clusters(const TFile *in, TFile *out); + Int_t Digits2Clusters(TTree *in, TTree *out); void Digits2Clusters(TFile *out); void FindClustersSPD(const TClonesArray *dig, TClonesArray *cls); void FindClustersSPD(AliITSRawStream* input, TClonesArray** clusters); @@ -34,7 +35,7 @@ public: void FindClustersSSD(AliITSRawStream* input, TClonesArray** clusters); void RecPoints2Clusters(const TClonesArray *p, Int_t idx, TClonesArray *c); - void Hits2Clusters(const TFile *in, TFile *out); + Int_t Hits2Clusters(TTree *in, TTree *out); private: class Ali1Dcluster { diff --git a/ITS/AliITSpidESD.cxx b/ITS/AliITSpidESD.cxx index d502949afd7..c4109bd1858 100644 --- a/ITS/AliITSpidESD.cxx +++ b/ITS/AliITSpidESD.cxx @@ -68,7 +68,7 @@ Int_t AliITSpidESD::MakePID(AliESD *event) Double_t bethe=Bethe(mom/mass); Double_t sigma=fRes*bethe; if (TMath::Abs(dedx-bethe) > fRange*sigma) { - p[j]=0.; + p[j]=TMath::Exp(-0.5*fRange*fRange); continue; } p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma)); diff --git a/ITS/AliITStestV2.C b/ITS/AliITStestV2.C index a1416f59de9..e2ef8353c8c 100644 --- a/ITS/AliITStestV2.C +++ b/ITS/AliITStestV2.C @@ -1,44 +1,45 @@ -Int_t AliITStestV2(Char_t SlowOrFast='s') { +/**************************************************************************** + * Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch * + ****************************************************************************/ + +Int_t AliITStestV2(Int_t nev=5, Char_t SlowOrFast='s') { Int_t rc=0; - if (gAlice) - { + + if (gAlice) { delete gAlice->GetRunLoader(); delete gAlice; gAlice=0; - } + } + rl = AliRunLoader::Open("galice.root"); - if (rl == 0x0) - { + if (rl == 0x0) { cerr<<"Can not open session"<LoadgAlice()) - { + } + + if (rl->LoadgAlice()) { cerr<<"Error occured while loading AliRun"<GetAliRun()->Field()->SolenoidField()); - delete rl; + } + AliKalmanTrack:: + SetConvConst(1000/0.299792458/rl->GetAliRun()->Field()->SolenoidField()); - - - cerr<<"Fast AliITSRecPoint(s) !\n"; - gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSHits2FastRecPoints.C"); - AliITSHits2FastRecPoints(); + delete rl; - cerr<<"Slow AliITSRecPoint(s) !\n"; - gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSHits2SDigits.C"); - AliITSHits2SDigits(); - gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSSDigits2Digits.C"); - AliITSSDigits2Digits(); - gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSDigits2RecPoints.C"); - AliITSDigits2RecPoints(); - + if (SlowOrFast=='f') { + cerr<<"Fast AliITSRecPoint(s) !\n"; + } else { + cerr<<"Slow AliITSRecPoint(s) !\n"; + gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSHits2SDigits.C"); + AliITSHits2SDigits(); + gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSSDigits2Digits.C"); + AliITSSDigits2Digits(); + } gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSFindClustersV2.C"); - if (rc=AliITSFindClustersV2(SlowOrFast)) return rc; + if (rc=AliITSFindClustersV2(nev,SlowOrFast)) return rc; gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSFindTracksV2.C"); - if (rc=AliITSFindTracksV2()) return rc; + if (rc=AliITSFindTracksV2(nev)) return rc; gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSComparisonV2.C"); if (rc=AliITSComparisonV2()) return rc; diff --git a/ITS/AliITStrackerV2.cxx b/ITS/AliITStrackerV2.cxx index c7301371b28..2e5bb19b9e4 100644 --- a/ITS/AliITStrackerV2.cxx +++ b/ITS/AliITStrackerV2.cxx @@ -27,20 +27,15 @@ #include "AliITSgeom.h" #include "AliITSRecPoint.h" #include "AliTPCtrack.h" +#include "AliESD.h" #include "AliITSclusterV2.h" #include "AliITStrackerV2.h" -#include "AliRunLoader.h" -#include "AliLoader.h" -#include "AliITSLoader.h" - ClassImp(AliITStrackerV2) AliITStrackerV2::AliITSlayer AliITStrackerV2::fLayers[kMaxLayer]; // ITS layers -AliITStrackerV2::AliITStrackerV2(const AliITSgeom *geom, Int_t eventn, const char* evfoldname): - AliTracker(), - fEvFolderName(evfoldname) { +AliITStrackerV2::AliITStrackerV2(const AliITSgeom *geom) : AliTracker() { //-------------------------------------------------------------------- //This is the AliITStrackerV2 constructor //-------------------------------------------------------------------- @@ -48,8 +43,7 @@ AliITStrackerV2::AliITStrackerV2(const AliITSgeom *geom, Int_t eventn, const ch Float_t x,y,z; Int_t i; - for (i=1; iGetNladders(i); Int_t ndet=g->GetNdetectors(i); @@ -81,11 +75,11 @@ AliITStrackerV2::AliITStrackerV2(const AliITSgeom *geom, Int_t eventn, const ch new(&det) AliITSdetector(r,phi); } } + } fI=kMaxLayer; - fPass=0; fConstraint[0]=1; fConstraint[1]=0; @@ -104,43 +98,13 @@ void AliITStrackerV2::SetLayersNotToSkip(Int_t *l) { for (Int_t i=0; iGetLoader("ITSLoader"); - if (itsl == 0x0) - { - Error("LoadClusters","Can not get ITS loader."); - return 1; - } - - Int_t retval; - TTree *cTree=itsl->TreeC(); - if (!cTree) { - retval = itsl->LoadRawClusters("read"); - if (retval) - { - Error("LoadClusters","LoadRawClusters returned error."); - return 1; - } - cTree=itsl->TreeC(); - } - - if (!cTree) { - Error("LoadClusters"," can't get cTree !\n"); - return 1; - } TBranch *branch=cTree->GetBranch("Clusters"); if (!branch) { - Error("LoadClusters"," can't get Clusters branch !\n"); + Error("LoadClusters"," can't get the branch !\n"); return 1; } @@ -162,8 +126,7 @@ Int_t AliITStrackerV2::LoadClusters() { } fLayers[i].ResetRoad(); //road defined by the cluster density } - - itsl->UnloadRawClusters(); + return 0; } @@ -208,44 +171,101 @@ static Int_t CorrectForDeadZoneMaterial(AliITStrackV2 *t) { return 0; } -Int_t AliITStrackerV2::Clusters2Tracks() { +Int_t AliITStrackerV2::Clusters2Tracks(AliESD *event) { //-------------------------------------------------------------------- - //This functions reconstructs ITS tracks + // This functions reconstructs ITS tracks + // The clusters must be already loaded ! //-------------------------------------------------------------------- + TObjArray itsTracks(15000); - - AliRunLoader* rl = AliRunLoader::GetRunLoader(fEvFolderName); - if (rl == 0x0) - { - Error("Clusters2Tracks","Can not get RL from specified folder %s",fEvFolderName.Data()); - return 1; - } - rl->GetEvent(GetEventNumber()); - //leave loading clusters here - than it is not necessary to GetEvent two times - if (LoadClusters()!=0) return 1; - - AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader"); - AliLoader* tpcl = rl->GetLoader("TPCLoader"); - if ( !tpcl || !itsl) - { - Error("Clusters2Tracks","Can not get loaders"); - return 2; - } - + {/* Read ESD tracks */ + Int_t nentr=event->GetNumberOfTracks(); + Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr); + while (nentr--) { + AliESDtrack *esd=event->GetTrack(nentr); + + if (esd->GetStatus() != AliESDtrack::kTPCin) continue; + + AliITStrackV2 *t=0; + try { + t=new AliITStrackV2(*esd); + } catch (const Char_t *msg) { + Warning("Clusters2Tracks",msg); + delete t; + continue; + } + if (TMath::Abs(t->GetD())>4) continue; + if (CorrectForDeadZoneMaterial(t)!=0) { + Warning("Clusters2Tracks", + "failed to correct for the material in the dead zone !\n"); + delete t; + continue; + } + itsTracks.AddLast(t); + } + } /* End Read ESD tracks */ + + itsTracks.Sort(); + Int_t nentr=itsTracks.GetEntriesFast(); + + Int_t ntrk=0; + for (fPass=0; fPass<2; fPass++) { + Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue; + for (Int_t i=0; iGetLabel(); //save the TPC track label + + ResetTrackToFollow(*t); + ResetBestTrack(); + + for (FollowProlongation(); fITreeT() == 0x0) tpcl->LoadTracks("read"); - TTree *tpcTree = tpcl->TreeT(); - if (!tpcTree) { - Error("Clusters2Tracks","can't get a tree with TPC tracks !\n"); - return 3; - } AliTPCtrack *itrack=new AliTPCtrack; + TBranch *branch=tpcTree->GetBranch("tracks"); + if (!branch) { + Error("Clusters2Tracks","Can't get the branch !"); + return 1; + } tpcTree->SetBranchAddress("tracks",&itrack); nentr=(Int_t)tpcTree->GetEntries(); + Info("Clusters2Tracks","Number of TPC tracks: %d\n",nentr); for (Int_t i=0; iTreeT() == 0x0) itsl->MakeTree("T"); - - TTree& itsTree = *itsl->TreeT(); AliITStrackV2 *otrack=&fBestTrack; - - itsTree.Branch("tracks","AliITStrackV2",&otrack,32000,0); + TBranch *branch=itsTree->GetBranch("tracks"); + if (!branch) itsTree->Branch("tracks","AliITStrackV2",&otrack,32000,3); + else branch->SetAddress(&otrack); for (fPass=0; fPass<2; fPass++) { Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue; @@ -291,65 +309,118 @@ Int_t AliITStrackerV2::Clusters2Tracks() { ResetTrackToFollow(*t); ResetBestTrack(); - for (FollowProlongation(); fIFill(); UseClusters(&fBestTrack); delete itsTracks.RemoveAt(i); } } - nentr=(Int_t)itsTree.GetEntries(); + nentr=(Int_t)itsTree->GetEntries(); Info("Clusters2Tracks","Number of prolonged tracks: %d\n",nentr); - itsl->WriteTracks("OVERWRITE"); - itsTracks.Delete(); - UnloadClusters(); + return 0; +} + +Int_t AliITStrackerV2::PropagateBack(AliESD *event) { + //-------------------------------------------------------------------- + // This functions propagates reconstructed ITS tracks back + // The clusters must be loaded ! + //-------------------------------------------------------------------- + Int_t nentr=event->GetNumberOfTracks(); + Info("PropagateBack", "Number of ESD tracks: %d\n", nentr); + + Int_t ntrk=0; + for (Int_t i=0; iGetTrack(i); + + if (esd->GetStatus()!=(AliESDtrack::kTPCin|AliESDtrack::kITSin)) continue; + + AliITStrackV2 *t=0; + try { + t=new AliITStrackV2(*esd); + } catch (const Char_t *msg) { + Warning("PropagateBack",msg); + delete t; + continue; + } + + ResetTrackToFollow(*t); + + // propagete to vertex [SR, GSI 17.02.2003] + fTrackToFollow.PropagateTo(3.,0.0028,65.19); + fTrackToFollow.PropagateToVertex(); + + // Start Time measurement [SR, GSI 17.02.2003] + fTrackToFollow.StartTimeIntegral(); + fTrackToFollow.PropagateTo(3.,-0.0028,65.19); + + fTrackToFollow.ResetCovariance(); fTrackToFollow.ResetClusters(); + if (RefitAt(49.,&fTrackToFollow,t)) { + if (CorrectForDeadZoneMaterial(&fTrackToFollow)!=0) { + Warning("PropagateBack", + "failed to correct for the material in the dead zone !\n"); + delete t; + continue; + } + fTrackToFollow.SetLabel(t->GetLabel()); + fTrackToFollow.CookdEdx(); + CookLabel(&fTrackToFollow,0.); //For comparison only + fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout); + UseClusters(&fTrackToFollow); + ntrk++; + } + delete t; + } + + Info("PropagateBack","Number of back propagated ITS tracks: %d\n",ntrk); return 0; } -Int_t AliITStrackerV2::PropagateBack() { +Int_t AliITStrackerV2::PropagateBack(TTree *inp, TTree *out) { //-------------------------------------------------------------------- //This functions propagates reconstructed ITS tracks back //-------------------------------------------------------------------- + Error("PropagateBack","This method is not converted to NewIO yet\n"); + return 1; + /* + TFile *in=(TFile*)inp; + TDirectory *savedir=gDirectory; + if (LoadClusters()!=0) return 1; - AliRunLoader* rl = AliRunLoader::GetRunLoader(fEvFolderName); - if (rl == 0x0) - { - Error("Clusters2Tracks","Can not get RL from specified folder %s",fEvFolderName.Data()); - return 1; - } - rl->GetEvent(GetEventNumber()); - AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader"); - if (itsl == 0x0) - { - Error("LoadClusters","Can not get ITS loader."); + if (!in->IsOpen()) { + Error("PropagateBack","file with ITS tracks is not open !\n"); return 1; - } - - if (itsl->TreeT() == 0x0) itsl->LoadTracks("read"); - TTree *itsTree=itsl->TreeT(); + } + + if (!out->IsOpen()) { + Error("PropagateBack","file for back propagated ITS tracks is not open !\n"); + return 2; + } + + in->cd(); + + Char_t tname[100]; + sprintf(tname,"TreeT_ITS_%d",GetEventNumber()); + TTree *itsTree=(TTree*)in->Get(tname); if (!itsTree) { Error("PropagateBack","can't get a tree with ITS tracks !\n"); return 3; @@ -357,15 +428,10 @@ Int_t AliITStrackerV2::PropagateBack() { AliITStrackV2 *itrack=new AliITStrackV2; itsTree->SetBranchAddress("tracks",&itrack); - - itsl->MakeTree("B");//nake tree for back propagated tracks - if (itsl->TreeB() == 0x0) - { - Error("PropagateBack","Can not create tree for back propagated tracks !\n"); - return 3; - } - TTree &backTree = *itsl->TreeB(); + out->cd(); + sprintf(tname,"TreeT_ITSb_%d",GetEventNumber()); + TTree backTree(tname,"Tree with back propagated ITS tracks"); AliTPCtrack *otrack=0; backTree.Branch("tracks","AliTPCtrack",&otrack,32000,2); @@ -400,42 +466,99 @@ Int_t AliITStrackerV2::PropagateBack() { UseClusters(&fTrackToFollow); } i=(Int_t)backTree.GetEntries(); - - - itsl->WriteBackTracks("OVERWRITE"); + backTree.Write(); Info("PropagateBack","Number of ITS tracks: %d\n",nentr); Info("PropagateBack","Number of back propagated ITS tracks: %d\n",i); delete itrack; - - itsl->UnloadTracks(); - itsl->UnloadBackTracks(); - + delete itsTree; //Thanks to Mariana Bondila + UnloadClusters(); + savedir->cd(); + return 0; + */ } -Int_t AliITStrackerV2::RefitInward() { +Int_t AliITStrackerV2::RefitInward(AliESD *event) { //-------------------------------------------------------------------- // This functions refits ITS tracks using the // "inward propagated" TPC tracks + // The clusters must be loaded ! //-------------------------------------------------------------------- + Int_t nentr=event->GetNumberOfTracks(); + Info("RefitInward", "Number of ESD tracks: %d\n", nentr); + + Int_t ntrk=0; + for (Int_t i=0; iGetTrack(i); + + ULong_t flags = AliESDtrack::kITSin | AliESDtrack::kTPCrefit; + + if ( (esd->GetStatus() & flags) != flags ) continue; + if ( esd->GetStatus() & AliESDtrack::kITSrefit) continue; + + AliITStrackV2 *t=0; + try { + t=new AliITStrackV2(*esd); + } catch (const Char_t *msg) { + Warning("RefitInward",msg); + delete t; + continue; + } + + if (CorrectForDeadZoneMaterial(t)!=0) { + Warning("RefitInward", + "failed to correct for the material in the dead zone !\n"); + delete t; + continue; + } + + ResetTrackToFollow(*t); + fTrackToFollow.ResetClusters(); + + //Refitting... + if (RefitAt(3.7, &fTrackToFollow, t)) { + fTrackToFollow.SetLabel(t->GetLabel()); + fTrackToFollow.CookdEdx(); + CookLabel(&fTrackToFollow,0.); //For comparison only + fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit); + UseClusters(&fTrackToFollow); + ntrk++; + } + delete t; + } + + Info("RefitInward","Number of refitted tracks: %d\n",ntrk); - //PH 19/05 This function has to be adapted to the NewIO - TFile * inp = 0x0; - TFile * out = 0x0; + return 0; +} +Int_t AliITStrackerV2::RefitInward(TTree *in, TTree *out) { + //-------------------------------------------------------------------- + // This functions refits ITS tracks using the + // "inward propagated" TPC tracks + //-------------------------------------------------------------------- + Error("RefitInward","This method is not converted to NewIO yet\n"); + return 1; + /* TFile *in=(TFile*)inp; + TDirectory *savedir=gDirectory; if (LoadClusters()!=0) return 1; - if (!in->IsOpen()) { + if (!inSeeds->IsOpen()) { Error("RefitInward","file with inward TPC tracks is not open !\n"); return 2; } + if (!in->IsOpen()) { + Error("RefitInward","file with ITS tracks is not open !\n"); + return 2; + } + if (!out->IsOpen()) { Error("RefitInward","file for inward ITS tracks is not open !\n"); return 3; @@ -451,9 +574,9 @@ Int_t AliITStrackerV2::RefitInward() { Char_t tname[100]; TObjArray itsTracks(15000); - {/* Read the ITS tracks */ + {// Read the ITS tracks sprintf(tname,"TreeT_ITS_%d",GetEventNumber()); - TTree *itsTree=(TTree*)out->Get(tname); + TTree *itsTree=(TTree*)in->Get(tname); if (!itsTree) { Error("RefitInward","can't get a tree with ITS tracks !\n"); return 3; @@ -480,17 +603,10 @@ Int_t AliITStrackerV2::RefitInward() { delete itrack; } - out->cd(); - - //Create the output tree - sprintf(tname,"TreeT_ITSinward_%d",GetEventNumber()); - TTree outTree(tname,"Tree with inward refitted ITS tracks"); - AliITStrackV2 *otrack=0; - outTree.Branch("tracks","AliITStrackV2",&otrack,32000,0); - - //Get the input tree + //Get the input seeds tree + inSeeds->cd(); sprintf(tname,"tracksTPC_%d",GetEventNumber()); - TTree *tpcTree=(TTree*)in->Get(tname); + TTree *tpcTree=(TTree*)inSeeds->Get(tname); if (!tpcTree) { Error("RefitInward","can't get a tree with TPC tracks !\n"); return 3; @@ -501,18 +617,24 @@ Int_t AliITStrackerV2::RefitInward() { Info("RefitInward","Number of TPC tracks: %d\n",ntpc); + //Create the output tree + out->cd(); + sprintf(tname,"TreeT_ITSinward_%d",GetEventNumber()); + TTree outTree(tname,"Tree with inward refitted ITS tracks"); + AliITStrackV2 *otrack=0; + outTree.Branch("tracks","AliITStrackV2",&otrack,32000,0); + for (i=0; iGetEvent(i); - AliITStrackV2 *t=0; try { - t=new AliITStrackV2(*itrack); + otrack=new AliITStrackV2(*itrack); } catch (const Char_t *msg) { Warning("RefitInward",msg); - delete t; + delete otrack; continue; } //check if this track was reconstructed in the ITS - Int_t lab=TMath::Abs(t->GetLabel()); + Int_t lab=TMath::Abs(otrack->GetLabel()); if (lab >= nLab) { Warning("RefitInward","Too big TPC track label: %d\n!",lab); continue; @@ -520,20 +642,22 @@ Int_t AliITStrackerV2::RefitInward() { Int_t idx=lut[lab]; if (idx<0) continue; //no prolongation in the ITS for this track - if (CorrectForDeadZoneMaterial(t)!=0) { + if (CorrectForDeadZoneMaterial(otrack)!=0) { Warning("RefitInward", "failed to correct for the material in the dead zone !\n"); continue; } //Refitting... - otrack=(AliITStrackV2*)itsTracks.UncheckedAt(idx); - if (!RefitAt(3.7, t, otrack)) continue; + { + AliITStrackV2 *ctrack=(AliITStrackV2*)itsTracks.UncheckedAt(idx); + if (!RefitAt(3.7, otrack, ctrack)) continue; + } otrack->SetLabel(itrack->GetLabel()); //For comparison only otrack->CookdEdx(); CookLabel(otrack,0.); //For comparison only outTree.Fill(); - delete t; + delete otrack; } i=(Int_t)outTree.GetEntries(); Info("RefitInward","Number of inward refitted ITS tracks: %d\n",i); @@ -541,10 +665,13 @@ Int_t AliITStrackerV2::RefitInward() { delete tpcTree; delete itrack; - + UnloadClusters(); itsTracks.Delete(); + savedir->cd(); + return 0; + */ } AliCluster *AliITStrackerV2::GetCluster(Int_t index) const { @@ -861,10 +988,6 @@ FindDetectorIndex(Double_t phi, Double_t z) const { if (nz>=fNdetectors) return -1; if (nz<0) return -1; -#ifdef DEBUG -cout<GetNumberOfClusters(); + Int_t nc=c->GetNumberOfClusters(); for (k=0; kGetClusterIndex(k),nl=(idx&0xf0000000)>>28; + Int_t idx=c->GetClusterIndex(k),nl=(idx&0xf0000000)>>28; index[nl]=idx; } - t->~AliITStrackV2(); new (t) AliITStrackV2(*s); Int_t from, to, step; if (xx > t->GetX()) { @@ -1098,8 +1220,12 @@ AliITStrackerV2::RefitAt(Double_t xx,const AliITStrackV2 *s,AliITStrackV2 *ot) { t->SetDetectorIndex(idet); } Double_t chi2=t->GetPredictedChi2(c); - if (chi2PropagateTo(xx,0.,0.)) return kFALSE; - ot->~AliITStrackV2(); new (ot) AliITStrackV2(*t); return kTRUE; } diff --git a/ITS/AliITStrackerV2.h b/ITS/AliITStrackerV2.h index b1fb48fcd2e..b5983255ed4 100644 --- a/ITS/AliITStrackerV2.h +++ b/ITS/AliITStrackerV2.h @@ -11,31 +11,30 @@ #include "AliTracker.h" #include "AliITSrecoV2.h" #include "AliITStrackV2.h" -#include "AliConfig.h" class AliITSclusterV2; +class AliESD; class AliITSgeom; -class TFile; +class TTree; //------------------------------------------------------------------------- class AliITStrackerV2 : public AliTracker { public: AliITStrackerV2():AliTracker(){} - - AliITStrackerV2(const AliITSgeom *geom, Int_t event=0, - const char* evfoldname = AliConfig::fgkDefaultEventFolderName); - + AliITStrackerV2(const AliITSgeom *geom); AliCluster *GetCluster(Int_t index) const; - Int_t LoadClusters(); + Int_t LoadClusters(TTree *cf); void UnloadClusters(); - Int_t Clusters2Tracks();// - Int_t PropagateBack();// - Int_t RefitInward();// - + Int_t Clusters2Tracks(TTree *in, TTree *out); + Int_t Clusters2Tracks(AliESD *event); + Int_t PropagateBack(TTree *in, TTree *out); + Int_t PropagateBack(AliESD *event); + Int_t RefitInward(TTree *in, TTree *out); + Int_t RefitInward(AliESD *event); + Bool_t RefitAt(Double_t x, AliITStrackV2 *seed, const AliITStrackV2 *t); void SetupFirstPass(Int_t *flags, Double_t *cuts=0); void SetupSecondPass(Int_t *flags, Double_t *cuts=0); - //PH Bool_t RefitAt(Double_t xx, AliITStrackV2 *t, Int_t *index); void SetLastLayerToTrackTo(Int_t l=0) {fLastLayerToTrackTo=l;} void SetLayersNotToSkip(Int_t *l); @@ -97,7 +96,6 @@ private: Double_t GetEffectiveThickness(Double_t y, Double_t z) const; void FollowProlongation(); Int_t TakeNextProlongation(); - Bool_t RefitAt(Double_t x, const AliITStrackV2 *t, AliITStrackV2 *tt); void ResetBestTrack() { fBestTrack.~AliITStrackV2(); new(&fBestTrack) AliITStrackV2(fTrackToFollow); @@ -113,12 +111,11 @@ private: AliITStrackV2 fTrackToFollow; // followed track Int_t fPass; // current pass through the data Int_t fConstraint[2]; // constraint flags - TString fEvFolderName; //event folder name Int_t fLayersNotToSkip[kMaxLayer]; // layer masks Int_t fLastLayerToTrackTo; // the innermost layer to track to - ClassDef(AliITStrackerV2,2) //ITS tracker V2 + ClassDef(AliITStrackerV2,1) //ITS tracker V2 }; diff --git a/ITS/ITSLinkDef.h b/ITS/ITSLinkDef.h index 10c95411d9b..46f8fc5980b 100644 --- a/ITS/ITSLinkDef.h +++ b/ITS/ITSLinkDef.h @@ -151,7 +151,7 @@ // This class will always be for ITS only #pragma link C++ class AliITSvtest+; -//PH #pragma link C++ class AliITSclustererV2+; +#pragma link C++ class AliITSclustererV2+; #pragma link C++ class AliITSclusterV2+; #pragma link C++ class AliITStrackV2+; #pragma link C++ class AliITStrackerV2+; diff --git a/ITS/libITS.pkg b/ITS/libITS.pkg index 08618499f99..857c3effc6f 100644 --- a/ITS/libITS.pkg +++ b/ITS/libITS.pkg @@ -31,6 +31,7 @@ SRCS = AliITS.cxx AliITSv1.cxx AliITSv5.cxx AliITSvSPD02.cxx \ AliITSTrackV1.cxx AliITSIOTrack.cxx \ AliITSRad.cxx AliITSgeoinfo.cxx AliITSTrackerV1.cxx\ AliITSvtest.cxx AliITSStrLine.cxx\ + AliITSclustererV2.cxx \ AliITSclusterV2.cxx AliITStrackV2.cxx AliITStrackerV2.cxx \ AliITSVertex.cxx AliITSVertexer.cxx AliITSVertexerIons.cxx \ AliITSVertexerPPZ.cxx AliITSVertexerTracks.cxx\ diff --git a/STEER/AliESDanalysis.C b/STEER/AliESDanalysis.C index 7593cb69587..7355d2a0ad9 100644 --- a/STEER/AliESDanalysis.C +++ b/STEER/AliESDanalysis.C @@ -5,45 +5,75 @@ // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch //******************************************************************** -#ifndef __CINT__ +#if !defined( __CINT__) || defined(__MAKECINT__) #include #include "TKey.h" #include "TFile.h" + #include "TH1F.h" #include "TH2F.h" #include "TCanvas.h" #include "TStopwatch.h" + #include "TParticle.h" + + #include "AliRun.h" + #include "AliRunLoader.h" + #include "AliLoader.h" #include "AliESD.h" #endif +extern AliRun *gAlice; + Int_t AliESDanalysis(Int_t nev=1) { TH2F *tpcHist= - new TH2F("tpcHist","TPC dE/dX vs momentum",100,0.,4.,100,0.,500.); - tpcHist->SetXTitle("p (GeV/c)"); tpcHist->SetYTitle("dE/dx (Arb. Units)"); + new TH2F("tpcHist","TPC dE/dX vs momentum",100,0.,3.,100,0.,500.); + tpcHist->SetXTitle("p (GeV/c)"); tpcHist->SetYTitle("dE/dx (Arb. Units)"); tpcHist->SetMarkerStyle(8); tpcHist->SetMarkerSize(0.3); - TH2F *elHist=new TH2F("elHist","dE/dX vs momentum",100,0.,4.,50,0.,500.); - elHist->SetMarkerStyle(8); - elHist->SetMarkerSize(0.3); - - TH2F *piHist=new TH2F("piHist","dE/dX vs momentum",100,0.,4.,50,0.,500.); - piHist->SetMarkerColor(2); - piHist->SetMarkerStyle(8); - piHist->SetMarkerSize(0.3); - - TH2F *kaHist=new TH2F("kaHist","dE/dX vs momentum",100,0.,4.,100,0.,500.); - kaHist->SetMarkerColor(4); - kaHist->SetMarkerStyle(8); - kaHist->SetMarkerSize(0.3); - - TH2F *prHist=new - TH2F("prHist","Classification into e, pi, K and p",100,0.,4.,100,0.,500.); - prHist->SetXTitle("p (GeV/c)"); prHist->SetYTitle("dE/dx (Arb. Units)"); - prHist->SetMarkerColor(6); - prHist->SetMarkerStyle(8); - prHist->SetMarkerSize(0.3); - + TH1F *piG=new TH1F("piG","",20,0.,3.); + TH1F *piR=new TH1F("piR","",20,0.,3.); + TH1F *piF=new TH1F("piF","",20,0.,3.); + TH1F *piGood=new TH1F("piGood","Combined PID for pions",20,0.,3.); + piGood->SetLineColor(4); piGood->SetXTitle("p (GeV/c)"); + TH1F *piFake=new TH1F("piFake","",20,0.,3.); piFake->SetLineColor(2); + + TH1F *kaG=new TH1F("kaG","",20,0.,3.); + TH1F *kaR=new TH1F("kaR","",20,0.,3.); + TH1F *kaF=new TH1F("kaF","",20,0.,3.); + TH1F *kaGood=new TH1F("kaGood","Combined PID for kaons",20,0.,3.); + kaGood->SetLineColor(4); kaGood->SetXTitle("p (GeV/c)"); + TH1F *kaFake=new TH1F("kaFake","",20,0.,3.); kaFake->SetLineColor(2); + + TH1F *prG=new TH1F("prG","",20,0.,3.); + TH1F *prR=new TH1F("prR","",20,0.,3.); + TH1F *prF=new TH1F("prF","",20,0.,3.); + TH1F *prGood=new TH1F("prGood","Combined PID for protons",20,0.,3.); + prGood->SetLineColor(4); prGood->SetXTitle("p (GeV/c)"); + TH1F *prFake=new TH1F("prFake","",20,0.,3.); prFake->SetLineColor(2); + + if (gAlice) { + delete gAlice->GetRunLoader(); + delete gAlice; + gAlice=0; + } + AliRunLoader *rl = AliRunLoader::Open("galice.root"); + if (rl == 0x0) { + cerr<<"Can not open session"<LoadgAlice()) { + cerr<<"LoadgAlice returned error"<LoadHeader()) { + cerr<<"LoadHeader returned error"<LoadKinematics(); + TFile *ef=TFile::Open("AliESDs.root"); if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;} @@ -53,16 +83,22 @@ Int_t AliESDanalysis(Int_t nev=1) { TIter next(ef->GetListOfKeys()); //****** Tentative particle type "concentrations" - Double_t c[5]={0.05, 0., 0.85, 0.10, 0.05}; + Double_t c[5]={0.01, 0.01, 0.85, 0.10, 0.05}; //******* The loop over events while ((key=(TKey*)next())!=0) { + rl->GetEvent(n); + cerr<<"Processing event number : "<ReadObj(); Int_t ntrk=event->GetNumberOfTracks(); cerr<<"Number of ESD tracks : "<GetTrack(ntrk); @@ -73,41 +109,91 @@ Int_t AliESDanalysis(Int_t nev=1) { tpcHist->Fill(p,dedx,1); } - if (t->GetStatus()&AliESDtrack::kESDpid) { - Double_t dedx=t->GetTPCsignal(); - Double_t r[10]; t->GetESDpid(r); - Double_t rc=0.; + UInt_t status=AliESDtrack::kESDpid; + status|=AliESDtrack::kTPCpid; + status|=AliESDtrack::kTOFpid; + + if ((t->GetStatus()&status) == status) { + nsel++; + + Int_t lab=TMath::Abs(t->GetLabel()); + TParticle *part=gAlice->Particle(lab); + Int_t code=part->GetPdgCode(); + + Double_t r[10]; t->GetESDpid(r); + + Double_t rcc=0.; Int_t i; - for (i=0; iw[3] && w[4]>w[2] && w[4]>w[1] && w[4]>w[0]) {//proton + prsel++; + prG->Fill(p); + if (TMath::Abs(code)==2212) prR->Fill(p); + else prF->Fill(p); + } + + if (w[3]>w[4] && w[3]>w[2] && w[3]>w[1] && w[3]>w[0]) {//kaon + kasel++; + kaG->Fill(p); + if (TMath::Abs(code)==321) kaR->Fill(p); + else kaF->Fill(p); + } + + if (w[2]>w[3] && w[2]>w[4] && w[2]>w[1] && w[2]>w[0]) {//pion + pisel++; + piG->Fill(p); + if (TMath::Abs(code)==211) piR->Fill(p); + else piF->Fill(p); + } - if (w[4]>w[3] && w[4]>w[2] && w[4]>w[0]) prHist->Fill(p,dedx,1); - if (w[3]>w[4] && w[3]>w[2] && w[3]>w[0]) kaHist->Fill(p,dedx,1); - if (w[2]>w[3] && w[2]>w[4] && w[2]>w[0]) piHist->Fill(p,dedx,1); - if (w[0]>w[3] && w[0]>w[2] && w[0]>w[4]) elHist->Fill(p,dedx,1); } } delete event; + cerr<<"Number of selected ESD tracks : "<Divide(1,2); + c1->Divide(1,4); c1->cd(1); tpcHist->Draw(); + c1->cd(2); - prHist->Draw(); - kaHist->Draw("same"); - piHist->Draw("same"); - elHist->Draw("same"); + piG->Sumw2(); piF->Sumw2(); piR->Sumw2(); + piGood->Divide(piR,piG,1,1,"b"); + piFake->Divide(piF,piG,1,1,"b"); + piGood->Draw("hist"); + piFake->Draw("same"); + + c1->cd(3); + kaG->Sumw2(); kaF->Sumw2(); kaR->Sumw2(); + kaGood->Divide(kaR,kaG,1,1,"b"); + kaFake->Divide(kaF,kaG,1,1,"b"); + kaGood->Draw("hist"); + kaFake->Draw("same"); + + c1->cd(4); + prG->Sumw2(); prF->Sumw2(); prR->Sumw2(); + prGood->Divide(prR,prG,1,1,"b"); + prFake->Divide(prF,prG,1,1,"b"); + prGood->Draw("hist"); + prFake->Draw("same"); ef->Close(); + delete rl; + return rc; } diff --git a/STEER/AliESDpid.cxx b/STEER/AliESDpid.cxx index 1b92e349a04..4fac2ff281d 100644 --- a/STEER/AliESDpid.cxx +++ b/STEER/AliESDpid.cxx @@ -47,6 +47,18 @@ Int_t AliESDpid::MakePID(AliESD *event) for (Int_t j=0; jGetStatus()&AliESDtrack::kTRDpid )!=0) { + Double_t d[10]; + t->GetTRDpid(d); + for (Int_t j=0; jGetStatus()&AliESDtrack::kTOFpid )!=0) { + Double_t d[10]; + t->GetTOFpid(d); + for (Int_t j=0; jSetESDpid(p); } return 0; diff --git a/STEER/AliESDtest.C b/STEER/AliESDtest.C index fc0bf081bc0..ffaa41e209f 100644 --- a/STEER/AliESDtest.C +++ b/STEER/AliESDtest.C @@ -1,125 +1,253 @@ //******************************************************************** // Example of the reconstruction that generates the ESD // Input files: -// a) AliTPCclusters.root containing the TPC clusters -// (the AliTPCFindClusters.C macro can be used to generate it) -// b) AliITSclustersV2.root containing the ITS clusters +// a) file containing the ITS clusters // (the AliITSFindClustersV2.C macro can be used to generate it) +// b) file containing the TPC clusters +// (the AliTPCFindClusters.C macro can be used to generate it) +// c) file containing the TRD clusters +// (the AliTRDdigits2cluster.C macro can be used to generate it) +// d) file containing the TOF digits +// (the AliTOFSDigits2Digits.C macro can be used to generate it) // Ouput file: // AliESDs.root containing the ESD events // // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch //******************************************************************** -#ifndef __CINT__ +#if !defined(__CINT__) || defined(__MAKECINT__) #include #include "TFile.h" + #include "TSystem.h" #include "TStopwatch.h" + #include "TGeant3.h" + + #include "AliMagF.h" + #include "AliRun.h" + #include "AliRunLoader.h" + #include "AliLoader.h" #include "AliESD.h" #include "AliESDpid.h" - #include "AliTPCpidESD.h" - #include "AliTPCParam.h" - #include "AliTPCtracker.h" + + #include "AliITS.h" #include "AliITSgeom.h" #include "AliITStrackerV2.h" + #include "AliITSpidESD.h" + #include "AliITSLoader.h" + + #include "AliTPCParam.h" + #include "AliTPCtracker.h" + #include "AliTPCpidESD.h" + #include "AliTPCLoader.h" + #include "AliTRDtracker.h" #include "AliTRDPartID.h" - #include "AliITSpidESD.h" + + #include "AliTOFpidESD.h" #endif -Int_t AliESDtest(Int_t nev=1, - const char* fileNameITSClusters = "its.clusters.root", - const char* fileNameTPCClusters = "tpc.clusters.root", - const char* fileNameTRDClusters = "trd.clusters.root") { +extern TSystem *gSystem; +extern AliRun *gAlice; +extern TFile *gFile; + +Int_t AliESDtest(Int_t nev=1) { + +/**** Initialization of the NewIO *******/ - //File with the TPC clusters - TFile *tpccf=TFile::Open(fileNameTPCClusters); - if (!tpccf->IsOpen()) { - cerr<<"Can't open "<GetRunLoader(); + delete gAlice; + gAlice=0; + } + + gSystem->Load("libgeant321"); // needed for the PID in TOF + new TGeant3(""); // must be re-done ! + + AliRunLoader *rl = AliRunLoader::Open("galice.root"); + if (rl == 0x0) { + cerr<<"Can not open session"<LoadgAlice(); + if (retval) { + cerr<<"AliESDtest.C : LoadgAlice returned error"<LoadHeader(); + if (retval) { + cerr<<"AliESDtest.C : LoadHeader returned error"<Get("75x40_100x60_150x60"); - if (!par) {cerr<<"Can't get TPC parameters !\n"; return 3;} + gAlice=rl->GetAliRun(); + - //An instance of the TPC tracker - AliTPCtracker tpcTracker(par); + AliKalmanTrack::SetConvConst( + 1000/0.299792458/gAlice->Field()->SolenoidField() + ); - //An instance of the TPC PID maker - Double_t parTPC[]={47.,0.1,3.}; - AliTPCpidESD tpcPID(parTPC); - //File with the ITS clusters - TFile *itscf=TFile::Open(fileNameITSClusters); - if (!itscf->IsOpen()) { - cerr<<"Can't open "<GetLoader("ITSLoader"); + if (itsl == 0x0) { + cerr<<"AliESDtest.C : Can not get the ITS loader"<LoadRecPoints("read"); + + AliITS *dITS = (AliITS*)gAlice->GetDetector("ITS"); + if (!dITS) { + cerr<<"AliESDtest.C : Can not find the ITS detector !"<Get("AliITSgeom"); - if (!geom) {cerr<<"Can't get AliITSgeom !\n"; return 5;} + AliITSgeom *geom = dITS->GetITSgeom(); //An instance of the ITS tracker AliITStrackerV2 itsTracker(geom); //An instance of the ITS PID maker - Double_t parITS[]={34.,0.12,3.}; + Double_t parITS[]={34.,0.15,10.}; AliITSpidESD itsPID(parITS); - //File with the TRD clusters - TFile *trdcf=TFile::Open(fileNameTRDClusters); - if (!trdcf->IsOpen()) { - cerr<<"Can't open "<GetLoader("TPCLoader"); + if (tpcl == 0x0) { + cerr<<"AliESDtest.C : can not get the TPC loader"<LoadRecPoints("read"); + + rl->CdGAFile(); + AliTPCParam *par=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60"); + if (!par) { + cerr<<"TPC parameters have not been found !\n"; return 6; } + + //An instance of the TPC tracker + AliTPCtracker tpcTracker(par); + + //An instance of the TPC PID maker + Double_t parTPC[]={47.,0.10,10.}; + AliTPCpidESD tpcPID(parTPC); + + +/**** The TRD corner ********************/ + + AliLoader* trdl = rl->GetLoader("TRDLoader"); + if (trdl == 0x0) { + cerr<<"AliESDtest.C : can not get the TRD loader"<LoadRecPoints("read"); //An instance of the TRD tracker - AliTRDtracker trdTracker(trdcf); + rl->CdGAFile(); + AliTRDtracker trdTracker(gFile); //galice.root file +/* //An instance of the TRD PID maker - AliTRDPartID* trdPID = AliTRDPartID::GetFromFile(); - if (!trdPID) return 7; + TFile* pidFile = TFile::Open("pid.root"); + if (!pidFile->IsOpen()) { + cerr << "Can't get pid.root !\n"; + return 7; + } + AliTRDPartID* trdPID = (AliTRDPartID*) pidFile->Get("AliTRDPartID"); + if (!trdPID) { + cerr << "Can't get PID object !\n"; + return 8; + } +*/ + + +/**** The TOF corner ********************/ + + AliLoader* tofl = rl->GetLoader("TOFLoader"); + if (tofl == 0x0) { + cerr<<"AliESDtest.C : can not get the TOF loader"<LoadDigits("read"); + + + //Instance of the TOF PID maker + Double_t parTOF[]={130.,5.}; + AliTOFpidESD tofPID(parTOF); + + + //rl->UnloadgAlice(); + TFile *ef=TFile::Open("AliESDs.root","RECREATE"); if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;} TStopwatch timer; Int_t rc=0; + if (nev>rl->GetNumberOfEvents()) nev=rl->GetNumberOfEvents(); //The loop over events for (Int_t i=0; iGetEvent(i); + + TTree *tpcTree=tpcl->TreeR(); + if (!tpcTree) { + cerr<<"Can't get the TPC cluster tree !\n"; + return 4; + } + tpcTracker.LoadClusters(tpcTree); rc+=tpcTracker.Clusters2Tracks(event); + + TTree *itsTree=itsl->TreeR(); + if (!itsTree) { + cerr<<"Can't get the TPC cluster tree !\n"; + return 4; + } + itsTracker.LoadClusters(itsTree); rc+=itsTracker.Clusters2Tracks(event); rc+=itsTracker.PropagateBack(event); itsTracker.UnloadClusters(); - itsPID.MakePID(event); rc+=tpcTracker.PropagateBack(event); tpcTracker.UnloadClusters(); - tpcPID.MakePID(event); - trdTracker.SetEventNumber(i); - trdcf->cd(); - trdTracker.LoadClusters(); + TTree *trdTree=trdl->TreeR(); + if (!trdTree) { + cerr<<"Can't get the TPC cluster tree !\n"; + return 4; + } + trdTracker.LoadClusters(trdTree); rc+=trdTracker.PropagateBack(event); trdTracker.UnloadClusters(); +/* for (Int_t iTrack = 0; iTrack < event->GetNumberOfTracks(); iTrack++) { AliESDtrack* track = event->GetTrack(iTrack); trdPID->MakePID(track); } +*/ + + TTree *tofTree=tofl->TreeD(); + if (!tofTree) { + cerr<<"Can't get the TOF cluster tree !\n"; + return 4; + } + tofPID.LoadClusters(tofTree); + tofPID.MakePID(event); + tofPID.UnloadClusters(); + //Here is the combined PID AliESDpid::MakePID(event); @@ -137,12 +265,13 @@ Int_t AliESDtest(Int_t nev=1, } timer.Stop(); timer.Print(); - trdcf->Close(); - delete geom; - itscf->Close(); + //pidFile->Close(); + delete par; - tpccf->Close(); + ef->Close(); + delete rl; + return rc; } diff --git a/STEER/AliESDtrack.cxx b/STEER/AliESDtrack.cxx index 8a4e2fbc794..3e6f9bc0c1c 100644 --- a/STEER/AliESDtrack.cxx +++ b/STEER/AliESDtrack.cxx @@ -156,6 +156,17 @@ void AliESDtrack::SetIntegratedTimes(const Double_t *times) { for (Int_t i=0; iAddHitList(fHits); + if (gAlice==0) { + Fatal("AliTOF","gAlice==0 !"); + } + if (gAlice->GetHitLists()) + gAlice->AddHitList(fHits); + else Error("AliTOF","gAlice->GetHitLists()==0"); + fIshunt = 0; fSDigits = new TClonesArray("AliTOFSDigit", 1000); fDigits = new TClonesArray("AliTOFdigit", 1000); diff --git a/TOF/AliTOFSDigitizer.cxx b/TOF/AliTOFSDigitizer.cxx index 79845a73f1e..ae0f922c2d3 100644 --- a/TOF/AliTOFSDigitizer.cxx +++ b/TOF/AliTOFSDigitizer.cxx @@ -178,6 +178,7 @@ void AliTOFSDigitizer::Exec(Option_t *verboseOption, Option_t *allEvents) { fRunLoader->LoadgAlice(); fRunLoader->LoadHeader(); + fRunLoader->LoadKinematics(); gAlice = fRunLoader->GetAliRun(); AliLoader* gime = fRunLoader->GetLoader("TOFLoader"); diff --git a/TOF/AliTOFSDigits2Digits.C b/TOF/AliTOFSDigits2Digits.C index 684fbcec11f..8de37e776da 100644 --- a/TOF/AliTOFSDigits2Digits.C +++ b/TOF/AliTOFSDigits2Digits.C @@ -1,111 +1,88 @@ -Int_t AliTOFSDigits2Digits(TString infileNameSDigits, Int_t firstEvent=0,Int_t nEvents=1, Int_t ndump=15) -{ - // Create TOF digits out of TOF sdigits (no merging implemented here). - // Input: infileNameSDigits --> TOF sdigits filename - // firstEvent --> first event to digitize - // nEvents --> number of events to digitize - // including the first and the last ones - // if nEvents==0 we sdigitize all events in infileNameSDigits - // ndump --> number of dumped digits - - // Author: F. Pierella (Bologna University) - // report problem to pierella@bo.infn.it - - // Use case: (start root) - //root [0] .L AliTOFSDigits2Digits.C - //root [1] AliTOFSDigits2Digits("fileWithTOFSdigits.root") +Int_t AliTOFSDigits2Digits(Int_t nev=5) { + // Adapted to the NewIO by I.Belikov (Jouri.Belikov@cern.ch) // number 3 is a legacy from AliDigit object - const Int_t kMAXDIGITS = 3; - - - // Dynamically link some shared libs - if (gClassTable->GetID("AliRun") < 0) { - gROOT->LoadMacro("loadlibs.C"); - loadlibs(); - } else { - delete gAlice; - gAlice = 0; - } - - Int_t rc=0; - - - TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(infileNameSDigits.Data()); - if(file){ - cout<<"headerFile already open \n"; - } - else { - // file is open in "update" mode - // in order to have a writable digits tree - if(!file)file=TFile::Open(infileNameSDigits.Data(),"update"); - } - - // Get AliRun object from file - if (!gAlice) { - gAlice = (AliRun*)file->Get("gAlice"); - if (gAlice) printf("AliRun object found on file\n"); - } - - if (nEvents == 0) nEvents = (Int_t) gAlice->TreeE()->GetEntries(); - - AliTOF * tof = (AliTOF *) gAlice->GetDetector("TOF") ; - - if (!tof) { - cout << " No TOF detector found" << endl; - rc = 2; - return rc; - } - - Int_t totndig=0; // total number of digits - Int_t tottracks=0; // total number of tracks contributing to totndig - Int_t upperBoundEvNum=firstEvent+nEvents; - for (Int_t ievent = firstEvent; ievent < upperBoundEvNum; ievent++) { - - gAlice->GetEvent(ievent); - if (gAlice->TreeD () == 0) - gAlice->MakeTree ("D"); - //Make branches - char branchname[20]; - sprintf (branchname, "%s", tof->GetName ()); - //Make branch for digits - tof->MakeBranch ("D"); - - - // get the TOF branch in TreeS for current event - char tname[100]; sprintf(tname,"TreeS%d",ievent); + if (gAlice) { + delete gAlice->GetRunLoader(); + delete gAlice;//if everything was OK here it is already NULL + gAlice = 0x0; + } + + AliRunLoader* rl = AliRunLoader::Open("galice.root"); + if (rl == 0x0) { + cerr<<"Can not open session"<LoadgAlice()) { + cerr<<"Error occured while loading gAlice"<GetLoader("TOFLoader"); + if (tofl == 0x0) { + cerr<<"Can not get the TOF Loader"<GetAliRun(); + if (!gAlice) { + cerr<<"Can't get gAlice !\n"; + return 1; + } + + tofl->LoadSDigits("read"); + tofl->LoadDigits("recreate"); + + Int_t totndig=0; // total number of digits + Int_t tottracks=0; // total number of tracks contributing to totndig + + if (nev>rl->GetNumberOfEvents()) nev=rl->GetNumberOfEvents(); + + TClonesArray *fSDigits=new TClonesArray("AliTOFSDigit", 1000); + TClonesArray *fDigits =new TClonesArray("AliTOFdigit", 1000); + TClonesArray &da=*fDigits; + + for (Int_t ievent = 0; ievent < nev; ievent++) { + rl->GetEvent(ievent); + + TTree *sTree=tofl->TreeS(); + if (sTree == 0) { + cerr<<"Can't get the sdigit tree !\n"; + return 1; + } + TBranch *branch=sTree->GetBranch("TOF"); + if (!branch) { + cerr<<"Cant' get the branch !\n"; + return 1; + } + branch->SetAddress(&fSDigits); - TTree *sdigitstree=(TTree*)file->Get(tname); - - TClonesArray * fSDigits= new TClonesArray("AliTOFSDigit", 1000); - sdigitstree->GetBranch("TOF")->SetAddress(&fSDigits); - - Int_t nEntries = sdigitstree->GetEntries(); - //cout << nEntries << endl; - - // Loop through all entries in the tree - Int_t nbytes; - + TTree *dTree=tofl->TreeD(); + if (dTree == 0) { + tofl->MakeTree("D"); + dTree=tofl->TreeD(); + } + branch=dTree->GetBranch("TOF"); + if (!branch) dTree->Branch("TOF",&fDigits); + else branch->SetAddress(&fDigits); + + Int_t nEntries = sTree->GetEntries(); for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) { + sTree->GetEvent(iEntry); - // Import the tree - nbytes += sdigitstree->GetEvent(iEntry); - - // Get the number of sdigits Int_t ndig = fSDigits->GetEntriesFast(); - cout << "------------------------------------" << endl; + cout << "--------------------------------" << endl; cout << "Found " << ndig << " TOF SDigits for event " << ievent << endl; - cout << "------------------------------------------------------------" << endl; + cout << "------------------------------------------------------" << endl; - // start loop on sdigits for (Int_t k = 0; k < ndig; k++) { - Int_t vol[5]; // location for a digit // Get the information for this digit - AliTOFSDigit *tofsdigit = (AliTOFSDigit *) fSDigits->UncheckedAt(k); + AliTOFSDigit *tofsdigit = (AliTOFSDigit *)fSDigits->UncheckedAt(k); Int_t nslot=tofsdigit->GetNDigits(); // get the number of slots // for current sdigit @@ -125,27 +102,22 @@ Int_t AliTOFSDigits2Digits(TString infileNameSDigits, Int_t firstEvent=0,Int_t n //--------------------- QA section ---------------------- // in the while, I perform QA - Bool_t isSDigitBad = (sector<1 || sector>18 || plate<1 || plate >5 || padz<1 || padz>2 || padx<1 || padx>48); + Bool_t isSDigitBad = (sector<1 || sector>18 || + plate<1 || plate >5 || + padz<1 || padz>2 || + padx<1 || padx>48); if (isSDigitBad) { cout << " strange sdigit found" << endl; - rc = 3; - return rc; + return 3; } //------------------------------------------------------- - //------------------- Dump section ---------------------- - if(kGetTdc(islot); digit[0]=tdc; Float_t adc=tofsdigit->GetAdc(islot); digit[1]=adc; @@ -164,32 +136,37 @@ Int_t AliTOFSDigits2Digits(TString infileNameSDigits, Int_t firstEvent=0,Int_t n } // adding a TOF digit for each slot - tof->AddDigit(tracknum, vol, digit); + { + Int_t ndigits=da.GetEntriesFast(); + //cerr<Clear(); } // end loop on entries - // free used memory - fSDigits->Clear(); - fSDigits=0; + dTree->Fill(); + tofl->WriteDigits("OVERWRITE"); - gAlice->TreeD()->Reset(); - gAlice->TreeD()->Fill(); - //gAlice->TreeS()->Write(0,TObject::kOverwrite) ; - gAlice->TreeD()->AutoSave(); + // free used memory + fDigits->Clear(); } // end loop on events + delete fSDigits; + delete fDigits; + cout << "----------------------------------------------------------" << endl; cout << " Summary" << endl; cout << "contributing tracks to " << totndig << " digits: " << tottracks << endl; cout << "----------------------------------------------------------" << endl; - return rc; + return 0; } diff --git a/TOF/AliTOFpidESD.cxx b/TOF/AliTOFpidESD.cxx new file mode 100644 index 00000000000..2506b790e41 --- /dev/null +++ b/TOF/AliTOFpidESD.cxx @@ -0,0 +1,451 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +//----------------------------------------------------------------- +// Implementation of the TOF PID class +// Very naive one... Should be made better by the detector experts... +// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch +//----------------------------------------------------------------- +#include "TFile.h" +#include "TTree.h" +#include "TClonesArray.h" + +#include "AliTOFpidESD.h" +#include "AliESD.h" +#include "AliESDtrack.h" +#include "AliTOFdigit.h" + +#include +#include +#include "../STRUCT/AliBODY.h" +#include "../STRUCT/AliFRAMEv2.h" +#include "AliTOFv2FHoles.h" + + +ClassImp(AliTOFpidESD) + +static Int_t InitGeo() { + //gSystem->Load("libgeant321"); + //new TGeant3("C++ Interface to Geant3"); + + AliBODY *BODY = new AliBODY("BODY", "Alice envelop"); + BODY->CreateGeometry(); + + AliFRAMEv2 *FRAME = new AliFRAMEv2("FRAME", "Space Frame"); + FRAME->SetHoles(1); + FRAME->CreateGeometry(); + + AliTOF *TOF = new AliTOFv2FHoles("TOF", "TOF with Holes"); + TOF->CreateGeometry(); + + return 0; +} + +//_________________________________________________________________________ +AliTOFpidESD::AliTOFpidESD(Double_t *param) throw (const Char_t *) { + // + // The main constructor + // + if (InitGeo()) throw "AliTOFpidESD: can not initialize the geomtry !\n"; + + fR=376.; fDy=2.5; fDz=3.5; fN=0; fEventN=0; + + fSigma=param[0]; + fRange=param[1]; + +} + +static Int_t DtoM(Int_t *dig, Float_t *g) { + const Int_t MAX=13; + Int_t lnam[MAX],lnum[MAX]; + + extern TVirtualMC *gMC; + + TGeant3 *geant3=(TGeant3*)gMC; + if (!geant3) {cerr<<"no geant3 found !\n"; return 1;} + + strncpy((Char_t*)(lnam+0),"ALIC",4); lnum[0]=1; + strncpy((Char_t*)(lnam+1),"B077",4); lnum[1]=1; + + //sector + switch (dig[0]) { + case 1: + strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=9; + strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1; + break; + case 2: + strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=10; + strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1; + break; + case 3: + strncpy((Char_t*)(lnam+2),"B074",4); lnum[2]=4; + strncpy((Char_t*)(lnam+3),"BTO2",4); lnum[3]=1; + if (dig[1]>3) lnum[3]=2; + break; + case 4: + strncpy((Char_t*)(lnam+2),"B074",4); lnum[2]=5; + strncpy((Char_t*)(lnam+3),"BTO2",4); lnum[3]=1; + if (dig[1]>3) lnum[3]=2; + break; + case 5: + strncpy((Char_t*)(lnam+2),"B074",4); lnum[2]=1; + strncpy((Char_t*)(lnam+3),"BTO2",4); lnum[3]=1; + if (dig[1]>3) lnum[3]=2; + break; + case 6: + strncpy((Char_t*)(lnam+2),"B074",4); lnum[2]=2; + strncpy((Char_t*)(lnam+3),"BTO2",4); lnum[3]=1; + if (dig[1]>3) lnum[3]=2; + break; + case 7: + strncpy((Char_t*)(lnam+2),"B074",4); lnum[2]=3; + strncpy((Char_t*)(lnam+3),"BTO2",4); lnum[3]=1; + if (dig[1]>3) lnum[3]=2; + break; + case 8: + strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=1; + strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1; + break; + case 9: + strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=2; + strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1; + break; + case 10: + strncpy((Char_t*)(lnam+2),"B075",4); lnum[2]=1; + strncpy((Char_t*)(lnam+3),"BTO3",4); lnum[3]=1; + if (dig[1]>4) lnum[3]=2; + break; + case 11: + strncpy((Char_t*)(lnam+2),"B075",4); lnum[2]=2; + strncpy((Char_t*)(lnam+3),"BTO3",4); lnum[3]=1; + if (dig[1]>4) lnum[3]=2; + break; + case 12: + strncpy((Char_t*)(lnam+2),"B075",4); lnum[2]=3; + strncpy((Char_t*)(lnam+3),"BTO3",4); lnum[3]=1; + if (dig[1]>4) lnum[3]=2; + break; + case 13: + strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=3; + strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1; + break; + case 14: + strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=4; + strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1; + break; + case 15: + strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=5; + strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1; + break; + case 16: + strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=6; + strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1; + break; + case 17: + strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=7; + strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1; + break; + case 18: + strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=8; + strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1; + break; + default: + cerr<<"Wrong sector number : "<=8 && dig[0]<=9) lnum[4]=2; + if (dig[0]>=13 && dig[0]<=18) lnum[4]=2; + strncpy((Char_t*)(lnam+5),"FLTC",4); lnum[5]=0; + break; + case 2: + strncpy((Char_t*)(lnam+4),"FTOB",4); lnum[4]=1; + if (dig[0]==1 || dig[0]==2 ) lnum[4]=2; + if (dig[0]>=8 && dig[0]<=9) lnum[4]=2; + if (dig[0]>=13 && dig[0]<=18) lnum[4]=2; + strncpy((Char_t*)(lnam+5),"FLTB",4); lnum[5]=0; + break; + case 3: + strncpy((Char_t*)(lnam+4),"FTOA",4); lnum[4]=0; + strncpy((Char_t*)(lnam+5),"FLTA",4); lnum[5]=0; + break; + case 4: + strncpy((Char_t*)(lnam+4),"FTOB",4); lnum[4]=1; + strncpy((Char_t*)(lnam+5),"FLTB",4); lnum[5]=0; + break; + case 5: + strncpy((Char_t*)(lnam+4),"FTOC",4); lnum[4]=1; + strncpy((Char_t*)(lnam+5),"FLTC",4); lnum[5]=0; + break; + default: + cerr<<"Wrong module number : "<19) { + cerr<<"Wrong strip number : "<2) { + cerr<<"Wrong z-division number : "<48) { + cerr<<"Wrong x-division number : "<Gcvolu(); gcvolu->nlevel=0; + Int_t err=geant3->Glvolu(11,lnam,lnum); //11-th level + if (err) {cout<Gdtom(l,g,1); + return 0; +} + +Int_t AliTOFpidESD::LoadClusters(const TFile *df) { + //-------------------------------------------------------------------- + //This function loads the TOF clusters + //-------------------------------------------------------------------- + if (!((TFile *)df)->IsOpen()) { + Error("LoadClusters","file with the TOF digits has not been open !\n"); + return 1; + } + + Char_t name[100]; + sprintf(name,"TreeD%d",GetEventNumber()); + TTree *dTree=(TTree*)((TFile *)df)->Get(name); + if (!dTree) { + Error("LoadClusters"," can't get the tree with digits !\n"); + return 1; + } + TBranch *branch=dTree->GetBranch("TOF"); + if (!branch) { + Error("LoadClusters"," can't get the branch with the TOF digits !\n"); + return 1; + } + + TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy; + branch->SetAddress(&digits); + + dTree->GetEvent(0); + Int_t nd=digits->GetEntriesFast(); + Info("LoadClusters","number of digits: %d",nd); + + for (Int_t i=0; iUncheckedAt(i); + Int_t dig[5]; Float_t g[3]; + dig[0]=d->GetSector(); + dig[1]=d->GetPlate(); + dig[2]=d->GetStrip(); + dig[3]=d->GetPadz(); + dig[4]=d->GetPadx(); + + DtoM(dig,g); + + Double_t h[5]; + h[0]=TMath::Sqrt(g[0]*g[0]+g[1]*g[1]); + h[1]=TMath::ATan2(g[1],g[0]); h[2]=g[2]; + h[3]=d->GetTdc(); h[4]=d->GetAdc(); + + AliTOFcluster *cl=new AliTOFcluster(h,d->GetTracks()); + InsertCluster(cl); + } + + delete dTree; + return 0; +} + +Int_t AliTOFpidESD::LoadClusters(TTree *dTree) { + //-------------------------------------------------------------------- + //This function loads the TOF clusters + //-------------------------------------------------------------------- + TBranch *branch=dTree->GetBranch("TOF"); + if (!branch) { + Error("LoadClusters"," can't get the branch with the TOF digits !\n"); + return 1; + } + + TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy; + branch->SetAddress(&digits); + + dTree->GetEvent(0); + Int_t nd=digits->GetEntriesFast(); + Info("LoadClusters","number of digits: %d",nd); + + for (Int_t i=0; iUncheckedAt(i); + Int_t dig[5]; Float_t g[3]; + dig[0]=d->GetSector(); + dig[1]=d->GetPlate(); + dig[2]=d->GetStrip(); + dig[3]=d->GetPadz(); + dig[4]=d->GetPadx(); + + DtoM(dig,g); + + Double_t h[5]; + h[0]=TMath::Sqrt(g[0]*g[0]+g[1]*g[1]); + h[1]=TMath::ATan2(g[1],g[0]); h[2]=g[2]; + h[3]=d->GetTdc(); h[4]=d->GetAdc(); + + AliTOFcluster *cl=new AliTOFcluster(h,d->GetTracks()); + InsertCluster(cl); + } + + return 0; +} + +void AliTOFpidESD::UnloadClusters() { + //-------------------------------------------------------------------- + //This function unloads TOF clusters + //-------------------------------------------------------------------- + for (Int_t i=0; iGetZ()); + memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTOFcluster*)); + fClusters[i]=c; fN++; + + return 0; +} + +Int_t AliTOFpidESD::FindClusterIndex(Double_t z) const { + //-------------------------------------------------------------------- + // This function returns the index of the nearest cluster + //-------------------------------------------------------------------- + if (fN==0) return 0; + if (z <= fClusters[0]->GetZ()) return 0; + if (z > fClusters[fN-1]->GetZ()) return fN; + Int_t b=0, e=fN-1, m=(b+e)/2; + for (; b fClusters[m]->GetZ()) b=m+1; + else e=m; + } + return m; +} + +//_________________________________________________________________________ +Int_t AliTOFpidESD::MakePID(AliESD *event) +{ + // + // This function calculates the "detector response" PID probabilities + // Just for a bare hint... + + static const Double_t masses[]={ + 0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613 + }; + + Int_t ntrk=event->GetNumberOfTracks(); + Int_t nmatch=0; + for (Int_t i=0; iGetTrack(i); + + if ((t->GetStatus()&AliESDtrack::kTRDout)==0) continue; + if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue; + + Double_t x,par[5]; t->GetExternalParameters(x,par); + Double_t cov[15]; t->GetExternalCovariance(cov); + + Double_t dphi2=3*3*(cov[0] + fDy*fDy/12.)/fR/fR; + Double_t dz=3*TMath::Sqrt(cov[2] + fDz*fDz/12.); + + Double_t phi=TMath::ATan2(par[0],x) + t->GetAlpha(); + if (phi<-TMath::Pi()) phi+=2*TMath::Pi(); + if (phi>=TMath::Pi()) phi-=2*TMath::Pi(); + Double_t z=par[1]; + + Double_t d2max=1000.; + Int_t index=-1; + for (Int_t k=FindClusterIndex(z-dz); kGetZ() > z+dz) break; + if (c->IsUsed()) continue; + + Double_t dphi=TMath::Abs(c->GetPhi()-phi); + if (dphi>TMath::Pi()) dphi-=TMath::Pi(); + if (dphi*dphi > dphi2) continue; + + Double_t d2=dphi*dphi*fR*fR + (c->GetZ()-z)*(c->GetZ()-z); + if (d2 > d2max) continue; + + d2max=d2; + index=k; + } + + if (index<0) { + //Info("MakePID","matching failed ! %d",TMath::Abs(t->GetLabel())); + continue; + } + + nmatch++; + + AliTOFcluster *c=fClusters[index]; + + Double_t time[10]; t->GetIntegratedTimes(time); + Double_t tof=50*c->GetTDC(); // in ps + + Double_t p[10]; + Double_t mom=t->GetP(); + for (Int_t j=0; j fRange*sigma) { + p[j]=TMath::Exp(-0.5*fRange*fRange); + continue; + } + p[j]=TMath::Exp(-0.5*(tof-time[j])*(tof-time[j])/(sigma*sigma)); + } + + //if (c->GetLabel(0)!=TMath::Abs(t->GetLabel())) continue; + + t->SetTOFsignal(tof); + t->SetTOFcluster(index); + t->SetTOFpid(p); + + c->Use(); + } + + Info("MakePID","Number of matched ESD track: %d",nmatch); + + return 0; +} diff --git a/TOF/AliTOFpidESD.h b/TOF/AliTOFpidESD.h new file mode 100644 index 00000000000..98f3ee67293 --- /dev/null +++ b/TOF/AliTOFpidESD.h @@ -0,0 +1,80 @@ +#ifndef ALITOFPIDESD_H +#define ALITOFPIDESD_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//------------------------------------------------------- +// TOF PID class +// A very naive design... Should be made better by the detector experts... +// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch +//------------------------------------------------------- + +#include + +class AliESD; +class TFile; +class TTree; + +const Int_t kMaxCluster=77777; //maximal number of the TOF clusters + +class AliTOFpidESD : public TObject { +public: + AliTOFpidESD(){fR=376.; fDy=2.5; fDz=3.5; fN=0; fEventN=0;} + AliTOFpidESD(Double_t *param) throw (const Char_t *); + ~AliTOFpidESD(){UnloadClusters();} + + Int_t MakePID(AliESD *event); + Int_t LoadClusters(const TFile *f); + Int_t LoadClusters(TTree *f); + void UnloadClusters(); + void SetEventNumber(Int_t n) {fEventN=n;} + + Int_t GetEventNumber() const {return fEventN;} + +public: + class AliTOFcluster { + public: + AliTOFcluster(Double_t *h, Int_t *l) { + fR=h[0]; fPhi=h[1]; fZ=h[2]; fTDC=h[3]; fADC=h[4]; + fLab[0]=l[0]; fLab[1]=l[1]; fLab[2]=l[2]; + } + void Use() {fADC=-fADC;} + + Double_t GetR() const {return fR;} + Double_t GetPhi() const {return fPhi;} + Double_t GetZ() const {return fZ;} + Double_t GetTDC() const {return fTDC;} + Double_t GetADC() const {return TMath::Abs(fADC);} + Int_t IsUsed() const {return (fADC<0) ? 1 : 0;} + Int_t GetLabel(Int_t n) const {return fLab[n];} + private: + Int_t fLab[3]; + Double_t fR; + Double_t fPhi; + Double_t fZ; + Double_t fTDC; + Double_t fADC; + }; + +private: + Int_t InsertCluster(AliTOFcluster*); + Int_t FindClusterIndex(Double_t z) const; + + Int_t fEventN; //event number + + Double_t fR; // mean readius of the TOF barrel + Double_t fDy; // size of the TOF cell in R*Phi + Double_t fDz; // size of the TOF cell in Z + + Double_t fSigma; // intrinsic TOF resolution + Double_t fRange; // one particle type PID range (in sigmas) + + Int_t fN; // number of the TOF clusters + AliTOFcluster *fClusters[kMaxCluster]; // pointers to the TOF clusters + + ClassDef(AliTOFpidESD,1) // TOF PID class +}; + +#endif + + diff --git a/TOF/TOFLinkDef.h b/TOF/TOFLinkDef.h index f930b36f693..c492470ba76 100644 --- a/TOF/TOFLinkDef.h +++ b/TOF/TOFLinkDef.h @@ -41,4 +41,6 @@ #pragma link C++ class AliTOFReconstructionerV2+; #pragma link C++ class AliTOFDigitMap+; +#pragma link C++ class AliTOFpidESD+; + #endif diff --git a/TOF/libTOF.pkg b/TOF/libTOF.pkg index ec4b8016e0c..7f71fe123f2 100644 --- a/TOF/libTOF.pkg +++ b/TOF/libTOF.pkg @@ -1,4 +1,4 @@ -SRCS = AliTOF.cxx AliTOFv0.cxx AliTOFv1.cxx AliTOFv2.cxx AliTOFv3.cxx AliTOFv4.cxx AliTOFv4T0.cxx AliTOFhit.cxx AliTOFhitT0.cxx AliTOFPID.cxx AliTOFT0.cxx AliTOFdigit.cxx AliTOFRawSector.cxx AliTOFRoc.cxx AliTOFRawDigit.cxx AliTOFDigitizer.cxx AliTOFSDigitizer.cxx AliTOFMerger.cxx AliTOFv2FHoles.cxx AliTOFSDigit.cxx AliTOFHitMap.cxx AliTOFConstants.cxx AliTOFPad.cxx AliTOFRecHit.cxx AliTOFTrack.cxx AliTOFReconstructioner.cxx AliTOFProb.cxx AliTOFTrackV2.cxx AliTOFReconstructionerV2.cxx AliTOFDigitMap.cxx +SRCS = AliTOF.cxx AliTOFv0.cxx AliTOFv1.cxx AliTOFv2.cxx AliTOFv3.cxx AliTOFv4.cxx AliTOFv4T0.cxx AliTOFhit.cxx AliTOFhitT0.cxx AliTOFPID.cxx AliTOFT0.cxx AliTOFdigit.cxx AliTOFRawSector.cxx AliTOFRoc.cxx AliTOFRawDigit.cxx AliTOFDigitizer.cxx AliTOFSDigitizer.cxx AliTOFMerger.cxx AliTOFv2FHoles.cxx AliTOFSDigit.cxx AliTOFHitMap.cxx AliTOFConstants.cxx AliTOFPad.cxx AliTOFRecHit.cxx AliTOFTrack.cxx AliTOFReconstructioner.cxx AliTOFProb.cxx AliTOFTrackV2.cxx AliTOFReconstructionerV2.cxx AliTOFDigitMap.cxx AliTOFpidESD.cxx HDRS:= $(SRCS:.cxx=.h) diff --git a/TPC/AliTPC.cxx b/TPC/AliTPC.cxx index f006ca43eca..4ab2d4e3092 100644 --- a/TPC/AliTPC.cxx +++ b/TPC/AliTPC.cxx @@ -73,8 +73,6 @@ #include "AliTPCTrackHits.h" #include "AliTPCTrackHitsV2.h" #include "AliTPCcluster.h" -#include "AliTPCclusterer.h" -#include "AliTPCtracker.h" #include "AliTrackReference.h" @@ -316,8 +314,8 @@ void AliTPC::Clusters2Tracks() //----------------------------------------------------------------- // This is a track finder. //----------------------------------------------------------------- - AliTPCtracker tracker(fTPCParam,0,fLoader->GetEventFolder()->GetName()); - tracker.Clusters2Tracks(); + Error("Clusters2Tracks", + "Dummy function ! Call AliTPCtracker::Clusters2Tracks(...) instead !"); } //_____________________________________________________________________________ @@ -839,13 +837,14 @@ void AliTPC::Digits2Clusters(Int_t eventnumber) //----------------------------------------------------------------- // This is a simple cluster finder. //----------------------------------------------------------------- - AliTPCclusterer::Digits2Clusters(fTPCParam, fLoader,eventnumber); + Error("Digits2Clusters", + "Dummy function ! Call AliTPCclusterer::Digits2Clusters(...) instead !"); } extern Double_t SigmaY2(Double_t, Double_t, Double_t); extern Double_t SigmaZ2(Double_t, Double_t); //_____________________________________________________________________________ -void AliTPC::Hits2Clusters(TFile *of, Int_t eventn) +void AliTPC::Hits2Clusters(Int_t eventn) { //-------------------------------------------------------- // TPC simple cluster generator from hits @@ -917,7 +916,7 @@ void AliTPC::Hits2Clusters(TFile *of, Int_t eventn) AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::fgkRunLoaderName); rl->CdGAFile(); - fTPCParam->Write(fTPCParam->GetTitle()); + //fTPCParam->Write(fTPCParam->GetTitle()); AliTPCClustersArray carray; carray.Setup(fTPCParam); diff --git a/TPC/AliTPC.h b/TPC/AliTPC.h index 3cdd502292b..03f6e10a711 100644 --- a/TPC/AliTPC.h +++ b/TPC/AliTPC.h @@ -69,7 +69,7 @@ public: virtual void BuildGeometry(); virtual void CreateGeometry() {} virtual void CreateMaterials(); - virtual void Hits2Clusters(TFile *of, Int_t eventn=0); + virtual void Hits2Clusters(Int_t eventn=0); virtual void Hits2ExactClustersSector(Int_t isec); // MI change calculate "exact" cluster position virtual void SDigits2Digits(){;} //MI change -cycling to production diff --git a/TPC/AliTPCComparison.C b/TPC/AliTPCComparison.C index 710497f7b7f..44091e134d5 100644 --- a/TPC/AliTPCComparison.C +++ b/TPC/AliTPCComparison.C @@ -8,11 +8,13 @@ * with several nice improvements by: M.Ivanov, GSI, m.ivanov@gsi.de * ****************************************************************************/ -#ifndef __CINT__ +#if !defined(__CINT__) || defined(__MAKECINT__) #include "alles.h" #include "AliTPCtracker.h" #include "AliRunLoader.h" + #include "AliTPCLoader.h" #include "AliMagF.h" + #include "AliRun.h" #endif struct GoodTrackTPC { @@ -22,12 +24,14 @@ struct GoodTrackTPC { Float_t x,y,z; }; -Int_t good_tracks_tpc(GoodTrackTPC *gt, const Int_t max, const char* evfoldname); +Int_t +good_tracks_tpc(GoodTrackTPC *gt, const Int_t max, const char* evfoldname); -Int_t AliTPCComparison(Int_t event=0) { +extern AliRun *gAlice; - if (gAlice) - { +Int_t AliTPCComparison() { + + if (gAlice) { delete gAlice->GetRunLoader(); delete gAlice;//if everything was OK here it is already NULL gAlice = 0x0; @@ -36,34 +40,36 @@ Int_t AliTPCComparison(Int_t event=0) { cerr<<"Doing comparison...\n"; const Int_t MAX=20000; - Int_t good_tracks_tpc(GoodTrackTPC *gt, const Int_t max, const char* evfoldname = AliConfig::fgkDefaultEventFolderName);//declaration only + Int_t good_tracks_tpc( + GoodTrackTPC *gt, const Int_t max, + const char* evfoldname = AliConfig::fgkDefaultEventFolderName + );//declaration only gBenchmark->Start("AliTPCComparison"); AliRunLoader *rl = AliRunLoader::Open("galice.root","COMPARISON"); - if (!rl) - { + if (!rl) { cerr<<"Can't start sesion !\n"; return 1; - } + } rl->LoadgAlice(); if (rl->GetAliRun()) - AliKalmanTrack::SetConvConst(1000/0.299792458/rl->GetAliRun()->Field()->SolenoidField()); - else - { + AliKalmanTrack::SetConvConst( + 1000/0.299792458/rl->GetAliRun()->Field()->SolenoidField() + ); + else { cerr<<"AliTPCComparison.C :Can't get AliRun !\n"; return 1; - } + } rl->UnloadgAlice(); - AliLoader * tpcl = rl->GetLoader("TPCLoader"); - if (tpcl == 0x0) - { + AliTPCLoader * tpcl = (AliTPCLoader *)rl->GetLoader("TPCLoader"); + if (tpcl == 0x0) { cerr<<"AliTPCComparison.C : Can not find TPCLoader\n"; delete rl; return 3; - } + } Int_t nentr=0,i=0; TObjArray tarray(MAX); { /*Load tracks*/ @@ -76,13 +82,12 @@ Int_t AliTPCComparison(Int_t event=0) { TBranch *tbranch=tracktree->GetBranch("tracks"); nentr=(Int_t)tracktree->GetEntries(); AliTPCtrack *iotrack=0; - for (i=0; iSetAddress(&iotrack); tracktree->GetEvent(i); tarray.AddLast(iotrack); - } + } tpcl->UnloadTracks(); } @@ -118,7 +123,6 @@ Int_t AliTPCComparison(Int_t event=0) { } - TH1F *hp=new TH1F("hp","PHI resolution",50,-20.,20.); hp->SetFillColor(4); TH1F *hl=new TH1F("hl","LAMBDA resolution",50,-20,20);hl->SetFillColor(4); TH1F *hpt=new TH1F("hpt","Relative Pt resolution",30,-10.,10.); @@ -311,23 +315,22 @@ Int_t AliTPCComparison(Int_t event=0) { } -Int_t good_tracks_tpc(GoodTrackTPC *gt, const Int_t max, const char* evfoldname) { +Int_t +good_tracks_tpc(GoodTrackTPC *gt, const Int_t max, const char* evfoldname) { Int_t nt=0; - AliRunLoader* rl = AliRunLoader::GetRunLoader(evfoldname); - if (rl == 0x0) - { + AliRunLoader *rl = AliRunLoader::GetRunLoader(evfoldname); + if (rl == 0x0) { ::Fatal("AliTPCComparison.C::good_tracks_tpc", "Can not find Run Loader in Folder Named %s", evfoldname); - } - AliLoader * tpcl = rl->GetLoader("TPCLoader"); - if (tpcl == 0x0) - { + } + AliTPCLoader *tpcl = (AliTPCLoader *)rl->GetLoader("TPCLoader"); + if (tpcl == 0x0) { cerr<<"AliTPCHits2Digits.C : Can not find TPCLoader\n"; delete rl; return 0; - } + } rl->LoadgAlice(); diff --git a/TPC/AliTPCFindClusters.C b/TPC/AliTPCFindClusters.C index 15267160720..a5dad6b5389 100644 --- a/TPC/AliTPCFindClusters.C +++ b/TPC/AliTPCFindClusters.C @@ -2,47 +2,46 @@ * Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch * ****************************************************************************/ -#ifndef __CINT__ +#if !defined(__CINT__) || defined(__MAKECINT__) #include + #include "AliRun.h" + #include "AliRunLoader.h" + #include "AliTPCLoader.h" #include "AliTPCv1.h" - #include "AliTPCv2.h" #include "AliTPCParam.h" + #include "AliTPCclusterer.h" - #include "TFile.h" + #include "TTree.h" #include "TStopwatch.h" #endif -Int_t AliTPCFindClusters(Int_t N=-1) - { - - if (gAlice) - { +extern AliRun *gAlice; + +Int_t AliTPCFindClusters(Int_t nev=5) { + + if (gAlice) { delete gAlice->GetRunLoader(); delete gAlice;//if everything was OK here it is already NULL gAlice = 0x0; - } + } AliRunLoader* rl = AliRunLoader::Open("galice.root"); - if (rl == 0x0) - { + if (rl == 0x0) { cerr<<"Can not open session"<LoadgAlice()) - { + if (rl->LoadgAlice()) { cerr<<"Error occured while l"<GetAliRun()->Field()->SolenoidField()); + } - tpcl = (AliTPCLoader*)rl->GetLoader("TPCLoader"); - if (tpcl == 0x0) - { + AliTPCLoader *tpcl = (AliTPCLoader*)rl->GetLoader("TPCLoader"); + if (tpcl == 0x0) { cerr<<"Can not get TPC Loader"<GetAliRun(); if (!gAlice) { @@ -56,66 +55,57 @@ Int_t AliTPCFindClusters(Int_t N=-1) rl->CdGAFile(); AliTPCParam *dig=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60"); - if (!dig) - { - dig=(AliTPCParam *)gDirectory->Get("75x40_100x60"); - if (!param) - { + if (!dig) { cerr<<"TPC parameters have not been found !\n"; return 1; - } - else - { - cout<<"TPC 75x40_100x60 geometry found"<GetNumberOfEvents(); - else n=N; + } + + if (nev>rl->GetNumberOfEvents()) nev=rl->GetNumberOfEvents(); + tpcl->LoadRecPoints("recreate"); + if (ver==1) tpcl->LoadHits("read"); + else tpcl->LoadDigits("read"); + TStopwatch timer; - switch (ver) { - case 1: + if (ver==1) { cerr<<"Making clusters...\n"; - { - AliTPCv1 &tpc=*((AliTPCv1*)TPC); - tpc.SetParam(dig); timer.Start(); cwd->cd(); - tpc.SetLoader(tpcl); - tpcl->LoadHits("read"); - tpcl->LoadRecPoints("recreate"); - for(Int_t i=0;iLoadKinematics(); + timer.Start(); + for(Int_t i=0;iGetEvent(i); - tpc.Hits2Clusters(out,i); - } + tpc.Hits2Clusters(i); } - break; - case 2: + } else if (ver==2) { cerr<<"Looking for clusters...\n"; - { - // delete gAlice; gAlice=0; - AliTPCv2 * tpc = new AliTPCv2(); - tpc->SetLoader(tpcl); - tpcl->LoadDigits("read"); - tpcl->LoadRecPoints("recreate"); - - tpc->SetParam(dig); timer.Start(); - for (Int_t i=0;iDigits2Clusters(i); - //AliTPCclusterer::Digits2Clusters(dig, out, i); - } - delete tpc; + rl->GetEvent(i); + + TTree *out=tpcl->TreeR(); + if (!out) { + tpcl->MakeTree("R"); + out=tpcl->TreeR(); + } + TTree *in=tpcl->TreeD(); + if (!in) { + cerr<<"Can't get digits tree !\n"; + return 4; + } + + clusterer.Digits2Clusters(in,out); + + tpcl->WriteRecPoints("OVERWRITE"); } - break; - default: + delete dummy; + delete dig; + } else { cerr<<"Invalid TPC version !\n"; delete rl; return 5; @@ -124,5 +114,6 @@ Int_t AliTPCFindClusters(Int_t N=-1) timer.Stop(); timer.Print(); delete rl; + return 0; } diff --git a/TPC/AliTPCFindTracks.C b/TPC/AliTPCFindTracks.C index 3c67d7e3591..aa4449b3cc2 100644 --- a/TPC/AliTPCFindTracks.C +++ b/TPC/AliTPCFindTracks.C @@ -2,90 +2,93 @@ * Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch * ****************************************************************************/ -#ifndef __CINT__ +#if !defined(__CINT__) || defined(__MAKECINT__) #include #include "AliTPCParam.h" #include "AliTPCtracker.h" + #include "AliRun.h" + #include "AliMagF.h" + #include "AliRunLoader.h" + #include "AliTPCLoader.h" #include "TFile.h" #include "TStopwatch.h" #endif -Int_t AliTPCFindTracks(Int_t N=-1) { +extern AliRun *gAlice; + +Int_t AliTPCFindTracks(Int_t nev=5) { cerr<<"Looking for tracks...\n"; - if (gAlice) - { + if (gAlice) { delete gAlice->GetRunLoader(); delete gAlice; gAlice = 0x0; - } - - rl = AliRunLoader::Open("galice.root"); - if (rl == 0x0) - { + } + + AliRunLoader *rl = AliRunLoader::Open("galice.root"); + if (rl == 0x0) { cerr<<"Can not open session"<GetLoader("TPCLoader"); - if (tpcl == 0x0) - { + } + + AliTPCLoader *tpcl = (AliTPCLoader*)rl->GetLoader("TPCLoader"); + if (tpcl == 0x0) { cerr<<"Can not get TPC Loader"<LoadgAlice()) - { - cerr<<"Error occured while l"<LoadgAlice()) { + cerr<<"Error occured while loading gAlice"<GetAliRun()->Field()->SolenoidField()); + } + AliKalmanTrack::SetConvConst( + 1000/0.299792458/rl->GetAliRun()->Field()->SolenoidField() + ); rl->CdGAFile(); AliTPCParam *dig=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60"); - if (!dig) - { - dig=(AliTPCParam *)gDirectory->Get("75x40_100x60"); - if (!param) - { + if (!dig) { cerr<<"TPC parameters have not been found !\n"; return 1; - } - else - { - cout<<"TPC 75x40_100x60 geometry found"<UnloadgAlice() - + } + + rl->UnloadgAlice(); + tpcl->LoadRecPoints("read"); tpcl->LoadTracks("recreate"); - Int_t eventn; - if (N<=0) - { - eventn = rl->GetNumberOfEvents(); - rl->UnloadHeader(); - } - else - eventn = N; + if (nev>rl->GetNumberOfEvents()) nev=rl->GetNumberOfEvents(); TStopwatch timer; Int_t rc=0; - for (Int_t i=0;iClusters2Tracks(); - delete tracker; + rl->GetEvent(i); + + TTree *in=tpcl->TreeR(); + if (!in) { + cerr<<"Can't get clusters tree !\n"; + return 4; + } + + TTree *out=tpcl->TreeT(); + if (!out) { + tpcl->MakeTree("T"); + out=tpcl->TreeT(); + } + + rc=tracker.Clusters2Tracks(in,out); + + tpcl->WriteTracks("OVERWRITE"); } + timer.Stop(); timer.Print(); delete dig; //Thanks to Mariana Bondila + delete rl; + return rc; } diff --git a/TPC/AliTPCHits2Digits.C b/TPC/AliTPCHits2Digits.C index 529900a3ea8..4ac823f3cfc 100644 --- a/TPC/AliTPCHits2Digits.C +++ b/TPC/AliTPCHits2Digits.C @@ -1,22 +1,29 @@ -Int_t AliTPCHits2Digits(Int_t nevent=1) -{ +#if !defined(__CINT__) || defined(__MAKECINT__) + #include - // new version by J.Belikov + #include "AliRun.h" + #include "AliRunLoader.h" + #include "AliLoader.h" + #include "AliTPC.h" + #include "TStopwatch.h" +#endif + +extern AliRun *gAlice; + +Int_t AliTPCHits2Digits(Int_t nev=5) { // Connect the Root Galice file containing Geometry, Kine and Hits - if (gAlice) - { + if (gAlice) { delete gAlice->GetRunLoader(); delete gAlice;//if everything was OK here it is already NULL gAlice = 0x0; - } + } AliRunLoader *rl = AliRunLoader::Open("galice.root","Event","update"); - if (!rl) - { - cerr<<"Can't load RunLoader from "<GetEvent(0); AliTPC *TPC = (AliTPC*)gAlice->GetDetector("TPC"); AliLoader * tpcl = rl->GetLoader("TPCLoader"); - if ((TPC == 0x0) || (tpcl == 0x0)) - { + if ((TPC == 0x0) || (tpcl == 0x0)) { cerr<<"AliTPCHits2Digits.C : Can not find TPC or TPCLoader\n"; -// delete rl; + delete rl; return 3; - } + } tpcl->LoadHits("READ"); tpcl->LoadDigits("recreate"); TStopwatch timer; timer.Start(); - // uncomment below lines to set sectors active + // uncomment below lines to set sectors active // Int_t sec[10]={0,1,2,3,4,5,6,7,8,9}; // TPC->SetActiveSectors(sec,10); - for(Int_t eventn =0;eventnGetEvent(eventn); + for (Int_t i=0; iGetEvent(i); TPC->SetActiveSectors(); // all sectors set active - for (Int_t i=0;i<72;i++) if (TPC->IsSectorActive(i)) printf("%d\t",i); - TPC->Hits2Digits(eventn); + TPC->Hits2Digits(i); } delete rl; @@ -63,5 +65,5 @@ Int_t AliTPCHits2Digits(Int_t nevent=1) timer.Print(); return 0; -}; +} diff --git a/TPC/AliTPCclusterer.cxx b/TPC/AliTPCclusterer.cxx index aeb503fd5cd..a064fcd3480 100644 --- a/TPC/AliTPCclusterer.cxx +++ b/TPC/AliTPCclusterer.cxx @@ -25,15 +25,21 @@ #include "AliTPCcluster.h" #include #include -#include "AliTPCClustersArray.h" #include "AliTPCClustersRow.h" #include "AliDigits.h" #include "AliSimDigits.h" #include "AliTPCParam.h" -#include #include -#include "AliRunLoader.h" -#include "AliLoader.h" + +ClassImp(AliTPCclusterer) + +AliTPCclusterer::AliTPCclusterer(const AliTPCParam *par) { +//------------------------------------------------------- +// The main constructor +//------------------------------------------------------- + fClusterArray.Setup(par); + fClusterArray.SetClusterType("AliTPCcluster"); +} void AliTPCclusterer::FindPeaks(Int_t k,Int_t max, AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) { @@ -87,65 +93,37 @@ AliTPCcluster &c) { } //_____________________________________________________________________________ -void AliTPCclusterer::Digits2Clusters(const AliTPCParam *par, AliLoader *of, Int_t eventn) -{ +Int_t AliTPCclusterer::Digits2Clusters(TTree *dTree, TTree *cTree) { //----------------------------------------------------------------- // This is a simple cluster finder. //----------------------------------------------------------------- - AliRunLoader* rl = (AliRunLoader*)of->GetEventFolder()->FindObject(AliRunLoader::fgkRunLoaderName); - rl->GetEvent(eventn); - TDirectory *savedir=gDirectory; - - if (of->TreeR() == 0x0) of->MakeTree("R"); - - - if (of == 0x0) - { - cerr<<"AliTPC::Digits2Clusters(): output file not open !\n"; - return; - } - - const Int_t kMAXZ=par->GetMaxTBin()+2; - - - TTree *t = (TTree *)of->TreeD(); - - if (!t) { - cerr<<"Input tree with digits not found"<GetBranch("Segment"); + if (!branch) { + Error("Digits2Cluster","Can't get the branch !"); + return 1; } - AliSimDigits digarr, *dummy=&digarr; - t->GetBranch("Segment")->SetAddress(&dummy); - Stat_t nentries = t->GetEntries(); - - cout<<"Got "<Write(par->GetTitle()); + branch->SetAddress(&dummy); - AliTPCClustersArray carray; - carray.Setup(par); - carray.SetClusterType("AliTPCcluster"); - - TTree* treeR = of->TreeR(); - carray.MakeTree(treeR); - + fClusterArray.MakeTree(cTree); + AliTPCParam *par=(AliTPCParam *)fClusterArray.GetParam(); + const Int_t kMAXZ=par->GetMaxTBin()+2; Int_t nclusters=0; - for (Int_t n=0; nGetEntries(); + for (Int_t n=0; nGetEvent(n); + dTree->GetEvent(n); if (!par->AdjustSectorRow(digarr.GetID(),sec,row)) { - cerr<<"AliTPC warning: invalid segment ID ! "<GetPadRowRadii(sec,row); @@ -254,23 +232,16 @@ void AliTPCclusterer::Digits2Clusters(const AliTPCParam *par, AliLoader *of, Int clrow->InsertCluster(&c); ncl++; } } - carray.StoreRow(sec,row); - carray.ClearRow(sec,row); - - //cerr<<"sector, row, compressed digits, clusters: " - //<WriteRecPoints("OVERWRITE"); - - savedir->cd(); + Info("Digits2Cluster","Number of found clusters : %d",nclusters); -// delete t; //Thanks to Mariana Bondila + return 0; } diff --git a/TPC/AliTPCclusterer.h b/TPC/AliTPCclusterer.h index 0b12df79e0c..bfe05583e43 100644 --- a/TPC/AliTPCclusterer.h +++ b/TPC/AliTPCclusterer.h @@ -6,41 +6,44 @@ /* $Id$ */ //------------------------------------------------------- -// TPC clusterer +// The TPC cluster finder // // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch //------------------------------------------------------- -#include +#include +#include "AliTPCClustersArray.h" -#define kMAXCLUSTER 2500 - -class TFile; +class TTree; class AliTPCParam; class AliTPCcluster; -class AliLoader; -class AliTPCclusterer { +class AliTPCclusterer : public TObject { public: - static void Digits2Clusters(const AliTPCParam *par, AliLoader *of, Int_t eventn=1); + AliTPCclusterer(const AliTPCParam *par); + Int_t Digits2Clusters(TTree *dig, TTree *clu); private: - class AliBin { - public: - UShort_t GetQ() const {return fQ;} - UInt_t GetMask() const {return fMask;} - void SetQ(UShort_t q) {fQ=q;} - void SetMask(UInt_t m) {fMask=m;} - private: - UShort_t fQ; //signal - UInt_t fMask; //peak mask - }; + class AliBin { + public: + UShort_t GetQ() const {return fQ;} + UInt_t GetMask() const {return fMask;} + void SetQ(UShort_t q) {fQ=q;} + void SetMask(UInt_t m) {fMask=m;} + private: + UShort_t fQ; //signal + UInt_t fMask; //peak mask + }; private: - static Bool_t IsMaximum(Int_t k, Int_t max, const AliBin *bins); + static Bool_t IsMaximum(Int_t k, Int_t max, const AliBin *bins); static void FindPeaks(Int_t k,Int_t m,AliBin*b,Int_t*idx,UInt_t*msk,Int_t&n); - static void MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m); - static void MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m, + static void MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m); + static void MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m, AliTPCcluster &c); + + AliTPCClustersArray fClusterArray; //! container managing the clusters + + ClassDef(AliTPCclusterer,1) // the TPC cluster finder }; diff --git a/TPC/AliTPCpidESD.cxx b/TPC/AliTPCpidESD.cxx index c84c84da3d5..fa7098cdbde 100644 --- a/TPC/AliTPCpidESD.cxx +++ b/TPC/AliTPCpidESD.cxx @@ -68,7 +68,7 @@ Int_t AliTPCpidESD::MakePID(AliESD *event) Double_t bethe=Bethe(mom/mass); Double_t sigma=fRes*bethe; if (TMath::Abs(dedx-bethe) > fRange*sigma) { - p[j]=0.; + p[j]=TMath::Exp(-0.5*fRange*fRange); continue; } p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma)); diff --git a/TPC/AliTPCtest.C b/TPC/AliTPCtest.C index 5242de10456..3b09fea3f49 100644 --- a/TPC/AliTPCtest.C +++ b/TPC/AliTPCtest.C @@ -10,19 +10,18 @@ Int_t AliTPCtest(Int_t n = 5) { grun(n); - AliKalmanTrack::SetConvConst(1000/0.299792458/gAlice->Field()->SolenoidField()); +AliKalmanTrack::SetConvConst(1000/0.299792458/gAlice->Field()->SolenoidField()); Int_t ver=gAlice->GetDetector("TPC")->IsVersion(); - AliRunLoader* rl = gAlice->GetRunLoader(); - if (rl == 0x0) - { + AliRunLoader *rl = gAlice->GetRunLoader(); + if (rl == 0x0) { cerr<<"Can not get run loader from gAlice"< Proceeding witg digitization \n\n\n"; if ((ver!=1)&&(ver!=2)) { diff --git a/TPC/AliTPCtrack.cxx b/TPC/AliTPCtrack.cxx index 94f170f5431..5cc00a5908c 100644 --- a/TPC/AliTPCtrack.cxx +++ b/TPC/AliTPCtrack.cxx @@ -255,6 +255,15 @@ Double_t AliTPCtrack::GetPredictedChi2(const AliCluster *c) const return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det; } +Double_t AliTPCtrack::GetYat(Double_t xk) const { +//----------------------------------------------------------------- +// This function calculates the Y-coordinate of a track at the plane x=xk. +//----------------------------------------------------------------- + Double_t c1=fP4*fX - fP2, r1=TMath::Sqrt(1.- c1*c1); + Double_t c2=fP4*xk - fP2, r2=TMath::Sqrt(1.- c2*c2); + return fP0 + (xk-fX)*(c1+c2)/(r1+r2); +} + //_____________________________________________________________________________ Int_t AliTPCtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho) { //----------------------------------------------------------------- diff --git a/TPC/AliTPCtrack.h b/TPC/AliTPCtrack.h index 5764a472cc1..ee117edfdfd 100644 --- a/TPC/AliTPCtrack.h +++ b/TPC/AliTPCtrack.h @@ -45,6 +45,7 @@ public: Int_t Rotate(Double_t angle); void SetdEdx(Double_t dedx) {fdEdx=dedx;} + Double_t GetYat(Double_t x) const ; Double_t GetX() const {return fX;} Double_t GetAlpha() const {return fAlpha;} Double_t GetdEdx() const {return fdEdx;} diff --git a/TPC/AliTPCtracker.cxx b/TPC/AliTPCtracker.cxx index dee9237690b..5e26f495c42 100644 --- a/TPC/AliTPCtracker.cxx +++ b/TPC/AliTPCtracker.cxx @@ -13,76 +13,149 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* $Id$ */ +/* +$Log$ +Revision 1.35.2.3 2003/07/15 09:58:03 hristov +Corrected back-propagation (Yu.Belikov) + +Revision 1.35.2.2 2003/07/14 09:19:33 hristov +TOF included in the combined PID (Yu.Belikov) + +Revision 1.35.2.1 2003/07/11 10:53:01 hristov +Inward refit for TPC and TRD in the ESD schema (T.Kuhr) + +Revision 1.35 2003/05/23 10:08:51 hristov +SetLabel replaced by SetNumber (Yu.Belikov) + +Revision 1.34 2003/05/22 13:57:48 hristov +First implementation of ESD classes (Yu.Belikov) + +Revision 1.32 2003/04/10 10:36:54 hristov +Code for unified TPC/TRD tracking (S.Radomski) + +Revision 1.31 2003/03/19 17:14:11 hristov +Load/UnloadClusters added to the base class and the derived classes changed correspondingly. Possibility to give 2 input files for ITS and TPC tracks in PropagateBack. TRD tracker uses fEventN from the base class (T.Kuhr) + +Revision 1.30 2003/02/28 16:13:32 hristov +Typos corrected + +Revision 1.29 2003/02/28 15:18:16 hristov +Corrections suggested by J.Chudoba + +Revision 1.28 2003/02/27 16:15:52 hristov +Code for inward refitting (S.Radomski) + +Revision 1.27 2003/02/25 16:47:58 hristov +allow back propagation for more than 1 event (J.Chudoba) + +Revision 1.26 2003/02/19 08:49:46 hristov +Track time measurement (S.Radomski) + +Revision 1.25 2003/01/28 16:43:35 hristov +Additional protection: to be discussed with the Root team (M.Ivanov) + +Revision 1.24 2002/11/19 16:13:24 hristov +stdlib.h included to declare exit() on HP + +Revision 1.23 2002/11/19 11:50:08 hristov +Removing CONTAINERS (Yu.Belikov) + +Revision 1.19 2002/07/19 07:31:40 kowal2 +Improvement in tracking by J. Belikov + +Revision 1.18 2002/05/13 07:33:52 kowal2 +Added protection in Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const +in the case of defined region of interests. + +Revision 1.17 2002/03/18 17:59:13 kowal2 +Chnges in the pad geometry - 3 pad lengths introduced. + +Revision 1.16 2001/11/08 16:39:03 hristov +Additional protection (M.Masera) + +Revision 1.15 2001/11/08 16:36:33 hristov +Updated V2 stream of tracking (Yu.Belikov). The new long waited features are: 1) Possibility to pass the primary vertex position to the trackers (both for the TPC and the ITS) 2) Possibility to specify the number of tracking passes together with applying (or not applying) the vertex constraint (ITS only) 3) Possibility to make some use of partial PID provided by the TPC when doing tracking in the ITS (ITS only) 4) V0 reconstruction with a helix minimisation of the DCA. (new macros: AliV0FindVertices.C and AliV0Comparison.C) 4a) ( Consequence of the 4) ) All the efficiencies and resolutions are from now on calculated including *secondary*particles* too. (Don't be surprised by the drop in efficiency etc) + +Revision 1.14 2001/10/21 19:04:55 hristov +Several patches were done to adapt the barel reconstruction to the multi-event case. Some memory leaks were corrected. (Yu.Belikov) + +Revision 1.13 2001/07/23 12:01:30 hristov +Initialisation added + +Revision 1.12 2001/07/20 14:32:44 kowal2 +Processing of many events possible now + +Revision 1.11 2001/05/23 08:50:10 hristov +Weird inline removed + +Revision 1.10 2001/05/16 14:57:25 alibrary +New files for folders and Stack + +Revision 1.9 2001/05/11 07:16:56 hristov +Fix needed on Sun and Alpha + +Revision 1.8 2001/05/08 15:00:15 hristov +Corrections for tracking in arbitrary magnenetic field. Changes towards a concept of global Alice track. Back propagation of reconstructed tracks (Yu.Belikov) + +Revision 1.5 2000/12/20 07:51:59 kowal2 +Changes suggested by Alessandra and Paolo to avoid overlapped +data fields in encapsulated classes. + +Revision 1.4 2000/11/02 07:27:16 kowal2 +code corrections + +Revision 1.2 2000/06/30 12:07:50 kowal2 +Updated from the TPC-PreRelease branch + +Revision 1.1.2.1 2000/06/25 08:53:55 kowal2 +Splitted from AliTPCtracking + +*/ //------------------------------------------------------- // Implementation of the TPC tracker // // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch //------------------------------------------------------- - #include +#include #include #include -#include -#include -#include +#include + +#include "AliESD.h" #include "AliTPCtracker.h" #include "AliTPCcluster.h" #include "AliTPCParam.h" -#include "AliTPCClustersRow.h" -#include "AliTPCcluster.h" +#include "AliClusters.h" + +ClassImp(AliTPCtracker) //_____________________________________________________________________________ -AliTPCtracker::AliTPCtracker(const AliTPCParam *par, Int_t eventn, const char* evfoldname): -AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2), -fEvFolderName(evfoldname) +AliTPCtracker::AliTPCtracker(const AliTPCParam *par): +AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2) { //--------------------------------------------------------------------- // The main TPC tracker constructor //--------------------------------------------------------------------- - cout<<"fkNIS = "<GetEvent(fEventN); + fParam = (AliTPCParam*) par; + fSeeds=0; - AliLoader* tpcl = rl->GetLoader("TPCLoader"); - if (tpcl == 0x0) - { - Error("AliTPCtracker","Can not get TPC Laoder from Run Loader"); - return; - } + // [SR 17.03.2003] + + fBarrelFile = 0; + fBarrelTree = 0; + fBarrelArray = 0; + fBarrelTrack = 0; - if (tpcl->TreeR() == 0x0) tpcl->LoadRecPoints("read"); - TTree* treeR = tpcl->TreeR(); - if (treeR == 0x0) - { - cout<<"Error: Can not get TreeR\n"; - } - fClustersArray.Setup(par); - fClustersArray.SetClusterType("AliTPCcluster"); - fClustersArray.ConnectTree(treeR); - - fSeeds=0; } //_____________________________________________________________________________ @@ -96,6 +169,96 @@ AliTPCtracker::~AliTPCtracker() { fSeeds->Delete(); delete fSeeds; } + + // [SR, 01.04.2003] + if (fBarrelFile) { + fBarrelFile->Close(); + delete fBarrelFile; + } +} + +//_____________________________________________________________________________ + +void AliTPCtracker::SetBarrelTree(const char *mode) { + // + // Creates a tree for BarrelTracks + // mode = "back" or "refit" + // + // [SR, 01.04.2003] + // + + if (!IsStoringBarrel()) return; + + TDirectory *sav = gDirectory; + if (!fBarrelFile) fBarrelFile = new TFile("AliBarrelTracks.root", "UPDATE"); + + char buff[40]; + sprintf(buff, "BarrelTPC_%d_%s", GetEventNumber(), mode); + + fBarrelFile->cd(); + fBarrelTree = new TTree(buff, "Barrel TPC tracks"); + + if (!fBarrelArray) fBarrelArray = new TClonesArray("AliBarrelTrack", 4); + for(Int_t i=0; i<4; i++) new((*fBarrelArray)[i]) AliBarrelTrack(); + + fBarrelTree->Branch("tracks", &fBarrelArray); + + sav->cd(); +} +//_____________________________________________________________________________ + +void AliTPCtracker::StoreBarrelTrack(AliTPCtrack *ps, Int_t refPlane, Int_t isIn) { + // + // Stores Track at a given reference plane + // + // refPlane: 1-4 + // isIn: 1 - backward, 2 - refit + // + + if (!IsStoringBarrel()) return; + if (refPlane < 0 || refPlane > 4) return; + if (isIn > 2) return; + + static Int_t nClusters; + static Int_t nWrong; + static Double_t chi2; + static Int_t index; + + Int_t newClusters, newWrong; + Double_t newChi2; + + if ( (refPlane == 1 && isIn == kTrackBack) || + (refPlane == 4 && isIn == kTrackRefit) ) { + + fBarrelArray->Clear(); + nClusters = nWrong = 0; + chi2 = 0.0; + index = 0; + } + + // propagate + Double_t refX = 0; + if (refPlane == 1) refX = fParam->GetInnerRadiusLow(); + if (refPlane == 2) refX = fParam->GetInnerRadiusUp(); + if (refPlane == 3) refX = fParam->GetOuterRadiusLow(); + if (refPlane == 4) refX = fParam->GetOuterRadiusUp(); + + ps->PropagateTo(refX); + + fBarrelTrack = (AliBarrelTrack*)(*fBarrelArray)[index++]; + ps->GetBarrelTrack(fBarrelTrack); + + newClusters = ps->GetNumberOfClusters() - nClusters; + newWrong = ps->GetNWrong() - nWrong; + newChi2 = ps->GetChi2() - chi2; + + nClusters = ps->GetNumberOfClusters(); + nWrong = ps->GetNWrong(); + chi2 = ps->GetChi2(); + + fBarrelTrack->SetNClusters(newClusters, newChi2); + fBarrelTrack->SetNWrongClusters(newWrong); + fBarrelTrack->SetRefPlane(refPlane, isIn); } //_____________________________________________________________________________ @@ -191,98 +354,75 @@ Double_t f3(Double_t x1,Double_t y1, } //_____________________________________________________________________________ -void AliTPCtracker::LoadOuterSectors() { +Int_t AliTPCtracker::LoadClusters(TTree *cTree) { //----------------------------------------------------------------- - // This function fills outer TPC sectors with clusters. + // This function loads TPC clusters. //----------------------------------------------------------------- - UInt_t index; - TTree* tree = fClustersArray.GetTree(); - if (tree == 0x0) - { - Error("LoadOuterSectors","Can not get tree from fClustersArray"); - return; - } - Int_t j=Int_t(tree->GetEntries()); - cout<<"fClustersArray.GetTree()->GetEntries() = "<AdjustSectorRow(s->GetID(),sec,row); - if (secGetArray()->GetEntriesFast(); - while (ncl--) { - AliTPCcluster *c=(AliTPCcluster*)(*clrow)[ncl]; - index=(((sec<<8)+row)<<16)+ncl; - fOuterSec[(sec-fkNIS*2)%fkNOS][row].InsertCluster(c,index); - } + TBranch *branch=cTree->GetBranch("Segment"); + if (!branch) { + Error("LoadClusters","Can't get the branch !"); + return 1; } - fN=fkNOS; - fSectors=fOuterSec; -} + AliClusters carray, *addr=&carray; + carray.SetClass("AliTPCcluster"); + carray.SetArray(0); + branch->SetAddress(&addr); -//_____________________________________________________________________________ -void AliTPCtracker::UnloadOuterSectors() { - //----------------------------------------------------------------- - // This function clears outer TPC sectors. - //----------------------------------------------------------------- - Int_t nup=fOuterSec->GetNRows(); - for (Int_t i=0; iGetEntries(); - fN=0; - fSectors=0; -} + for (Int_t i=0; iGetEvent(i); -//_____________________________________________________________________________ -void AliTPCtracker::LoadInnerSectors() { - //----------------------------------------------------------------- - // This function fills inner TPC sectors with clusters. - //----------------------------------------------------------------- - UInt_t index; - Int_t j=Int_t(fClustersArray.GetTree()->GetEntries()); - for (Int_t i=0; iAdjustSectorRow(s->GetID(),sec,row); - if (sec>=fkNIS*2) continue; - AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row); - Int_t ncl=clrow->GetArray()->GetEntriesFast(); - while (ncl--) { - AliTPCcluster *c=(AliTPCcluster*)(*clrow)[ncl]; - index=(((sec<<8)+row)<<16)+ncl; - fInnerSec[sec%fkNIS][row].InsertCluster(c,index); + Int_t ncl=carray.GetArray()->GetEntriesFast(); + + Int_t nir=fInnerSec->GetNRows(), nor=fOuterSec->GetNRows(); + Int_t id=carray.GetID(); + if ((id<0) || (id>2*(fkNIS*nir + fkNOS*nor))) { + Error("LoadClusters","Wrong index !"); + exit(1); } + Int_t outindex = 2*fkNIS*nir; + if (idClear(); } - fN=fkNIS; - fSectors=fInnerSec; + return 0; } //_____________________________________________________________________________ -void AliTPCtracker::UnloadInnerSectors() { +void AliTPCtracker::UnloadClusters() { //----------------------------------------------------------------- - // This function clears inner TPC sectors. + // This function unloads TPC clusters. //----------------------------------------------------------------- - Int_t nlow=fInnerSec->GetNRows(); - for (Int_t i=0; iGetNRows(); + for (Int_t n=0; nGetNRows(); + for (Int_t n=0; nkMaxROAD) { if (t.GetNumberOfClusters()>4) - cerr<GetY() > y+road) break; - if (c->IsUsed()) continue; - if ((c->GetZ()-z)*(c->GetZ()-z) > 16.*(t.GetSigmaZ2()+sz2)) continue; - Double_t chi2=t.GetPredictedChi2(c); - if (chi2 > maxchi2) continue; - maxchi2=chi2; - cl=c; + AliTPCcluster *c=(AliTPCcluster*)(krow[i]); + if (c->GetY() > y+road) break; + if (c->IsUsed()) continue; + if ((c->GetZ()-z)*(c->GetZ()-z) > 16.*(t.GetSigmaZ2()+sz2)) continue; + Double_t chi2=t.GetPredictedChi2(c); + if (chi2 > maxchi2) continue; + maxchi2=chi2; + cl=c; index=krow.GetIndex(i); } } @@ -356,6 +495,57 @@ Int_t AliTPCtracker::FollowProlongation(AliTPCseed& t, Int_t rf) { return 1; } +//_____________________________________________________________________________ + +Int_t AliTPCtracker::FollowRefitInward(AliTPCseed *seed, AliTPCtrack *track) { + // + // This function propagates seed inward TPC using old clusters + // from the track. + // + // Sylwester Radomski, GSI + // 26.02.2003 + // + + // loop over rows + + Int_t nRows = fSectors->GetNRows(); + for (Int_t iRow = nRows; iRow >= 0; iRow--) { + + Double_t x = fSectors->GetX(iRow); + if (!seed->PropagateTo(x)) return 0; + + // try to find an assigned cluster in this row + + AliTPCcluster* cluster = NULL; + Int_t idx = -1; + Int_t sec = -1; + for (Int_t iCluster = 0; iCluster < track->GetNumberOfClusters(); iCluster++) { + idx = track->GetClusterIndex(iCluster); + sec = (idx&0xff000000)>>24; + Int_t row = (idx&0x00ff0000)>>16; + if (((fSectors == fInnerSec) && (sec >= fkNIS)) || + ((fSectors == fOuterSec) && (sec < fkNIS))) continue; + if (row == iRow) { + cluster = (AliTPCcluster*) GetCluster(idx); + break; + } + } + + // update the track seed with the found cluster + + if (cluster) { + Double_t dAlpha = fParam->GetAngle(sec) - seed->GetAlpha(); + if (TMath::Abs(dAlpha) > 0.0001) { + if (!seed->Rotate(dAlpha)) return 0; + if (!seed->PropagateTo(x)) return 0; + } + + seed->Update(cluster, seed->GetPredictedChi2(cluster), idx); + } + } + + return 1; +} //_____________________________________________________________________________ Int_t AliTPCtracker::FollowBackProlongation @@ -369,21 +559,20 @@ Int_t AliTPCtracker::FollowBackProlongation Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN; Int_t idx=-1, sec=-1, row=-1; - Int_t nc=seed.GetLabel(); //index of the cluster to start with + Int_t nc=seed.GetNumber(); + if (nc--) { idx=track.GetClusterIndex(nc); sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16; } - if (fSectors==fInnerSec) { if (sec >= 2*fkNIS) row=-1; } - else { if (sec < 2*fkNIS) row=-1; } + if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; } + else { if (sec < fkNIS) row=-1; } Int_t nr=fSectors->GetNRows(); for (Int_t i=0; iGetX(i), ymax=fSectors->GetMaxY(i); - - if (!seed.PropagateTo(x)) return 0; - - Double_t y=seed.GetY(); + Double_t y=seed.GetYat(x); + if (y > ymax) { s = (s+1) % fN; if (!seed.Rotate(fSectors->GetAlpha())) return 0; @@ -392,6 +581,8 @@ Int_t AliTPCtracker::FollowBackProlongation if (!seed.Rotate(-fSectors->GetAlpha())) return 0; } + if (!seed.PropagateTo(x)) return 0; + AliTPCcluster *cl=0; Int_t index=0; Double_t maxchi2=kMaxCHI2; @@ -400,12 +591,10 @@ Int_t AliTPCtracker::FollowBackProlongation Double_t sz2=SigmaZ2(seed.GetX(),seed.GetTgl()); Double_t road=4.*sqrt(seed.GetSigmaY2() + sy2), z=seed.GetZ(); if (road>kMaxROAD) { - cerr<>24; row=(idx&0x00ff0000)>>16; } - if (fSectors==fInnerSec) { if (sec >= 2*fkNIS) row=-1; } - else { if (sec < 2*fkNIS) row=-1; } + if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; } + else { if (sec < fkNIS) row=-1; } } if (!cl) { @@ -429,14 +618,14 @@ Int_t AliTPCtracker::FollowBackProlongation if (accepted>27) if (krow) { for (Int_t i=krow.Find(y-road); iGetY() > y+road) break; - if (c->IsUsed()) continue; - if ((c->GetZ()-z)*(c->GetZ()-z)>16.*(seed.GetSigmaZ2()+sz2)) continue; - Double_t chi2=seed.GetPredictedChi2(c); - if (chi2 > maxchi2) continue; - maxchi2=chi2; - cl=c; + AliTPCcluster *c=(AliTPCcluster*)(krow[i]); + if (c->GetY() > y+road) break; + if (c->IsUsed()) continue; + if ((c->GetZ()-z)*(c->GetZ()-z)>16.*(seed.GetSigmaZ2()+sz2)) continue; + Double_t chi2=seed.GetPredictedChi2(c); + if (chi2 > maxchi2) continue; + maxchi2=chi2; + cl=c; index=krow.GetIndex(i); } } @@ -451,7 +640,7 @@ Int_t AliTPCtracker::FollowBackProlongation } - seed.SetLabel(nc); + seed.SetNumber(nc); return 1; } @@ -461,81 +650,80 @@ void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) { //----------------------------------------------------------------- // This function creates track seeds. //----------------------------------------------------------------- - cout<<" Making Seeds i1="<GetAlpha(), shift=fOuterSec->GetAlphaShift(); + Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift(); Double_t cs=cos(alpha), sn=sin(alpha); - Double_t x1 =fOuterSec->GetX(i1); - Double_t xx2=fOuterSec->GetX(i2); + Double_t x1 =fSectors->GetX(i1); + Double_t xx2=fSectors->GetX(i2); - for (Int_t ns=0; nsGetY(), z1=kr1[is]->GetZ(); for (Int_t js=0; js < nl+nm+nu; js++) { - const AliTPCcluster *kcl; + const AliTPCcluster *kcl; Double_t x2, y2, z2; Double_t x3=GetX(), y3=GetY(), z3=GetZ(); - if (jsGetY(); z2=kcl->GetZ(); x2= xx2*cs+y2*sn; y2=-xx2*sn+y2*cs; - } else - if (jsGetY(); z2=kcl->GetZ(); - } else { - const AliTPCRow& kr2=fOuterSec[(ns+1)%fkNOS][i2]; - kcl=kr2[js-nl-nm]; + } else { + const AliTPCRow& kr2=fSectors[(ns+1)%fN][i2]; + kcl=kr2[js-nl-nm]; y2=kcl->GetY(); z2=kcl->GetZ(); x2=xx2*cs-y2*sn; y2=xx2*sn+y2*cs; - } + } Double_t zz=z1 - (z1-z3)/(x1-x3)*(x1-x2); if (TMath::Abs(zz-z2)>5.) continue; Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1); - if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;} - - x[0]=y1; - x[1]=z1; - x[4]=f1(x1,y1,x2,y2,x3,y3); - if (TMath::Abs(x[4]) >= 0.0066) continue; - x[2]=f2(x1,y1,x2,y2,x3,y3); - //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue; - x[3]=f3(x1,y1,x2,y2,z1,z2); - if (TMath::Abs(x[3]) > 1.2) continue; - Double_t a=asin(x[2]); - Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2])); - if (TMath::Abs(zv-z3)>10.) continue; + if (d==0.) { + Warning("MakeSeeds","Straight seed !"); + continue; + } + x[0]=y1; + x[1]=z1; + x[4]=f1(x1,y1,x2,y2,x3,y3); + if (TMath::Abs(x[4]) >= 0.0066) continue; + x[2]=f2(x1,y1,x2,y2,x3,y3); + //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue; + x[3]=f3(x1,y1,x2,y2,z1,z2); + if (TMath::Abs(x[3]) > 1.2) continue; + Double_t a=asin(x[2]); + Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2])); + if (TMath::Abs(zv-z3)>10.) continue; Double_t sy1=kr1[is]->GetSigmaY2(), sz1=kr1[is]->GetSigmaZ2(); Double_t sy2=kcl->GetSigmaY2(), sz2=kcl->GetSigmaZ2(); //Double_t sy3=400*3./12., sy=0.1, sz=0.1; Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1; - Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy; - Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy; - Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy; - Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy; - Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy; - Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy; - Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy; - Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz; - Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy; - Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz; + Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy; + Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy; + Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy; + Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy; + Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy; + Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy; + Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy; + Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz; + Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy; + Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz; c[0]=sy1; c[1]=0.; c[2]=sz1; @@ -547,17 +735,13 @@ void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) { c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43; UInt_t index=kr1.GetIndex(is); - AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift); - Float_t l=fOuterSec->GetPadPitchWidth(); + AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift); + Float_t l=fSectors->GetPadPitchWidth(); track->SetSampledEdx(kr1[is]->GetQ()/l,0); Int_t rc=FollowProlongation(*track, i2); if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/2) delete track; - else - { - fSeeds->AddLast(track); - cout<<"Adding seed "<GetEntries()<<"\r"; - } + else fSeeds->AddLast(track); } } } @@ -572,15 +756,14 @@ Int_t AliTPCtracker::ReadSeeds(const TFile *inp) { TFile *in=(TFile*)inp; if (!in->IsOpen()) { - cerr<<"AliTPCtracker::ReadSeeds(): input file is not open !\n"; + Error("ReadSeeds","Input file has not been open !"); return 1; } in->cd(); TTree *seedTree=(TTree*)in->Get("Seeds"); if (!seedTree) { - cerr<<"AliTPCtracker::ReadSeeds(): "; - cerr<<"can't get a tree with track seeds !\n"; + Error("ReadSeeds","Can't get a tree with track seeds !"); return 2; } AliTPCtrack *seed=new AliTPCtrack; @@ -604,189 +787,451 @@ Int_t AliTPCtracker::ReadSeeds(const TFile *inp) { } //_____________________________________________________________________________ -Int_t AliTPCtracker::Clusters2Tracks() -{ +Int_t AliTPCtracker::Clusters2Tracks(AliESD *event) { //----------------------------------------------------------------- // This is a track finder. + // The clusters must be already loaded ! //----------------------------------------------------------------- - Int_t retval = 0; - - AliRunLoader* rl = AliRunLoader::GetRunLoader(fEvFolderName); - if (rl == 0x0) - { - Error("Clusters2Tracks","Can not get RL from specified folder %s",fEvFolderName.Data()); - return 1; - } - - retval = rl->GetEvent(fEventN); - if (retval) - { - Error("Clusters2Tracks","Error while getting event %d",fEventN); - return 1; - } - - AliLoader* tpcl = rl->GetLoader("TPCLoader"); - if (tpcl == 0x0) - { - Error("Clusters2Tracks","Can not get TPC Laoder from Run Loader"); - return 1; - } - - if ( tpcl->TreeR() == 0x0) - { - retval = tpcl->LoadRecPoints("READ"); - if (retval) - { - Error("Clusters2Tracks","Error while loading Reconstructed Points"); - return 1; - } - } - - if (tpcl->TreeT() == 0x0) tpcl->MakeTree("T"); - - TTree &tracktree = *(tpcl->TreeT()); - - TBranch* br= tracktree.GetBranch("tracks"); - if (br) - { - Error("Clusters2Tracks","Branch \"tracks\" already exists in TreeT for TPC"); - return 1; - } - - AliTPCtrack *iotrack=0; - tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,0); - - LoadOuterSectors(); //find track seeds Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows(); Int_t nrows=nlow+nup; if (fSeeds==0) { Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap); + fSectors=fOuterSec; fN=fkNOS; + fSeeds=new TObjArray(15000); MakeSeeds(nup-1, nup-1-gap); MakeSeeds(nup-1-shift, nup-1-shift-gap); - } + } fSeeds->Sort(); - //tracking in outer sectors Int_t nseed=fSeeds->GetEntriesFast(); - cout<<"nseed="<UncheckedAt(i), &t=*pt; - if (FollowProlongation(t)) { - UseClusters(&t); + if (!FollowProlongation(t)) { + delete fSeeds->RemoveAt(i); continue; } + + //tracking in the inner sectors + fSectors=fInnerSec; fN=fkNIS; + + Double_t alpha=t.GetAlpha() - fInnerSec->GetAlphaShift(); + if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi(); + if (alpha < 0. ) alpha += 2.*TMath::Pi(); + Int_t ns=Int_t(alpha/fInnerSec->GetAlpha())%fkNIS; + + alpha=ns*fInnerSec->GetAlpha()+fInnerSec->GetAlphaShift()-t.GetAlpha(); + + if (t.Rotate(alpha)) { + if (FollowProlongation(t)) { + if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) { + t.CookdEdx(); + CookLabel(pt,0.1); //For comparison only + pt->PropagateTo(fParam->GetInnerRadiusLow()); + AliESDtrack iotrack; + iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin); + + event->AddTrack(&iotrack); + + UseClusters(&t); + } + } + } delete fSeeds->RemoveAt(i); - } - //UnloadOuterSectors(); + } + + Info("Clusters2Tracks","Number of found tracks : %d", + event->GetNumberOfTracks()); + + fSeeds->Clear(); delete fSeeds; fSeeds=0; + + return 0; +} + +//_____________________________________________________________________________ +Int_t AliTPCtracker::Clusters2Tracks(TTree *cTree, TTree *tTree) { + //----------------------------------------------------------------- + // This is a track finder. + //----------------------------------------------------------------- + if (LoadClusters(cTree)) { + Error("Clusters2Tracks","Problem with loading the clusters..."); + return 1; + } + + AliTPCtrack *iotrack=0; + TBranch *branch=tTree->GetBranch("tracks"); + if (!branch) tTree->Branch("tracks","AliTPCtrack",&iotrack,32000,3); + else branch->SetAddress(&iotrack); + + //find track seeds + Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows(); + Int_t nrows=nlow+nup; + if (fSeeds==0) { + Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap); + fSectors=fOuterSec; fN=fkNOS; + fSeeds=new TObjArray(15000); + MakeSeeds(nup-1, nup-1-gap); + MakeSeeds(nup-1-shift, nup-1-shift-gap); + } + fSeeds->Sort(); - //tracking in inner sectors - LoadInnerSectors(); Int_t found=0; - for (i=0; iGetEntriesFast(); + + for (Int_t i=0; iUncheckedAt(i), &t=*pt; - if (!pt) continue; - Int_t nc=t.GetNumberOfClusters(); + if (!FollowProlongation(t)) { + delete fSeeds->RemoveAt(i); + continue; + } + + //tracking in the inner sectors + fSectors=fInnerSec; fN=fkNIS; Double_t alpha=t.GetAlpha() - fInnerSec->GetAlphaShift(); if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi(); if (alpha < 0. ) alpha += 2.*TMath::Pi(); Int_t ns=Int_t(alpha/fInnerSec->GetAlpha())%fkNIS; - alpha=ns*fInnerSec->GetAlpha() + fInnerSec->GetAlphaShift() - t.GetAlpha(); + alpha=ns*fInnerSec->GetAlpha()+fInnerSec->GetAlphaShift()-t.GetAlpha(); if (t.Rotate(alpha)) { - if (FollowProlongation(t)) { - if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) { - t.CookdEdx(); - CookLabel(pt,0.1); //For comparison only - iotrack=pt; - tracktree.Fill(); - UseClusters(&t,nc); - cerr<= Int_t(0.4*nrows)) { + t.CookdEdx(); + CookLabel(pt,0.1); //For comparison only + pt->PropagateTo(fParam->GetInnerRadiusLow()); + iotrack=pt; + tTree->Fill(); + UseClusters(&t); + found++; + } + } } delete fSeeds->RemoveAt(i); - } - UnloadInnerSectors(); - UnloadOuterSectors(); + } + + Info("Clusters2Tracks","Number of found tracks : %d",found); + + UnloadClusters(); + + fSeeds->Clear(); delete fSeeds; fSeeds=0; + + return 0; +} + +//_____________________________________________________________________________ +Int_t AliTPCtracker::RefitInward(AliESD* event) { + // + // The function propagates tracks throught TPC inward + // using already associated clusters. + // The clusters must be already loaded ! + // + + Int_t nTracks = event->GetNumberOfTracks(); + Int_t nRefited = 0; + + for (Int_t i = 0; i < nTracks; i++) { + AliESDtrack* track = event->GetTrack(i); + ULong_t status = track->GetStatus(); + + if ( (status & AliESDtrack::kTPCout ) == 0 ) continue; + if ( (status & AliESDtrack::kTPCrefit) != 0 ) continue; + + AliTPCtrack* tpcTrack = new AliTPCtrack(*track); + AliTPCseed* seed = new AliTPCseed(*tpcTrack, tpcTrack->GetAlpha()); - tpcl->WriteTracks("OVERWRITE"); + fSectors = fOuterSec; + + Int_t res = FollowRefitInward(seed, tpcTrack); + UseClusters(seed); + Int_t nc = seed->GetNumberOfClusters(); + + fSectors = fInnerSec; + + res = FollowRefitInward(seed, tpcTrack); + UseClusters(seed, nc); + + if (res) { + seed->PropagateTo(fParam->GetInnerRadiusLow()); + seed->SetLabel(tpcTrack->GetLabel()); + seed->SetdEdx(tpcTrack->GetdEdx()); + track->UpdateTrackParams(seed, AliESDtrack::kTPCrefit); + nRefited++; + } + + delete seed; + delete tpcTrack; + } + + Info("RefitInward","Number of refitted tracks : %d",nRefited); + + return 0; +} + +//_____________________________________________________________________________ +Int_t AliTPCtracker::RefitInward(TTree *in, TTree *out) { + // + // The function propagates tracks throught TPC inward + // using already associated clusters. + // + Error("RefitInward","This method is not converted to NewIO yet\n"); + return 1; + /* + if (!inSeeds->IsOpen()) { + cout << "Input File with seeds not open !\n" << endl; + return 1; + } + + if (!inTracks->IsOpen()) { + cout << "Input File not open !\n" << endl; + return 2; + } - cerr<<"Number of found tracks : "<IsOpen()) { + cout << "Output File not open !\n" << endl; + return 3; + } + + TDirectory *savedir = gDirectory; + LoadClusters(); + + // Connect to input seeds tree + inSeeds->cd(); + char tname[100]; + sprintf(tname, "seedTRDtoTPC_%d", GetEventNumber()); + TTree *seedsTree = (TTree*)inSeeds->Get(tname); + TBranch *inSeedBranch = seedsTree->GetBranch("tracks"); + AliTPCtrack *inSeedTrack = 0; + inSeedBranch->SetAddress(&inSeedTrack); + + Int_t nSeeds = (Int_t)seedsTree->GetEntries(); + + // Connect to input tree + inTracks->cd(); + sprintf(tname,"TreeT_TPC_%d",GetEventNumber()); +// sprintf(tname,"seedsTPCtoTRD_%d",GetEventNumber()); + TTree *inputTree = (TTree*)inTracks->Get(tname); + TBranch *inBranch = inputTree->GetBranch("tracks"); + AliTPCtrack *inTrack = 0; + inBranch->SetAddress(&inTrack); + + Int_t nTracks = (Int_t)inputTree->GetEntries(); + + // Create output tree + outTracks->cd(); + AliTPCtrack *outTrack = new AliTPCtrack(); + sprintf(tname, "tracksTPC_%d", GetEventNumber()); + TTree *outputTree = new TTree(tname,"Refited TPC tracks"); + outputTree->Branch("tracks", "AliTPCtrack", &outTrack, 32000, 3); + + Int_t nRefited = 0; + + // create table of track labels + Int_t* inLab = new Int_t[nTracks]; + for (Int_t i = 0; i < nTracks; i++) inLab[i] = -1; + + // [SR, 01.04.2003] - Barrel tracks + if (IsStoringBarrel()) SetBarrelTree("refit"); + + for(Int_t t=0; t < nSeeds; t++) { + + seedsTree->GetEntry(t); + // find TPC track for seed + Int_t lab = TMath::Abs(inSeedTrack->GetLabel()); + for(Int_t i=0; i < nTracks; i++) { + if (inLab[i] == lab) { + inputTree->GetEntry(i); + break; + } else if (inLab[i] < 0) { + inputTree->GetEntry(i); + inLab[i] = TMath::Abs(inTrack->GetLabel()); + if (inLab[i] == lab) break; + } + } + if (TMath::Abs(inTrack->GetLabel()) != lab) continue; + AliTPCseed *seed = new AliTPCseed(*inSeedTrack, inTrack->GetAlpha()); + if (IsStoringBarrel()) StoreBarrelTrack(seed, 4, 2); + + fSectors = fOuterSec; + + Int_t res = FollowRefitInward(seed, inTrack); + UseClusters(seed); + Int_t nc = seed->GetNumberOfClusters(); + + if (IsStoringBarrel()) StoreBarrelTrack(seed, 3, 2); + if (IsStoringBarrel()) StoreBarrelTrack(seed, 2, 2); + + fSectors = fInnerSec; + + res = FollowRefitInward(seed, inTrack); + UseClusters(seed, nc); + + if (IsStoringBarrel()) StoreBarrelTrack(seed, 1, 2); + + if (res) { + seed->PropagateTo(fParam->GetInnerRadiusLow()); + outTrack = (AliTPCtrack*)seed; + outTrack->SetLabel(inTrack->GetLabel()); + outTrack->SetdEdx(inTrack->GetdEdx()); + outputTree->Fill(); + nRefited++; + } + + if (IsStoringBarrel()) fBarrelTree->Fill(); + delete seed; + } + + cout << "Refitted " << nRefited << " tracks." << endl; + + outTracks->cd(); + outputTree->Write(); + + delete[] inLab; + + // [SR, 01.04.2003] + if (IsStoringBarrel()) { + fBarrelFile->cd(); + fBarrelTree->Write(); + fBarrelFile->Flush(); + + delete fBarrelArray; + delete fBarrelTree; + } + + if (inputTree) delete inputTree; + if (outputTree) delete outputTree; + + savedir->cd(); + UnloadClusters(); + + return 0; + */ +} + +Int_t AliTPCtracker::PropagateBack(AliESD *event) { + //----------------------------------------------------------------- + // This function propagates tracks back through the TPC. + // The clusters must be already loaded ! + //----------------------------------------------------------------- + Int_t nentr=event->GetNumberOfTracks(); + Info("PropagateBack", "Number of ESD tracks: %d\n", nentr); + + Int_t ntrk=0; + for (Int_t i=0; iGetTrack(i); + ULong_t status=esd->GetStatus(); + + if ( (status & AliESDtrack::kTPCin ) == 0 ) continue; + if ( (status & AliESDtrack::kTPCout) != 0 ) continue; + + const AliTPCtrack t(*esd); + AliTPCseed s(t,t.GetAlpha()); + + if ( (status & AliESDtrack::kITSout) == 0 ) s.ResetCovariance(); + + s.ResetNWrong(); + s.ResetNRotation(); + + Int_t nc=t.GetNumberOfClusters(); + s.SetNumber(nc); //set number of the cluster to start with + + //inner sectors + fSectors=fInnerSec; fN=fkNIS; + + Double_t alpha=s.GetAlpha() - fSectors->GetAlphaShift(); + if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi(); + if (alpha < 0. ) alpha += 2.*TMath::Pi(); + Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN; + alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift(); + alpha-=s.GetAlpha(); + + if (!s.Rotate(alpha)) continue; + if (!FollowBackProlongation(s,t)) continue; + + UseClusters(&s); + + //outer sectors + fSectors=fOuterSec; fN=fkNOS; + + nc=s.GetNumberOfClusters(); + + alpha=s.GetAlpha() - fSectors->GetAlphaShift(); + if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi(); + if (alpha < 0. ) alpha += 2.*TMath::Pi(); + ns=Int_t(alpha/fSectors->GetAlpha())%fN; + + alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift(); + alpha-=s.GetAlpha(); + + if (!s.Rotate(alpha)) continue; + if (!FollowBackProlongation(s,t)) continue; + { + Int_t nrows=fOuterSec->GetNRows()+fInnerSec->GetNRows(); + if (s.GetNumberOfClusters() < Int_t(0.4*nrows)) continue; + } + s.PropagateTo(fParam->GetOuterRadiusUp()); + s.CookdEdx(); + CookLabel(&s,0.1); //For comparison only + UseClusters(&s,nc); + esd->UpdateTrackParams(&s,AliESDtrack::kTPCout); + ntrk++; + } + Info("PropagateBack","Number of back propagated tracks: %d",ntrk); return 0; } //_____________________________________________________________________________ -Int_t AliTPCtracker::PropagateBack() - { +Int_t AliTPCtracker::PropagateBack(TTree *in, TTree *out) { //----------------------------------------------------------------- // This function propagates tracks back through the TPC. //----------------------------------------------------------------- - - cout<<"This method is not converted to NewIO yet\n"; + Error("PropagateBack","This method is not converted to NewIO yet\n"); return 1; - + /* fSeeds=new TObjArray(15000); - Int_t retval = 0; - AliRunLoader* rl = AliRunLoader::GetRunLoader(fEvFolderName); - if (rl == 0x0) - { - Error("Clusters2Tracks","Can not get RL from specified folder %s",fEvFolderName.Data()); - return 1; - } - - retval = rl->GetEvent(fEventN); - if (retval) - { - Error("Clusters2Tracks","Error while getting event %d",fEventN); - return 1; - } - - AliLoader* tpcl = rl->GetLoader("TPCLoader"); - if (tpcl == 0x0) - { - Error("Clusters2Tracks","Can not get TPC Laoder from Run Loader"); - return 1; - } + TFile *in=(TFile*)inp; + TFile *in2=(TFile*)inp2; + TDirectory *savedir=gDirectory; - AliLoader* itsl = rl->GetLoader("ITSLoader"); - if (tpcl == 0x0) - { - Error("Clusters2Tracks","Can not get ITS Laoder from Run Loader"); + if (!in->IsOpen()) { + cerr<<"AliTPCtracker::PropagateBack(): "; + cerr<<"file with TPC (or back propagated ITS) tracks is not open !\n"; return 1; - } - - if (itsl->TreeT() == 0x0) itsl->LoadTracks(); + } - TTree *bckTree=itsl->TreeT(); - if (!bckTree) { - cerr<<"AliTPCtracker::PropagateBack() "; - cerr<<"can't get a tree with back propagated ITS tracks !\n"; - return 3; + if (!out->IsOpen()) { + cerr<<"AliTPCtracker::PropagateBack(): "; + cerr<<"file for back propagated TPC tracks is not open !\n"; + return 2; } - - AliTPCtrack *bckTrack=new AliTPCtrack; - bckTree->SetBranchAddress("tracks",&bckTrack); + LoadClusters(); - TFile* in = 0x0; - TFile* out = 0x0; - cout<<"And NOW there will be a segmentation violation!!!!\n"; - bckTree=(TTree*)in->Get("TreeT_ITSb_0"); + in->cd(); + char tName[100]; + sprintf(tName,"TreeT_ITSb_%d",GetEventNumber()); + TTree *bckTree=(TTree*)in->Get(tName); + if (!bckTree && inp2) bckTree=(TTree*)in2->Get(tName); if (!bckTree) { cerr<<"AliTPCtracker::PropagateBack() "; cerr<<"can't get a tree with back propagated ITS tracks !\n"; - return 3; + //return 3; } + AliTPCtrack *bckTrack=new AliTPCtrack; + if (bckTree) bckTree->SetBranchAddress("tracks",&bckTrack); - - TTree *tpcTree=(TTree*)in->Get("TreeT_TPC_0"); + sprintf(tName,"TreeT_TPC_%d",GetEventNumber()); + TTree *tpcTree=(TTree*)in->Get(tName); if (!tpcTree) { cerr<<"AliTPCtracker::PropagateBack() "; cerr<<"can't get a tree with TPC tracks !\n"; @@ -795,61 +1240,128 @@ Int_t AliTPCtracker::PropagateBack() AliTPCtrack *tpcTrack=new AliTPCtrack; tpcTree->SetBranchAddress("tracks",&tpcTrack); -//*** Prepare an array of tracks to be back propagated + // [SR, 01.04.2003] - Barrel tracks + if (IsStoringBarrel()) SetBarrelTree("back"); + + // Prepare an array of tracks to be back propagated Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows(); Int_t nrows=nlow+nup; + // + // Match ITS tracks with old TPC tracks + // the tracks do not have to be sorted [SR, GSI, 18.02.2003] + // + // the algorithm is linear and uses LUT for sorting + // cost of the algorithm: 2 * number of tracks + // + TObjArray tracks(15000); - Int_t i=0,j=0; - Int_t tpcN=(Int_t)tpcTree->GetEntries(); - Int_t bckN=(Int_t)bckTree->GetEntries(); - if (jGetEvent(j++); - for (i=0; iGetEvent(i); - Double_t alpha=tpcTrack->GetAlpha(); - if (jGetLabel())==TMath::Abs(bckTrack->GetLabel())) { - if (!bckTrack->Rotate(alpha-bckTrack->GetAlpha())) continue; - fSeeds->AddLast(new AliTPCseed(*bckTrack,bckTrack->GetAlpha())); - bckTree->GetEvent(j++); - } else { - tpcTrack->ResetCovariance(); - fSeeds->AddLast(new AliTPCseed(*tpcTrack,alpha)); - } - tracks.AddLast(new AliTPCtrack(*tpcTrack)); + Int_t i=0; + Int_t tpcN= (Int_t)tpcTree->GetEntries(); + Int_t bckN= (bckTree)? (Int_t)bckTree->GetEntries() : 0; + + // look up table + const Int_t nLab = 200000; // limit on number of primaries (arbitrary) + Int_t lut[nLab]; + for(Int_t i=0; iGetEvent(i); + Int_t lab = TMath::Abs(bckTrack->GetLabel()); + if (lab < nLab) lut[lab] = i; + else { + cerr << "AliTPCtracker: The size of the LUT is too small\n"; + } + } + } + + for (Int_t i=0; iGetEvent(i); + Double_t alpha=tpcTrack->GetAlpha(); + + if (!bckTree) { + + // No ITS - use TPC track only + + tpcTrack->ResetCovariance(); + AliTPCseed * seed = new AliTPCseed(*tpcTrack, tpcTrack->GetAlpha()); + + fSeeds->AddLast(seed); + tracks.AddLast(new AliTPCtrack(*tpcTrack)); + + } else { + + // with ITS + // discard not prolongated tracks (to be discussed) + + Int_t lab = TMath::Abs(tpcTrack->GetLabel()); + if (lab > nLab || lut[lab] < 0) continue; + + bckTree->GetEvent(lut[lab]); + bckTrack->Rotate(alpha-bckTrack->GetAlpha()); + + fSeeds->AddLast(new AliTPCseed(*bckTrack,bckTrack->GetAlpha())); + tracks.AddLast(new AliTPCtrack(*tpcTrack)); + } } + // end of matching + out->cd(); - TTree backTree("TreeT_TPCb_0","Tree with back propagated TPC tracks"); + + // tree name seedsTPCtoTRD as expected by TRD and as + // discussed and decided in Strasbourg (May 2002) + // [SR, GSI, 18.02.2003] + + sprintf(tName,"seedsTPCtoTRD_%d",GetEventNumber()); + TTree backTree(tName,"Tree with back propagated TPC tracks"); + AliTPCtrack *otrack=0; backTree.Branch("tracks","AliTPCtrack",&otrack,32000,0); -//*** Back propagation through inner sectors - LoadInnerSectors(); + // + Int_t nRefPlane; + Int_t found=0; Int_t nseed=fSeeds->GetEntriesFast(); + + // loop changed [SR, 01.04.2003] + for (i=0; iClear(); + AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps; const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt; + ps->ResetNWrong(); + ps->ResetNRotation(); + + if (IsStoringBarrel()) StoreBarrelTrack(ps, 1, 1); + + // Load outer sectors + fSectors=fInnerSec; + fN=fkNIS; + Int_t nc=t.GetNumberOfClusters(); - s.SetLabel(nc-1); //set number of the cluster to start with + //s.SetLabel(nc-1); //set number of the cluster to start with + s.SetNumber(nc); - if (FollowBackProlongation(s,t)) { - UseClusters(&s); + if (FollowBackProlongation(s,t)) UseClusters(&s); + else { + fSeeds->RemoveAt(i); continue; } - delete fSeeds->RemoveAt(i); - } - UnloadInnerSectors(); -//*** Back propagation through outer sectors - LoadOuterSectors(); - Int_t found=0; - for (i=0; iUncheckedAt(i), &s=*ps; - if (!ps) continue; - Int_t nc=s.GetNumberOfClusters(); - const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt; + if (IsStoringBarrel()) StoreBarrelTrack(ps, 2, 1); + if (IsStoringBarrel()) StoreBarrelTrack(ps, 3, 1); + + // Load outer sectors + fSectors=fOuterSec; + fN=fkNOS; + + nc=s.GetNumberOfClusters(); Double_t alpha=s.GetAlpha() - fSectors->GetAlphaShift(); if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi(); @@ -860,24 +1372,44 @@ Int_t AliTPCtracker::PropagateBack() alpha-=s.GetAlpha(); if (s.Rotate(alpha)) { - if (FollowBackProlongation(s,t)) { - if (s.GetNumberOfClusters() >= Int_t(0.4*nrows)) { - s.CookdEdx(); - s.SetLabel(t.GetLabel()); - UseClusters(&s,nc); - otrack=ps; - backTree.Fill(); - cerr<= Int_t(0.4*nrows)) { + s.CookdEdx(); + s.SetLabel(t.GetLabel()); + UseClusters(&s,nc); + + // Propagate to outer reference plane for comparison reasons + // reason for keeping fParam object [SR, GSI, 18.02.2003] + + ps->PropagateTo(fParam->GetOuterRadiusUp()); + otrack=ps; + backTree.Fill(); + + if (IsStoringBarrel()) StoreBarrelTrack(ps, 4, 1); + cerr<Fill(); delete fSeeds->RemoveAt(i); } - UnloadOuterSectors(); + + out->cd(); backTree.Write(); + // [SR, 01.04.2003] + if (IsStoringBarrel()) { + fBarrelFile->cd(); + fBarrelTree->Write(); + fBarrelFile->Flush(); + + delete fBarrelArray; + delete fBarrelTree; + } + + savedir->cd(); cerr<<"Number of seeds: "<>16; Int_t ncl=(index&0x0000ffff)>>00; - AliTPCClustersRow *clrow=((AliTPCtracker *) this)->fClustersArray.GetRow(sec,row); - return (AliCluster*)(*clrow)[ncl]; + const AliTPCcluster *cl=0; + + if (secGetNRowLow(); -// cout<<"par->GetNRowLow() = "<GetNRowLow()<<" fN = "<GetPadRowRadiiLow(i)); } else { @@ -993,19 +1535,42 @@ void AliTPCtracker::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) { } //_________________________________________________________________________ -void -AliTPCtracker::AliTPCRow::InsertCluster(const AliTPCcluster* c, UInt_t index) { +void AliTPCtracker:: +AliTPCRow::InsertCluster(const AliTPCcluster* c, Int_t sec, Int_t row) { //----------------------------------------------------------------------- // Insert a cluster into this pad row in accordence with its y-coordinate //----------------------------------------------------------------------- if (fN==kMaxClusterPerRow) { - cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return; + ::Error("InsertCluster","Too many clusters !"); + return; + } + + Int_t index=(((sec<<8)+row)<<16)+fN; + + if (fN==0) { + fSize=kMaxClusterPerRow/8; + fClusterArray=new AliTPCcluster[fSize]; + fIndex[0]=index; + fClusterArray[0]=*c; fClusters[fN++]=fClusterArray; + return; } - if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;} + + if (fN==fSize) { + Int_t size=fSize*2; + AliTPCcluster *buff=new AliTPCcluster[size]; + memcpy(buff,fClusterArray,fSize*sizeof(AliTPCcluster)); + for (Int_t i=0; iGetY()); memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCcluster*)); memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t)); - fIndex[i]=index; fClusters[i]=c; fN++; + fIndex[i]=index; + fClusters[i]=fClusterArray+fN; fClusterArray[fN++]=*c; } //___________________________________________________________________ @@ -1031,17 +1596,24 @@ void AliTPCtracker::AliTPCseed::CookdEdx(Double_t low, Double_t up) { //----------------------------------------------------------------- Int_t i; Int_t nc=GetNumberOfClusters(); - Int_t * index = new Int_t[nc]; - TMath::Sort(nc, fdEdxSample,index,kFALSE); + + Int_t swap;//stupid sorting + do { + swap=0; + for (i=0; i= fN)) return NULL; + return fClusterArray+i; + } UInt_t GetIndex(Int_t i) const {return fIndex[i];} Int_t Find(Double_t y) const; void SetX(Double_t x) {fX=x;} @@ -62,6 +74,8 @@ public: private: Int_t fN; //number of clusters const AliTPCcluster *fClusters[kMaxClusterPerRow]; //pointers to clusters + Int_t fSize; //size of array of clusters + AliTPCcluster *fClusterArray; //array of clusters UInt_t fIndex[kMaxClusterPerRow]; //indeces of clusters Double_t fX; //X-coordinate of this row @@ -73,7 +87,7 @@ public: //**************** Internal tracker class ********************** class AliTPCSector { public: - AliTPCSector() { fN=0; fRow = 0;} + AliTPCSector() { fN=0; fRow = 0; } ~AliTPCSector() { delete[] fRow; } AliTPCRow& operator[](Int_t i) const { return *(fRow+i); } Int_t GetNRows() const { return fN; } @@ -113,8 +127,8 @@ public: Double_t fAlpha; //opening angle Double_t fAlphaShift; //shift angle; Double_t fPadPitchWidth; //pad pitch width - Double_t f1PadPitchLength; //pad pitch length - Double_t f2PadPitchLength; //pad pitch length + Double_t f1PadPitchLength; //pad pitch length + Double_t f2PadPitchLength; //pad pitch length private: AliTPCSector(const AliTPCSector &s); //dummy copy contructor AliTPCSector& operator=(const AliTPCSector &s);//dummy assignment operator @@ -141,25 +155,38 @@ public: }; private: + void MakeSeeds(Int_t i1, Int_t i2); Int_t FollowProlongation(AliTPCseed& t, Int_t rf=0); Int_t FollowBackProlongation(AliTPCseed &s, const AliTPCtrack &t); + Int_t FollowRefitInward(AliTPCseed *seed, AliTPCtrack *track); AliTPCtracker(const AliTPCtracker& r); //dummy copy constructor AliTPCtracker &operator=(const AliTPCtracker& r);//dummy assignment operator const Int_t fkNIS; //number of inner sectors - AliTPCSector *fInnerSec; //array of inner sectors; + AliTPCSector *fInnerSec; //array of inner sectors const Int_t fkNOS; //number of outer sectors - AliTPCSector *fOuterSec; //array of outer sectors; - - Int_t fN; //number of loaded sectors - AliTPCSector *fSectors; //pointer to loaded sectors; - - Int_t fEventN; //event number - AliTPCClustersArray fClustersArray; //array of TPC clusters - TObjArray *fSeeds; //array of track seeds - TString fEvFolderName;//! name of data folder + AliTPCSector *fOuterSec; //array of outer sectors + + Int_t fN; //number of "active" sectors + AliTPCSector *fSectors; //pointer to "active" sectors; + + AliTPCParam *fParam; //! TPC parameters for outer reference plane [SR, GSI, 18.02.2003] + TObjArray *fSeeds; //array of track seeds + + // [SR, 01.04.2003] + void SetBarrelTree(const char *mode); + void StoreBarrelTrack(AliTPCtrack *ps, Int_t refPlane, Int_t isIn); + + // [SR, 01.04.2003] + TFile *fBarrelFile; // file with "barrel" tracks + TTree *fBarrelTree; // tree with "barrel" tracks + TBranch *fBarrelBranch; + TClonesArray *fBarrelArray; + AliBarrelTrack *fBarrelTrack; + + ClassDef(AliTPCtracker,1) // Time Projection Chamber tracker }; #endif diff --git a/TPC/TPCLinkDef.h b/TPC/TPCLinkDef.h index ea252a50b6d..7f5e277e4e8 100644 --- a/TPC/TPCLinkDef.h +++ b/TPC/TPCLinkDef.h @@ -17,6 +17,7 @@ #pragma link C++ class AliTPChit+; #pragma link C++ class AliTPCdigit+; #pragma link C++ class AliTPCcluster+; +#pragma link C++ class AliTPCclusterer+; #pragma link C++ class AliTPCtrack+; #pragma link C++ class AliTPCtracker+; @@ -55,7 +56,6 @@ #pragma link C++ class AliHitInfo+; -#pragma link C++ class AliTPCclusterer-; #pragma link C++ class AliTPCDigitizer; #pragma link C++ class AliTPCtrackerParam; diff --git a/TRD/AliTRDclusterizer.cxx b/TRD/AliTRDclusterizer.cxx index e3f018a2859..58fcab22928 100644 --- a/TRD/AliTRDclusterizer.cxx +++ b/TRD/AliTRDclusterizer.cxx @@ -239,21 +239,23 @@ Bool_t AliTRDclusterizer::WriteClusters(Int_t det) printf("AliTRDclusterizer::WriteClusters -- "); printf("Writing the cluster tree %-18s for event %d.\n" ,fClusterTree->GetName(),fEvent); - + /* fClusterTree->Write(); - AliTRDgeometry *geo = fTRD->GetGeometry(); geo->SetName("TRDgeometry"); geo->Write(); - fPar->Write(); - + fPar->Write(); + */ + AliLoader* loader = fRunLoader->GetLoader("TRDLoader"); + loader->WriteRecPoints("OVERWRITE"); + return kTRUE; } - + /* AliLoader* loader = fRunLoader->GetLoader("TRDLoader"); loader->WriteDigits("OVERWRITE"); - + */ printf("AliTRDclusterizer::WriteClusters -- "); printf("Unexpected detector index %d.\n",det); diff --git a/TRD/AliTRDdigits2cluster.C b/TRD/AliTRDdigits2cluster.C index 2c36e1b3e55..bba873248fa 100644 --- a/TRD/AliTRDdigits2cluster.C +++ b/TRD/AliTRDdigits2cluster.C @@ -32,8 +32,8 @@ void AliTRDdigits2cluster() clusterizer->SetVerbose(1); // Open the AliRoot file - // clusterizer->Open(infile,0); - clusterizer->Open(infile,outfile,0); + clusterizer->Open(infile,0); + //clusterizer->Open(infile,outfile,0); // Load the digits diff --git a/TRD/AliTRDtrack.cxx b/TRD/AliTRDtrack.cxx index 32d36d370cf..60b2e33bd32 100644 --- a/TRD/AliTRDtrack.cxx +++ b/TRD/AliTRDtrack.cxx @@ -13,7 +13,66 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* $Id$ */ +/* +$Log$ +Revision 1.20.2.1 2003/07/14 09:19:33 hristov +TOF included in the combined PID (Yu.Belikov) + +Revision 1.20 2003/05/27 17:46:13 hristov +TRD PID included in the ESD schema (T.Kuhr) + +Revision 1.19 2003/05/22 10:46:46 hristov +Using access methods instead of data members + +Revision 1.18 2003/04/10 10:36:54 hristov +Code for unified TPC/TRD tracking (S.Radomski) + +Revision 1.17 2003/02/19 09:02:28 hristov +Track time measurement (S.Radomski) + +Revision 1.16 2003/02/10 14:06:10 cblume +Add tracking without tilted pads as option + +Revision 1.15 2003/01/27 16:34:49 cblume +Update of tracking by Sergei and Chuncheng + +Revision 1.14 2002/11/07 15:52:09 cblume +Update of tracking code for tilted pads + +Revision 1.13 2002/10/22 15:53:08 alibrary +Introducing Riostream.h + +Revision 1.12 2002/10/14 14:57:44 hristov +Merging the VirtualMC branch to the main development branch (HEAD) + +Revision 1.8.10.2 2002/07/24 10:09:31 alibrary +Updating VirtualMC + +RRevision 1.11 2002/06/13 12:09:58 hristov +Minor corrections + +Revision 1.10 2002/06/12 09:54:35 cblume +Update of tracking code provided by Sergei + +Revision 1.8 2001/05/30 12:17:47 hristov +Loop variables declared once + +Revision 1.7 2001/05/28 17:07:58 hristov +Last minute changes; ExB correction in AliTRDclusterizerV1; taking into account of material in G10 TEC frames and material between TEC planes (C.Blume,S.Sedykh) + +Revision 1.4 2000/12/08 16:07:02 cblume +Update of the tracking by Sergei + +Revision 1.3 2000/10/15 23:40:01 cblume +Remove AliTRDconst + +Revision 1.2 2000/10/06 16:49:46 cblume +Made Getters const + +Revision 1.1.2.1 2000/09/22 14:47:52 cblume +Add the tracking code + +*/ #include #include @@ -226,6 +285,12 @@ AliTRDtrack::AliTRDtrack(const AliESDtrack& t) fdQdl[i] = 0; fIndex[i] = 0; } + + if ((t.GetStatus()&AliESDtrack::kTIME) == 0) return; + StartTimeIntegral(); + Double_t times[10]; t.GetIntegratedTimes(times); SetIntegratedTimes(times); + SetIntegratedLength(t.GetIntegratedLength()); + } //_____________________________________________________________________________ diff --git a/TRD/AliTRDtracker.cxx b/TRD/AliTRDtracker.cxx index 12b8d483211..bf239bb7e81 100644 --- a/TRD/AliTRDtracker.cxx +++ b/TRD/AliTRDtracker.cxx @@ -13,82 +13,147 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* $Id$ */ +/* +$Log$ +Revision 1.27.2.2 2003/07/14 09:19:33 hristov +TOF included in the combined PID (Yu.Belikov) + +Revision 1.27.2.1 2003/07/11 10:53:01 hristov +Inward refit for TPC and TRD in the ESD schema (T.Kuhr) + +Revision 1.27 2003/05/27 17:46:13 hristov +TRD PID included in the ESD schema (T.Kuhr) + +Revision 1.26 2003/04/10 10:36:54 hristov +Code for unified TPC/TRD tracking (S.Radomski) + +Revision 1.25 2003/03/19 17:14:11 hristov +Load/UnloadClusters added to the base class and the derived classes changed correspondingly. Possibility to give 2 input files for ITS and TPC tracks in PropagateBack. TRD tracker uses fEventN from the base class (T.Kuhr) + +Revision 1.24 2003/02/19 09:02:28 hristov +Track time measurement (S.Radomski) + +Revision 1.23 2003/02/10 14:06:10 cblume +Add tracking without tilted pads as option + +Revision 1.22 2003/01/30 15:19:58 cblume +New set of parameters + +Revision 1.21 2003/01/27 16:34:49 cblume +Update of tracking by Sergei and Chuncheng + +Revision 1.20 2002/11/07 15:52:09 cblume +Update of tracking code for tilted pads + +Revision 1.19 2002/10/22 15:53:08 alibrary +Introducing Riostream.h + +Revision 1.18 2002/10/14 14:57:44 hristov +Merging the VirtualMC branch to the main development branch (HEAD) + +Revision 1.14.6.2 2002/07/24 10:09:31 alibrary +Updating VirtualMC + +Revision 1.17 2002/06/13 12:09:58 hristov +Minor corrections + +Revision 1.16 2002/06/12 09:54:36 cblume +Update of tracking code provided by Sergei + +Revision 1.14 2001/11/14 10:50:46 cblume +Changes in digits IO. Add merging of summable digits + +Revision 1.13 2001/05/30 12:17:47 hristov +Loop variables declared once + +Revision 1.12 2001/05/28 17:07:58 hristov +Last minute changes; ExB correction in AliTRDclusterizerV1; taking into account of material in G10 TEC frames and material between TEC planes (C.Blume,S.Sedykh) + +Revision 1.8 2000/12/20 13:00:44 cblume +Modifications for the HP-compiler + +Revision 1.7 2000/12/08 16:07:02 cblume +Update of the tracking by Sergei + +Revision 1.6 2000/11/30 17:38:08 cblume +Changes to get in line with new STEER and EVGEN + +Revision 1.5 2000/11/14 14:40:27 cblume +Correction for the Sun compiler (kTRUE and kFALSE) + +Revision 1.4 2000/11/10 14:57:52 cblume +Changes in the geometry constants for the DEC compiler + +Revision 1.3 2000/10/15 23:40:01 cblume +Remove AliTRDconst + +Revision 1.2 2000/10/06 16:49:46 cblume +Made Getters const + +Revision 1.1.2.2 2000/10/04 16:34:58 cblume +Replace include files by forward declarations + +Revision 1.1.2.1 2000/09/22 14:47:52 cblume +Add the tracking code + +*/ #include #include #include #include -#include -#include +#include #include "AliTRDgeometry.h" #include "AliTRDparameter.h" #include "AliTRDgeometryDetail.h" #include "AliTRDcluster.h" #include "AliTRDtrack.h" +#include "AliTRDPartID.h" #include "../TPC/AliTPCtrack.h" #include "AliTRDtracker.h" ClassImp(AliTRDtracker) - const Float_t AliTRDtracker::fgkSeedDepth = 0.5; - const Float_t AliTRDtracker::fgkSeedStep = 0.10; - const Float_t AliTRDtracker::fgkSeedGap = 0.25; + const Float_t AliTRDtracker::fSeedDepth = 0.5; + const Float_t AliTRDtracker::fSeedStep = 0.10; + const Float_t AliTRDtracker::fSeedGap = 0.25; - const Float_t AliTRDtracker::fgkMaxSeedDeltaZ12 = 40.; - const Float_t AliTRDtracker::fgkMaxSeedDeltaZ = 25.; - const Float_t AliTRDtracker::fgkMaxSeedC = 0.0052; - const Float_t AliTRDtracker::fgkMaxSeedTan = 1.2; - const Float_t AliTRDtracker::fgkMaxSeedVertexZ = 150.; + const Float_t AliTRDtracker::fMaxSeedDeltaZ12 = 40.; + const Float_t AliTRDtracker::fMaxSeedDeltaZ = 25.; + const Float_t AliTRDtracker::fMaxSeedC = 0.0052; + const Float_t AliTRDtracker::fMaxSeedTan = 1.2; + const Float_t AliTRDtracker::fMaxSeedVertexZ = 150.; - const Double_t AliTRDtracker::fgkSeedErrorSY = 0.2; - const Double_t AliTRDtracker::fgkSeedErrorSY3 = 2.5; - const Double_t AliTRDtracker::fgkSeedErrorSZ = 0.1; + const Double_t AliTRDtracker::fSeedErrorSY = 0.2; + const Double_t AliTRDtracker::fSeedErrorSY3 = 2.5; + const Double_t AliTRDtracker::fSeedErrorSZ = 0.1; - const Float_t AliTRDtracker::fgkMinClustersInSeed = 0.7; + const Float_t AliTRDtracker::fMinClustersInSeed = 0.7; - const Float_t AliTRDtracker::fgkMinClustersInTrack = 0.5; - const Float_t AliTRDtracker::fgkMinFractionOfFoundClusters = 0.8; - - const Float_t AliTRDtracker::fgkSkipDepth = 0.3; - const Float_t AliTRDtracker::fgkLabelFraction = 0.8; - const Float_t AliTRDtracker::fgkWideRoad = 20.; + const Float_t AliTRDtracker::fMinClustersInTrack = 0.5; + const Float_t AliTRDtracker::fMinFractionOfFoundClusters = 0.8; - const Double_t AliTRDtracker::fgkMaxChi2 = 12.; + const Float_t AliTRDtracker::fSkipDepth = 0.3; + const Float_t AliTRDtracker::fLabelFraction = 0.8; + const Float_t AliTRDtracker::fWideRoad = 20.; - const Int_t AliTRDtracker::kFirstPlane = 5; - const Int_t AliTRDtracker::kLastPlane = 17; -//____________________________________________________________________ -Int_t AliTRDtracker::PropagateBack() -{ -//Overrides pure virtual methods in AliTracker -//Left for responsible to make it compatible with NewIO + const Double_t AliTRDtracker::fMaxChi2 = 12.; - Error("PropagateBack","Not yet NewIO-ed"); - return 0; -} -//____________________________________________________________________ +const Int_t AliTRDtracker::kFirstPlane = 5; +const Int_t AliTRDtracker::kLastPlane = 17; -Int_t AliTRDtracker::Clusters2Tracks() -{ -//Overrides pure virtual methods in AliTracker -//Left for responsible to make it compatible with NewIO - Error("PropagateBack","Not yet NewIO-ed"); - return 0; -} //____________________________________________________________________ - -AliTRDtracker::AliTRDtracker(const TFile *geomfile) +AliTRDtracker::AliTRDtracker(const TFile *geomfile):AliTracker() { // // Main constructor // - Float_t fTzero = 0; + //Float_t fTzero = 0; fAddTRDseeds = kFALSE; fGeom = NULL; @@ -102,10 +167,10 @@ AliTRDtracker::AliTRDtracker(const TFile *geomfile) } else { in->cd(); - in->ls(); +// in->ls(); fGeom = (AliTRDgeometry*) in->Get("TRDgeometry"); fPar = (AliTRDparameter*) in->Get("TRDparameter"); - fGeom->Dump(); +// fGeom->Dump(); } if(fGeom) { @@ -113,20 +178,23 @@ AliTRDtracker::AliTRDtracker(const TFile *geomfile) printf("Found geometry version %d on file \n", fGeom->IsVersion()); } else { - printf("AliTRDtracker::AliTRDtracker(): cann't find TRD geometry!\n"); - printf(" DETAIL TRD geometry and DEFAULT TRD parameter will be used\n"); - fGeom = new AliTRDgeometryDetail(); + printf("AliTRDtracker::AliTRDtracker(): can't find TRD geometry!\n"); + printf("The DETAIL TRD geometry will be used\n"); + fGeom = new AliTRDgeometryDetail(); + } + + if (!fPar) { + printf("AliTRDtracker::AliTRDtracker(): can't find TRD parameter!\n"); + printf("The DEFAULT TRD parameter will be used\n"); fPar = new AliTRDparameter(); } + fPar->ReInit(); savedir->cd(); // fGeom->SetT0(fTzero); - // fEvent = 0; - AliTracker::SetEventNumber(0); - fNclusters = 0; fClusters = new TObjArray(2000); fNseeds = 0; @@ -166,7 +234,7 @@ AliTRDtracker::AliTRDtracker(const TFile *geomfile) tbAmp = TMath::Min(tbAmp,maxAmp); fTimeBinsPerPlane = tbAmp + tbDrift; - fMaxGap = (Int_t) (fTimeBinsPerPlane * fGeom->Nplan() * fgkSkipDepth); + fMaxGap = (Int_t) (fTimeBinsPerPlane * fGeom->Nplan() * fSkipDepth); fVocal = kFALSE; @@ -281,7 +349,7 @@ Bool_t AliTRDtracker::AdjustSector(AliTRDtrack *track) { Double_t y = track->GetY(); Double_t ymax = track->GetX()*TMath::Tan(0.5*alpha); - Int_t ns = AliTRDgeometry::kNsect; + //Int_t ns = AliTRDgeometry::kNsect; //Int_t s=Int_t(track->GetAlpha()/alpha)%ns; if (y > ymax) { @@ -312,7 +380,6 @@ inline Double_t f1trd(Double_t x1,Double_t y1, Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b); return -xr*yr/sqrt(xr*xr+yr*yr); - } //_____________________________________________________________________ @@ -334,7 +401,6 @@ inline Double_t f2trd(Double_t x1,Double_t y1, Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b); return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr); - } //_____________________________________________________________________ @@ -347,7 +413,6 @@ inline Double_t f3trd(Double_t x1,Double_t y1, // return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)); - } //___________________________________________________________________ @@ -383,12 +448,13 @@ Int_t AliTRDtracker::Clusters2Tracks(const TFile *inp, TFile *out) trd_tree.Branch("tracks","AliTRDtrack",&iotrack_trd,32000,0); Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins(); - Float_t foundMin = fgkMinClustersInTrack * timeBins; + Float_t foundMin = fMinClustersInTrack * timeBins; if (inp) { TFile *in=(TFile*)inp; if (!in->IsOpen()) { - cerr<<"AliTRDtracker::Clusters2Tracks(): file with seeds is not open !\n"; + cerr<< + "AliTRDtracker::Clusters2Tracks(): file with seeds is not open !\n"; cerr<<" ... going for seeds finding inside the TRD\n"; } else { @@ -430,7 +496,7 @@ Int_t AliTRDtracker::Clusters2Tracks(const TFile *inp, TFile *out) FollowProlongation(t, innerTB); if (t.GetNumberOfClusters() >= foundMin) { UseClusters(&t); - CookLabel(pt, 1-fgkLabelFraction); + CookLabel(pt, 1-fLabelFraction); // t.CookdEdx(); } iotrack_trd = pt; @@ -458,15 +524,15 @@ Int_t AliTRDtracker::Clusters2Tracks(const TFile *inp, TFile *out) // Find tracks for the seeds in the TRD Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins(); - Int_t nSteps = (Int_t) (fgkSeedDepth / fgkSeedStep); - Int_t gap = (Int_t) (timeBins * fgkSeedGap); - Int_t step = (Int_t) (timeBins * fgkSeedStep); + Int_t nSteps = (Int_t) (fSeedDepth / fSeedStep); + Int_t gap = (Int_t) (timeBins * fSeedGap); + Int_t step = (Int_t) (timeBins * fSeedStep); // make a first turn with tight cut on initial curvature for(Int_t turn = 1; turn <= 2; turn++) { if(turn == 2) { - nSteps = (Int_t) (fgkSeedDepth / (3*fgkSeedStep)); - step = (Int_t) (timeBins * (3*fgkSeedStep)); + nSteps = (Int_t) (fSeedDepth / (3*fSeedStep)); + step = (Int_t) (timeBins * (3*fSeedStep)); } for(Int_t i=0; i= foundMin) { UseClusters(&t); - CookLabel(pt, 1-fgkLabelFraction); + CookLabel(pt, 1-fLabelFraction); t.CookdEdx(); - found++; + found++; // cout<GetNumberOfTimeBins(); + Float_t foundMin = fMinClustersInTrack * timeBins; + Int_t nseed = 0; + Int_t found = 0; + Int_t innerTB = fTrSec[0]->GetInnerTimeBin(); + + Int_t n = event->GetNumberOfTracks(); + for (Int_t i=0; iGetTrack(i); + ULong_t status=seed->GetStatus(); + if ( (status & AliESDtrack::kTRDout ) == 0 ) continue; + if ( (status & AliESDtrack::kTRDin) != 0 ) continue; + nseed++; + + AliTRDtrack* seed2 = new AliTRDtrack(*seed); + seed2->ResetCovariance(); + AliTRDtrack *pt = new AliTRDtrack(*seed2,seed2->GetAlpha()); + AliTRDtrack &t=*pt; + FollowProlongation(t, innerTB); + if (t.GetNumberOfClusters() >= foundMin) { + UseClusters(&t); + CookLabel(pt, 1-fLabelFraction); + // t.CookdEdx(); + } + found++; +// cout<UpdateTrackParams(pt, AliESDtrack::kTRDin); + } + delete seed2; + delete pt; + } + + cout<<"Number of loaded seeds: "<GetNumberOfTimeBins(); + + Int_t nSteps = (Int_t) (fSeedDepth / fSeedStep); + Int_t gap = (Int_t) (timeBins * fSeedGap); + Int_t step = (Int_t) (timeBins * fSeedStep); + + // make a first turn with tight cut on initial curvature + for(Int_t turn = 1; turn <= 2; turn++) { + if(turn == 2) { + nSteps = (Int_t) (fSeedDepth / (3*fSeedStep)); + step = (Int_t) (timeBins * (3*fSeedStep)); + } + for(Int_t i=0; iGetEntriesFast(); + + MakeSeeds(inner, outer, turn); + + nseed=fSeeds->GetEntriesFast(); + printf("\n turn %d, step %d: number of seeds for TRD inward %d\n", + turn, i, nseed); + + for (Int_t i=0; iUncheckedAt(i), &t=*pt; + FollowProlongation(t,innerTB); + if (t.GetNumberOfClusters() >= foundMin) { + UseClusters(&t); + CookLabel(pt, 1-fLabelFraction); + t.CookdEdx(); + found++; +// cout<AddTrack(&track); + } + } + delete fSeeds->RemoveAt(i); + fNseeds--; + } + } + } + } + + cout<<"Total number of found tracks: "<GetEntriesFast(); - // Float_t foundMin = fgkMinClustersInTrack * fTimeBinsPerPlane * fGeom->Nplan(); + // Float_t foundMin = fMinClustersInTrack * fTimeBinsPerPlane * fGeom->Nplan(); Float_t foundMin = 40; Int_t outermost_tb = fTrSec[0]->GetOuterTimeBin(); @@ -628,7 +799,7 @@ Int_t AliTRDtracker::PropagateBack(const TFile *inp, TFile *out) { if (foundClr >= foundMin) { if(foundClr >= 2) { s.CookdEdx(); - CookLabel(ps, 1-fgkLabelFraction); + CookLabel(ps, 1-fLabelFraction); UseClusters(ps); } @@ -643,7 +814,7 @@ Int_t AliTRDtracker::PropagateBack(const TFile *inp, TFile *out) { if(((expectedClr < 10) && (last_tb == outermost_tb)) || ((expectedClr >= 10) && (((Float_t) foundClr) / ((Float_t) expectedClr) >= - fgkMinFractionOfFoundClusters) && (last_tb == outermost_tb))) { + fMinFractionOfFoundClusters) && (last_tb == outermost_tb))) { Double_t x_tof = 375.5; @@ -707,6 +878,63 @@ Int_t AliTRDtracker::PropagateBack(const TFile *inp, TFile *out) { } +//_____________________________________________________________________________ +Int_t AliTRDtracker::PropagateBack(AliESD* event) { + // + // Gets seeds from ESD event. The seeds are AliTPCtrack's found and + // backpropagated by the TPC tracker. Each seed is first propagated + // to the TRD, and then its prolongation is searched in the TRD. + // If sufficiently long continuation of the track is found in the TRD + // the track is updated, otherwise it's stored as originaly defined + // by the TPC tracker. + // + + Int_t found=0; + Float_t foundMin = 40; + + Int_t n = event->GetNumberOfTracks(); + for (Int_t i=0; iGetTrack(i); + ULong_t status=seed->GetStatus(); + if ( (status & AliESDtrack::kTPCout ) == 0 ) continue; + if ( (status & AliESDtrack::kTRDout) != 0 ) continue; + + Int_t lbl = seed->GetLabel(); + AliTRDtrack *track = new AliTRDtrack(*seed); + track->SetSeedLabel(lbl); + fNseeds++; + + /*Int_t expectedClr = */FollowBackProlongation(*track); + + Int_t foundClr = track->GetNumberOfClusters(); + if (foundClr >= foundMin) { + if(foundClr >= 2) { + track->CookdEdx(); +// CookLabel(track, 1-fLabelFraction); + UseClusters(track); + } + + // Propagate to outer reference plane [SR, GSI, 18.02.2003] +// track->PropagateTo(364.8); why? + + //seed->UpdateTrackParams(track, AliESDtrack::kTRDout); + //found++; + } + + if (track->PropagateTo(376.)) { //Propagation to the TOF (I.Belikov) + seed->UpdateTrackParams(track, AliESDtrack::kTRDout); + found++; + } + + } + + cerr<<"Number of seeds: "<fgkWideRoad) { + if (road>fWideRoad) { if (t.GetNumberOfClusters()>4) cerr<GetLayer(nr)->IsSensitive() != fTrSec[s]->GetLayer(nr+1)->IsSensitive() ) { - if (IsStoringBarrel()) StoreBarrelTrack(&t, nRefPlane++, kTrackBack); +// if (IsStoringBarrel()) StoreBarrelTrack(&t, nRefPlane++, kTrackBack); } if (fTrSec[s]->GetLayer(nr-1)->IsSensitive() && @@ -1053,7 +1281,7 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t) wSigmaZ2 = (Float_t) t.GetSigmaZ2(); wChi2 = -1; - if (road>fgkWideRoad) { + if (road>fWideRoad) { if (t.GetNumberOfClusters()>4) cerr< 0) max_chi2 = 3 * max_chi2; @@ -1274,15 +1502,15 @@ Int_t AliTRDtracker::PropagateToTPC(AliTRDtrack& t) // padrow of the TPC. // Returns 1 if track reaches the TPC, and 0 otherwise - Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5); + //Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5); Double_t alpha=t.GetAlpha(); - TVector2::Phi_0_2pi(alpha); + alpha = TVector2::Phi_0_2pi(alpha); Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect; Bool_t lookForCluster; - Double_t rad_length, rho, x, dx, y, ymax, z; + Double_t rad_length, rho, x, dx, y, /*ymax,*/ z; x = t.GetX(); @@ -1316,7 +1544,6 @@ Int_t AliTRDtracker::PropagateToTPC(AliTRDtrack& t) return 1; } -//_____________________________________________________________________________ void AliTRDtracker::LoadEvent() { // Fills clusters into TRD tracking_sectors @@ -1348,6 +1575,42 @@ void AliTRDtracker::LoadEvent() } +//_____________________________________________________________________________ +Int_t AliTRDtracker::LoadClusters(TTree *cTree) +{ + // Fills clusters into TRD tracking_sectors + // Note that the numbering scheme for the TRD tracking_sectors + // differs from that of TRD sectors + + if (ReadClusters(fClusters,cTree)) { + Error("LoadClusters","Problem with reading the clusters !"); + return 1; + } + Int_t ncl=fClusters->GetEntriesFast(); + cout<<"\n LoadSectors: sorting "<UncheckedAt(ncl); + Int_t detector=c->GetDetector(), local_time_bin=c->GetLocalTimeBin(); + Int_t sector=fGeom->GetSector(detector); + Int_t plane=fGeom->GetPlane(detector); + + Int_t tracking_sector = CookSectorIndex(sector); + + Int_t gtb = fTrSec[tracking_sector]->CookTimeBinIndex(plane,local_time_bin); + if(gtb < 0) continue; + Int_t layer = fTrSec[tracking_sector]->GetLayerNumber(gtb); + + index=ncl; + fTrSec[tracking_sector]->GetLayer(layer)->InsertCluster(c,index); + } + printf("\r\n"); + + return 0; +} + //_____________________________________________________________________________ void AliTRDtracker::UnloadEvent() { @@ -1463,11 +1726,11 @@ void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer, Int_t turn) y2=xx2*sn2+y2*cs2; } - if(TMath::Abs(z1-z2) > fgkMaxSeedDeltaZ12) continue; + if(TMath::Abs(z1-z2) > fMaxSeedDeltaZ12) continue; Double_t zz=z1 - z1/x1*(x1-x2); - if (TMath::Abs(zz-z2)>fgkMaxSeedDeltaZ) continue; + if (TMath::Abs(zz-z2)>fMaxSeedDeltaZ) continue; Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1); if (d==0.) {cerr<<"TRD MakeSeeds: Straight seed !\n"; continue;} @@ -1476,7 +1739,7 @@ void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer, Int_t turn) x[1]=z1; x[4]=f1trd(x1,y1,x2,y2,x3,y3); - if (TMath::Abs(x[4]) > fgkMaxSeedC) continue; + if (TMath::Abs(x[4]) > fMaxSeedC) continue; x[2]=f2trd(x1,y1,x2,y2,x3,y3); @@ -1484,16 +1747,16 @@ void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer, Int_t turn) x[3]=f3trd(x1,y1,x2,y2,z1,z2); - if (TMath::Abs(x[3]) > fgkMaxSeedTan) continue; + if (TMath::Abs(x[3]) > fMaxSeedTan) continue; Double_t a=asin(x[2]); Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2])); - if (TMath::Abs(zv)>fgkMaxSeedVertexZ) continue; + if (TMath::Abs(zv)>fMaxSeedVertexZ) continue; Double_t sy1=r1[is]->GetSigmaY2(), sz1=r1[is]->GetSigmaZ2(); Double_t sy2=cl->GetSigmaY2(), sz2=cl->GetSigmaZ2(); - Double_t sy3=fgkSeedErrorSY3, sy=fgkSeedErrorSY, sz=fgkSeedErrorSZ; + Double_t sy3=fSeedErrorSY3, sy=fSeedErrorSY, sz=fSeedErrorSZ; // Tilt changes Double_t h01 = GetTiltFactor(r1[is]); @@ -1537,7 +1800,7 @@ void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer, Int_t turn) if ((rc < 1) || (track->GetNumberOfClusters() < - (outer-inner)*fgkMinClustersInSeed)) delete track; + (outer-inner)*fMinClustersInSeed)) delete track; else { fSeeds->AddLast(track); fNseeds++; // cerr<<"\r found seed "<= 0) or AliTRDrecPoints (option < 0) // from the file. The names of the cluster tree and branches // should match the ones used in AliTRDclusterizer::WriteClusters() // - - TDirectory *savedir=gDirectory; - - if (inp) { - TFile *in=(TFile*)inp; - if (!in->IsOpen()) { - cerr<<"AliTRDtracker::ReadClusters(): input file is not open !\n"; - return; - } - else{ - in->cd(); - } - } - - Char_t treeName[12]; - sprintf(treeName,"TreeR%d_TRD",GetEventNumber()); - TTree *ClusterTree = (TTree*) gDirectory->Get(treeName); - TObjArray *ClusterArray = new TObjArray(400); - ClusterTree->GetBranch("TRDcluster")->SetAddress(&ClusterArray); + TBranch *branch=ClusterTree->GetBranch("TRDcluster"); + if (!branch) { + Error("ReadClusters","Can't get the branch !"); + return 1; + } + branch->SetAddress(&ClusterArray); Int_t nEntries = (Int_t) ClusterTree->GetEntries(); printf("found %d entries in %s.\n",nEntries,ClusterTree->GetName()); @@ -1608,8 +1858,8 @@ void AliTRDtracker::ReadClusters(TObjArray *array, const TFile *inp) } delete ClusterArray; - savedir->cd(); + return 0; } //______________________________________________________________________ @@ -1681,6 +1931,69 @@ void AliTRDtracker::ReadClusters(TObjArray *array, const Char_t *filename) } +void AliTRDtracker::ReadClusters(TObjArray *array, const TFile *inp) +{ + // + // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0) + // from the file. The names of the cluster tree and branches + // should match the ones used in AliTRDclusterizer::WriteClusters() + // + + TDirectory *savedir=gDirectory; + + if (inp) { + TFile *in=(TFile*)inp; + if (!in->IsOpen()) { + cerr<<"AliTRDtracker::ReadClusters(): input file is not open !\n"; + return; + } + else{ + in->cd(); + } + } + + Char_t treeName[12]; + sprintf(treeName,"TreeR%d_TRD",GetEventNumber()); + TTree *ClusterTree = (TTree*) gDirectory->Get(treeName); + + TObjArray *ClusterArray = new TObjArray(400); + + ClusterTree->GetBranch("TRDcluster")->SetAddress(&ClusterArray); + + Int_t nEntries = (Int_t) ClusterTree->GetEntries(); + printf("found %d entries in %s.\n",nEntries,ClusterTree->GetName()); + + // Loop through all entries in the tree + Int_t nbytes; + AliTRDcluster *c = 0; + printf("\n"); + + for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) { + + // Import the tree + nbytes += ClusterTree->GetEvent(iEntry); + + // Get the number of points in the detector + Int_t nCluster = ClusterArray->GetEntriesFast(); +// printf("\r Read %d clusters from entry %d", nCluster, iEntry); + + // Loop through all TRD digits + for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { + c = (AliTRDcluster*)ClusterArray->UncheckedAt(iCluster); + AliTRDcluster *co = new AliTRDcluster(*c); + co->SetSigmaY2(c->GetSigmaY2() * fSY2corr); + Int_t ltb = co->GetLocalTimeBin(); + if(ltb == 19) co->SetSigmaZ2(c->GetSigmaZ2()); + else if(fNoTilt) co->SetSigmaZ2(c->GetSigmaZ2() * fSZ2corr); + array->AddLast(co); + delete ClusterArray->RemoveAt(iCluster); + } + } + + delete ClusterArray; + savedir->cd(); + +} //__________________________________________________________________ void AliTRDtracker::CookLabel(AliKalmanTrack* pt, Float_t wrong) const { @@ -1750,6 +2063,7 @@ void AliTRDtracker::CookLabel(AliKalmanTrack* pt, Float_t wrong) const { } + //__________________________________________________________________ void AliTRDtracker::UseClusters(const AliKalmanTrack* t, Int_t from) const { Int_t ncl=t->GetNumberOfClusters(); @@ -1794,6 +2108,7 @@ Double_t AliTRDtracker::GetX(Int_t sector, Int_t plane, Int_t local_tb) const } + //_______________________________________________________ AliTRDtracker::AliTRDpropagationLayer::AliTRDpropagationLayer(Double_t x, Double_t dx, Double_t rho, Double_t rad_length, Int_t tb_index) @@ -2440,3 +2755,8 @@ Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster* c) { return h01; } + + + + + diff --git a/TRD/AliTRDtracker.h b/TRD/AliTRDtracker.h index b9fecd1ca07..7a8e18748d8 100644 --- a/TRD/AliTRDtracker.h +++ b/TRD/AliTRDtracker.h @@ -4,16 +4,13 @@ /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * See cxx source for full Copyright notice */ -/* $Id$ */ - #include "AliTracker.h" #include "TObjArray.h" #include "AliBarrelTrack.h" - -#include -#include +#include "AliESD.h" class TFile; +class TTree; class TParticle; class TParticlePDG; @@ -29,7 +26,7 @@ const unsigned kMAX_CLUSTER_PER_TIME_BIN = 7000; const unsigned kZONES = 5; const Int_t kTRACKING_SECTORS = 18; -class AliTRDtracker : public AliTracker { +class AliTRDtracker : public AliTracker { public: @@ -38,13 +35,14 @@ class AliTRDtracker : public AliTracker { ~AliTRDtracker(); Int_t Clusters2Tracks(const TFile *in, TFile *out); + Int_t Clusters2Tracks(AliESD* event); Int_t PropagateBack(const TFile *in, TFile *out); - Int_t PropagateBack();//Almost empty see .cxx for more comments - Int_t Clusters2Tracks();//Almost empty see .cxx for more comments + Int_t PropagateBack(AliESD* event); + Int_t RefitInward(AliESD* event) {return 0;} - Int_t LoadClusters() {LoadEvent(); return 0;}; - void UnloadClusters() {UnloadEvent();}; - AliCluster *GetCluster(Int_t index) const { return (AliCluster*) fClusters->UncheckedAt(index); }; + Int_t LoadClusters(TTree *cTree); + void UnloadClusters(){UnloadEvent();} + AliCluster *GetCluster(Int_t index) const { if (index >= fNclusters) return NULL; return (AliCluster*) fClusters->UncheckedAt(index); }; virtual void CookLabel(AliKalmanTrack *t,Float_t wrong) const; virtual void UseClusters(const AliKalmanTrack *t, Int_t from=0) const; @@ -54,28 +52,29 @@ class AliTRDtracker : public AliTracker { Double_t GetTiltFactor(const AliTRDcluster* c); void ReadClusters(TObjArray *array, const Char_t *filename); + Int_t ReadClusters(TObjArray *array, TTree *in); Int_t CookSectorIndex(Int_t gs) { return kTRACKING_SECTORS - 1 - gs; } - Float_t GetSeedGap() const {return fgkSeedGap;} + Float_t GetSeedGap() const {return fSeedGap;} Int_t GetMaxGap() const {return fMaxGap;} Int_t GetTimeBinsPerPlane() const {return fTimeBinsPerPlane;} - Float_t GetSeedStep() const {return fgkSeedStep;} - Float_t GetSeedDepth() const {return fgkSeedDepth;} - Float_t GetSkipDepth() const {return fgkSkipDepth;} - Double_t GetMaxChi2() const {return fgkMaxChi2;} - Float_t GetMaxSeedC() const {return fgkMaxSeedC;} - Float_t GetMaxSeedTan() const {return fgkMaxSeedTan;} - Double_t GetSeedErrorSY() const {return fgkSeedErrorSY;} - Double_t GetSeedErrorSY3() const {return fgkSeedErrorSY3;} - Double_t GetSeedErrorSZ() const {return fgkSeedErrorSZ;} - Float_t GetLabelFraction() const {return fgkLabelFraction;} - Float_t GetWideRoad() const {return fgkWideRoad;} - - Float_t GetMinClustersInTrack() const {return fgkMinClustersInTrack;} - Float_t GetMinClustersInSeed() const {return fgkMinClustersInSeed;} - Float_t GetMaxSeedDeltaZ() const {return fgkMaxSeedDeltaZ;} - Float_t GetMaxSeedVertexZ() const {return fgkMaxSeedVertexZ;} + Float_t GetSeedStep() const {return fSeedStep;} + Float_t GetSeedDepth() const {return fSeedDepth;} + Float_t GetSkipDepth() const {return fSkipDepth;} + Double_t GetMaxChi2() const {return fMaxChi2;} + Float_t GetMaxSeedC() const {return fMaxSeedC;} + Float_t GetMaxSeedTan() const {return fMaxSeedTan;} + Double_t GetSeedErrorSY() const {return fSeedErrorSY;} + Double_t GetSeedErrorSY3() const {return fSeedErrorSY3;} + Double_t GetSeedErrorSZ() const {return fSeedErrorSZ;} + Float_t GetLabelFraction() const {return fLabelFraction;} + Float_t GetWideRoad() const {return fWideRoad;} + + Float_t GetMinClustersInTrack() const {return fMinClustersInTrack;} + Float_t GetMinClustersInSeed() const {return fMinClustersInSeed;} + Float_t GetMaxSeedDeltaZ() const {return fMaxSeedDeltaZ;} + Float_t GetMaxSeedVertexZ() const {return fMaxSeedVertexZ;} // x <-> timebin conversions useful in analysis macros Double_t GetX(Int_t sec, Int_t plane, Int_t local_tb) const; @@ -191,7 +190,7 @@ class AliTRDtracker : public AliTracker { Int_t PropagateToOuterPlane(AliTRDtrack& t, Double_t x); Int_t WriteTracks(); - void ReadClusters(TObjArray *array, const TFile *in = 0); + void ReadClusters(TObjArray *array, const TFile *in=0); void SetSY2corr(Float_t w) {fSY2corr = w;} void SetSZ2corr(Float_t w) {fSZ2corr = w;} @@ -222,41 +221,41 @@ class AliTRDtracker : public AliTracker { Float_t fSZ2corr; // Correction coefficient for // cluster SigmaZ2 - static const Float_t fgkSeedGap; // Distance between inner and outer + static const Float_t fSeedGap; // Distance between inner and outer // time bin in seeding // (fraction of all time bins) - static const Float_t fgkSeedStep; // Step in iterations - static const Float_t fgkSeedDepth; // Fraction of TRD allocated fgkor seeding - static const Float_t fgkSkipDepth; // Fraction of TRD which can be skipped - // in track prolongation + static const Float_t fSeedStep; // Step in iterations + static const Float_t fSeedDepth; // Fraction of TRD allocated for seeding + static const Float_t fSkipDepth; // Fraction of TRD which can be skipped + // in track prolongation Int_t fTimeBinsPerPlane; // number of sensitive timebins per plane Int_t fMaxGap; // max gap (in time bins) in the track - // in track prolongation + // in track prolongation - static const Double_t fgkMaxChi2; // max increment in track chi2 - - static const Float_t fgkMinClustersInTrack; // min number of clusters in track + static const Double_t fMaxChi2; // max increment in track chi2 + + static const Float_t fMinClustersInTrack; // min number of clusters in track // out of total timebins - static const Float_t fgkMinFractionOfFoundClusters; // min fgkound clusters + static const Float_t fMinFractionOfFoundClusters; // min found clusters // out of expected - static const Float_t fgkMinClustersInSeed; // min fgkraction of clusters in seed - static const Float_t fgkMaxSeedDeltaZ; // max dZ in MakeSeeds - static const Float_t fgkMaxSeedDeltaZ12; // max abs(z1-z2) in MakeSeeds - static const Float_t fgkMaxSeedC; // max initial curvature in MakeSeeds - static const Float_t fgkMaxSeedTan; // max initial Tangens(lambda) in MakeSeeds - static const Float_t fgkMaxSeedVertexZ; // max vertex Z in MakeSeeds - static const Double_t fgkSeedErrorSY; // sy parameter in MakeSeeds - static const Double_t fgkSeedErrorSY3; // sy3 parameter in MakeSeeds - static const Double_t fgkSeedErrorSZ; // sz parameter in MakeSeeds - static const Float_t fgkLabelFraction; // min fgkraction of same label - static const Float_t fgkWideRoad; // max road width in FindProlongation + static const Float_t fMinClustersInSeed; // min fraction of clusters in seed + static const Float_t fMaxSeedDeltaZ; // max dZ in MakeSeeds + static const Float_t fMaxSeedDeltaZ12; // max abs(z1-z2) in MakeSeeds + static const Float_t fMaxSeedC; // max initial curvature in MakeSeeds + static const Float_t fMaxSeedTan; // max initial Tangens(lambda) in MakeSeeds + static const Float_t fMaxSeedVertexZ; // max vertex Z in MakeSeeds + static const Double_t fSeedErrorSY; // sy parameter in MakeSeeds + static const Double_t fSeedErrorSY3; // sy3 parameter in MakeSeeds + static const Double_t fSeedErrorSZ; // sz parameter in MakeSeeds + static const Float_t fLabelFraction; // min fraction of same label + static const Float_t fWideRoad; // max road width in FindProlongation Bool_t fVocal; Bool_t fAddTRDseeds; - + Bool_t fNoTilt; -- 2.43.0