X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=TPC%2FAliTPCComparison.C;h=e14f7cadace427589c90b6c0a94973e152bfb980;hp=5bef125d55c1dacfa3bdbc4dd9671a75db615408;hb=162637e45086576ffcf0568d28fa43359819765a;hpb=f38c8ae57f453128ee780dc2a89a2c66a73e9985 diff --git a/TPC/AliTPCComparison.C b/TPC/AliTPCComparison.C index 5bef125d55c..e14f7cadace 100644 --- a/TPC/AliTPCComparison.C +++ b/TPC/AliTPCComparison.C @@ -1,289 +1,324 @@ -#ifndef __CINT__ - #include "alles.h" - #include "AliTPCtracker.h" +/**************************************************************************** + * Very important, delicate and rather obscure macro. * + * * + * Creates list of "trackable" tracks, * + * calculates efficiency, resolutions etc. * + * There is a possibility to run this macro over several events. * + * * + * Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch * + * with several nice improvements by: M.Ivanov, GSI, m.ivanov@gsi.de * + ****************************************************************************/ + +#if !defined(__CINT__) || defined(__MAKECINT__) + #include + #include + #include + #include + #include + #include + #include + #include + #include + #include + #include + #include + #include + #include + + #include "AliStack.h" + #include "AliHeader.h" + #include "AliTrackReference.h" + #include "AliRunLoader.h" + #include "AliRun.h" + #include "AliESDEvent.h" + #include "AliESDtrack.h" + + #include "AliSimDigits.h" + #include "AliTPC.h" + #include "AliTPCParamSR.h" + #include "AliTPCClustersArray.h" + #include "AliTPCClustersRow.h" + #include "AliTPCcluster.h" + #include "AliTPCLoader.h" + + #include "AliCDBManager.h" + #include "AliTPCcalibDB.h" #endif -struct GoodTrackTPC { - Int_t fEventN; //event number - Int_t lab; - Int_t code; - Float_t px,py,pz; - Float_t x,y,z; -}; +Int_t GoodTracksTPC(const Char_t *dir="."); -Int_t AliTPCComparison(Int_t eventn=1) { - Int_t good_tracks_tpc(GoodTrackTPC *gt, Int_t max, Int_t eventn=1); +extern TBenchmark *gBenchmark; +extern TROOT *gROOT; - cerr<<"Doing comparison...\n"; - Int_t i; +static Int_t allgood=0; +static Int_t allselected=0; +static Int_t allfound=0; + +Int_t AliTPCComparison +(Float_t ptcutl=0.2, Float_t ptcuth=10., const Char_t *dir=".") { gBenchmark->Start("AliTPCComparison"); - TFile *cf=TFile::Open("AliTPCclusters.root"); - if (!cf->IsOpen()) {cerr<<"Can't open AliTPCclusters.root !\n"; return 1;} - AliTPCParam *digp= (AliTPCParam*)cf->Get("75x40_100x60"); - if (!digp) { cerr<<"TPC parameters have not been found !\n"; return 2; } - - - - /////////// - AliTPCtracker *tracker =0; - TObjArray tarray(2000); - AliTPCtrack *iotrack=0; - Int_t nentr= 0; - Int_t eventptr[1000]; - TFile *tf=TFile::Open("AliTPCtracks.root"); - TTree *tracktree=0; - // Load clusters - eventptr[0]=0; - eventptr[1]=0; - char tname[100]; //Y.B. - for (Int_t event=0;eventcd(); - tracker = new AliTPCtracker(digp,event); - tracker->LoadInnerSectors(); - tracker->LoadOuterSectors(); - //Y.B. char tname[100]; - if (eventn==-1) { - sprintf(tname,"TreeT_TPC"); - } - else { - sprintf(tname,"TreeT_TPC_%d",event); - } + ::Info("AliTPCComparison.C","Doing comparison..."); - // Load tracks - if (!tf->IsOpen()) {cerr<<"Can't open AliTPCtracks.root !\n"; return 3;} - tracktree=(TTree*)tf->Get(tname); - if (!tracktree) {cerr<<"Can't get a tree with TPC tracks !\n"; return 4;} - TBranch *tbranch=tracktree->GetBranch("tracks"); - Int_t nentr0=(Int_t)tracktree->GetEntries(); - nentr+=nentr0; - for (i=0; iSetAddress(&iotrack); - tracktree->GetEvent(i); - tracker->CookLabel(iotrack,0.1); - tarray.AddLast(iotrack); - } - eventptr[event+1] = nentr; //store end of the event - delete tracker; - delete tracktree; //Thanks to Mariana Bondila - } - tf->Close(); - delete digp; //Thanks to Mariana Bondila - cf->Close(); - -///////////////////////////////////////////////////////////////////////// - const Int_t MAX=45000; - GoodTrackTPC gt[MAX]; - - Int_t ngood=0; - ifstream in("good_tracks_tpc"); - if (in) { - cerr<<"Reading good tracks...\n"; - while (in>>gt[ngood].fEventN>>gt[ngood].lab>>gt[ngood].code>> - gt[ngood].px>>gt[ngood].py>>gt[ngood].pz>> - gt[ngood].x >>gt[ngood].y >>gt[ngood].z) { - ngood++; - cerr<Branch("tracks","AliTPCtrack",&iotrack,32000,0); - for (i=0; iFill(); - } - tracktree->Write(); - tf->Close(); - } - cerr<<"Number of good tracks : "<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.); - hpt->SetFillColor(2); - TH1F *hmpt=new TH1F("hmpt","Relative Pt resolution (pt>4GeV/c)",30,-60,60); + TH1F *hp=(TH1F*)gROOT->FindObject("hp"); + if (!hp) hp=new TH1F("hp","PHI resolution",50,-20.,20.); + hp->SetFillColor(4); + + TH1F *hl=(TH1F*)gROOT->FindObject("hl"); + if (!hl) hl=new TH1F("hl","LAMBDA resolution",50,-20,20); + hl->SetFillColor(4); + + TH1F *hpt=(TH1F*)gROOT->FindObject("hpt"); + if (!hpt) hpt=new TH1F("hpt","Relative Pt resolution",30,-10.,10.); + hpt->SetFillColor(2); + + TH1F *hmpt=(TH1F*)gROOT->FindObject("hmpt"); + if (!hmpt) + hmpt=new TH1F("hmpt","Relative Pt resolution (pt>4GeV/c)",30,-60,60); hmpt->SetFillColor(6); - AliTPCtrack *trk=(AliTPCtrack*)tarray.UncheckedAt(0); - Double_t pmin=0.1*(100/0.299792458/0.2/trk->GetConvConst()); - Double_t pmax=6.0+pmin; - TH1F *hgood=new TH1F("hgood","Good tracks",30,pmin,pmax); - TH1F *hfound=new TH1F("hfound","Found tracks",30,pmin,pmax); - TH1F *hfake=new TH1F("hfake","Fake tracks",30,pmin,pmax); - TH1F *hg=new TH1F("hg","",30,pmin,pmax); //efficiency for good tracks + + TH1F *hgood=(TH1F*)gROOT->FindObject("hgood"); + if (!hgood) hgood=new TH1F("hgood","Good tracks",30,0.2,6.1); + + TH1F *hfound=(TH1F*)gROOT->FindObject("hfound"); + if (!hfound) hfound=new TH1F("hfound","Found tracks",30,0.2,6.1); + + TH1F *hfake=(TH1F*)gROOT->FindObject("hfake"); + if (!hfake) hfake=new TH1F("hfake","Fake tracks",30,0.2,6.1); + + TH1F *hg=(TH1F*)gROOT->FindObject("hg"); + if (!hg) hg=new TH1F("hg","Efficiency for good tracks",30,0.2,6.1); hg->SetLineColor(4); hg->SetLineWidth(2); - TH1F *hf=new TH1F("hf","Efficiency for fake tracks",30,pmin,pmax); + + TH1F *hf=(TH1F*)gROOT->FindObject("hf"); + if (!hf) hf=new TH1F("hf","Efficiency for fake tracks",30,0.2,6.1); hf->SetFillColor(1); hf->SetFillStyle(3013); hf->SetLineWidth(2); - TH1F *he =new TH1F("he","dE/dX for pions with 0.4FindObject("he"); + if (!he) + he =new TH1F("he","dE/dX for pions with 0.4FindObject("hep"); + if (!hep) hep=new TH2F("hep","dE/dX vs momentum",50,0.,2.,50,0.,400.); hep->SetMarkerStyle(8); hep->SetMarkerSize(0.4); - Int_t mingood = ngood; //MI change - Int_t * track_notfound = new Int_t[ngood]; - Int_t itrack_notfound =0; - Int_t * track_fake = new Int_t[ngood]; - Int_t itrack_fake = 0; - Int_t * track_multifound = new Int_t[ngood]; - Int_t * track_multifound_n = new Int_t[ngood]; - Int_t itrack_multifound =0; - while (ngood--) { - Int_t lab=gt[ngood].lab,tlab=-1; - Float_t ptg= - TMath::Sqrt(gt[ngood].px*gt[ngood].px + gt[ngood].py*gt[ngood].py); + Char_t fname[100]; + sprintf(fname,"%s/GoodTracksTPC.root",dir); + + TFile *refFile=TFile::Open(fname,"old"); + if (!refFile || !refFile->IsOpen()) { + ::Info("AliTPCComparison.C","Marking good tracks (will take a while)..."); + if (GoodTracksTPC(dir)) { + ::Error("AliTPCComparison.C","Can't generate the reference file !"); + return 1; + } + } + refFile=TFile::Open(fname,"old"); + if (!refFile || !refFile->IsOpen()) { + ::Error("AliTPCComparison.C","Can't open the reference file !"); + return 2; + } + + TTree *tpcTree=(TTree*)refFile->Get("tpcTree"); + if (!tpcTree) { + ::Error("AliTPCComparison.C","Can't get the reference tree !"); + return 3; + } + TBranch *branch=tpcTree->GetBranch("TPC"); + if (!branch) { + ::Error("AliTPCComparison.C","Can't get the TPC branch !"); + return 4; + } + TClonesArray dummy("AliTrackReference",1000), *refs=&dummy; + branch->SetAddress(&refs); + + + sprintf(fname,"%s/AliESDs.root",dir); + TFile *ef=TFile::Open(fname); + if ((!ef)||(!ef->IsOpen())) { + sprintf(fname,"%s/AliESDtpc.root",dir); + ef=TFile::Open(fname); + if ((!ef)||(!ef->IsOpen())) { + ::Error("AliTPCComparison.C","Can't open AliESDtpc.root !"); + return 5; + } + } + AliESDEvent* event = new AliESDEvent(); + TTree* esdTree = (TTree*) ef->Get("esdTree"); + if (!esdTree) { + ::Error("AliTPCComparison.C", "no ESD tree found"); + return 6; + } + event->ReadFromTree(esdTree); - if (ptgFill(ptg); - Int_t ievent = gt[ngood].fEventN; + //******* Loop over events ********* + Int_t e=0; + while (esdTree->GetEvent(e)) { + cout<GetNumberOfTracks(); + allfound+=nentr; - track=(AliTPCtrack*)tarray.UncheckedAt(i); - tlab=track->GetLabel(); - if (lab==TMath::Abs(tlab)) break; + if (tpcTree->GetEvent(e++)==0) { + cerr<<"No reconstructable tracks !\n"; + continue; } - // if (i==nentr) { - if (i==eventptr[ievent+1]) { - //cerr<<"Track "<GetLabel())) micount++; - } - if (micount>1) { - //cout<<"Track no. "<GetEntriesFast(); + allgood+=ngood; - // - Double_t xk=gt[ngood].x;//digp->GetPadRowRadii(0,0); - track->PropagateTo(xk); + const Int_t MAX=15000; + //MI change + Int_t track_notfound[MAX], itrack_notfound=0; + Int_t track_fake[MAX], itrack_fake=0; + Int_t track_multifound[MAX],track_multifound_n[MAX],itrack_multifound=0; - if (lab==tlab) hfound->Fill(ptg); - else { - track_fake[itrack_fake++]=lab; - hfake->Fill(ptg); - //cerr<UncheckedAt(k); + Int_t lab=ref->Label(), tlab=-1; + Float_t ptg=TMath::Sqrt(ref->Px()*ref->Px() + ref->Py()*ref->Py()); - Double_t par[5]; track->GetExternalParameters(xk,par); - Float_t phi=TMath::ASin(par[2]) + track->GetAlpha(); - if (phi<-TMath::Pi()) phi+=2*TMath::Pi(); - if (phi>=TMath::Pi()) phi-=2*TMath::Pi(); - Float_t lam=TMath::ATan(par[3]); - Float_t pt_1=TMath::Abs(par[4]); + if (ptg<1e-33) continue; // for those not crossing 0 pad row + + if (ptgptcuth) continue; - if (TMath::Abs(gt[ngood].code)==11 && ptg>4.) { - hmpt->Fill((pt_1 - 1/ptg)/(1/ptg)*100.); - } else { - Float_t phig=TMath::ATan2(gt[ngood].py,gt[ngood].px); - hp->Fill((phi - phig)*1000.); + allselected++; - Float_t lamg=TMath::ATan2(gt[ngood].pz,ptg); - hl->Fill((lam - lamg)*1000.); + hgood->Fill(ptg); - hpt->Fill((pt_1 - 1/ptg)/(1/ptg)*100.); + AliESDtrack *track=0; + for (i=0; iGetTrack(i); + tlab=track->GetTPCLabel(); + if (lab==TMath::Abs(tlab)) break; + } + if (i==nentr) { + track_notfound[itrack_notfound++]=lab; + continue; + } + + //MI change - addition + Int_t micount=0; + Int_t mi; + AliESDtrack * mitrack; + for (mi=0; miGetTrack(mi); + if (lab==TMath::Abs(mitrack->GetTPCLabel())) micount++; + } + if (micount>1) { + track_multifound[itrack_multifound]=lab; + track_multifound_n[itrack_multifound]=micount; + itrack_multifound++; + } + if ((track->GetStatus()&AliESDtrack::kTPCrefit)==0) continue; + if (lab==tlab) hfound->Fill(ptg); + else { + track_fake[itrack_fake++]=lab; + hfake->Fill(ptg); + } + + Double_t pxpypz[3]; track->GetInnerPxPyPz(pxpypz); + Float_t phi=TMath::ATan2(pxpypz[1],pxpypz[0]); + if (phi<-TMath::Pi()) phi+=2*TMath::Pi(); + if (phi>=TMath::Pi()) phi-=2*TMath::Pi(); + Double_t pt=TMath::Sqrt(pxpypz[0]*pxpypz[0]+pxpypz[1]*pxpypz[1]); + Float_t lam=TMath::ATan2(pxpypz[2],pt); + Float_t pt_1=1/pt; + + Int_t pdg=(Int_t)ref->GetLength(); //this is particle's PDG ! + + if (TMath::Abs(pdg)==11 && ptg>4.) {//high momentum electrons + hmpt->Fill((pt_1 - 1/ptg)/(1/ptg)*100.); + } else { + Float_t phig=TMath::ATan2(ref->Py(),ref->Px()); + hp->Fill((phi - phig)*1000.); + + Float_t lamg=TMath::ATan2(ref->Pz(),ptg); + hl->Fill((lam - lamg)*1000.); + + hpt->Fill((pt_1 - 1/ptg)/(1/ptg)*100.); + } + + Float_t mom=pt/TMath::Cos(lam); + Float_t dedx=track->GetTPCsignal(); + hep->Fill(mom,dedx,1.); + if (TMath::Abs(pdg)==211) //pions + if (mom>0.4 && mom<0.5) { + he->Fill(dedx,1.); + } } - Float_t mom=1./(pt_1*TMath::Cos(lam)); - Float_t dedx=track->GetdEdx(); - hep->Fill(mom,dedx,1.); - if (TMath::Abs(gt[ngood].code)==211) - if (mom>0.4 && mom<0.5) { - he->Fill(dedx,1.); - } + cout<<"\nList of Not found tracks :\n"; + for ( i = 0; i< itrack_notfound; i++){ + cout<GetEntries(); - if (ng!=0) - cerr<<"\n\n\nIntegral efficiency is about "<Clear(); + }// ***** End of the loop over events + delete event; + delete esdTree; + ef->Close(); - cout<<"\nList of Not found tracks :\n"; - // Int_t i; - for ( i = 0; i< itrack_notfound; i++){ - cout<Close(); + + Stat_t ng=hgood->GetEntries(), nf=hfound->GetEntries(); + if (ng!=0) cout<<"\n\nIntegral efficiency is about "<Draw(); p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10); - hp->SetFillColor(4); hp->SetXTitle("(mrad)"); hp->Fit("gaus"); c1->cd(); + hp->SetFillColor(4); hp->SetXTitle("(mrad)"); + if (hp->GetEntries()Draw(); else hp->Fit("gaus"); c1->cd(); TPad *p2=new TPad("p2","",0.5,.3,1,.6); p2->Draw(); p2->cd(); p2->SetFillColor(42); p2->SetFrameFillColor(10); - hl->SetXTitle("(mrad)"); hl->Fit("gaus"); c1->cd(); + hl->SetXTitle("(mrad)"); + if (hl->GetEntries()Draw(); else hl->Fit("gaus"); c1->cd(); TPad *p3=new TPad("p3","",0,0,0.5,0.3); p3->Draw(); p3->cd(); p3->SetFillColor(42); p3->SetFrameFillColor(10); - hpt->SetXTitle("(%)"); hpt->Fit("gaus"); c1->cd(); + hpt->SetXTitle("(%)"); + if (hpt->GetEntries()Draw(); else hpt->Fit("gaus"); c1->cd(); TPad *p4=new TPad("p4","",0.5,0,1,0.3); p4->Draw(); p4->cd(); p4->SetFillColor(42); p4->SetFrameFillColor(10); - hmpt->SetXTitle("(%)"); hmpt->Fit("gaus"); c1->cd(); + hmpt->SetXTitle("(%)"); + if (hmpt->GetEntries()Draw(); else hmpt->Fit("gaus"); c1->cd(); TPad *p5=new TPad("p5","",0,0.6,1,1); p5->Draw(); p5->cd(); p5->SetFillColor(41); p5->SetFrameFillColor(10); @@ -295,9 +330,9 @@ Int_t AliTPCComparison(Int_t eventn=1) { hg->SetXTitle("Pt (GeV/c)"); hg->Draw(); - TLine *line1 = new TLine(pmin,1.0,pmax,1.0); line1->SetLineStyle(4); + TLine *line1 = new TLine(0.1,1.0,6.1,1.0); line1->SetLineStyle(4); line1->Draw("same"); - TLine *line2 = new TLine(pmin,0.9,pmax,0.9); line2->SetLineStyle(4); + TLine *line2 = new TLine(0.1,0.9,6.1,0.9); line2->SetLineStyle(4); line2->Draw("same"); hf->SetFillColor(1); @@ -320,195 +355,268 @@ Int_t AliTPCComparison(Int_t eventn=1) { p6->cd(); p6->SetFillColor(42); p6->SetFrameFillColor(10); he->SetFillColor(2); he->SetFillStyle(3005); he->SetXTitle("Arbitrary Units"); - he->Fit("gaus"); c2->cd(); + if (he->GetEntries()Draw(); else 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(); + TFile fc("AliTPCComparison.root","RECREATE"); + c1->Write(); + c2->Write(); + fc.Close(); + gBenchmark->Stop("AliTPCComparison"); gBenchmark->Show("AliTPCComparison"); - return 0; } -Int_t good_tracks_tpc(GoodTrackTPC *gt, Int_t max, Int_t eventn) { - //eventn - number of events in file +Int_t GoodTracksTPC(const Char_t *dir) { + Char_t fname[100]; + sprintf(fname,"%s/galice.root",dir); - TFile *file=TFile::Open("galice.root"); - if (!file->IsOpen()) {cerr<<"Can't open galice.root !\n"; exit(4);} - - if (!(gAlice=(AliRun*)file->Get("gAlice"))) { - cerr<<"gAlice have not been found on galice.root !\n"; - exit(5); + AliRunLoader *rl = AliRunLoader::Open(fname,"COMPARISON"); + if (!rl) { + ::Error("GoodTracksTPC","Can't start session !"); + return 1; } - // Int_t np=gAlice->GetEvent(0); //MI change + rl->LoadgAlice(); + rl->LoadHeader(); + rl->LoadKinematics(); + rl->LoadTrackRefs(); - AliTPC *TPC=(AliTPC*)gAlice->GetDetector("TPC"); + AliTPCLoader *tpcl = (AliTPCLoader *)rl->GetLoader("TPCLoader"); + if (tpcl == 0x0) { + ::Error("GoodTracksTPC","Can not find TPCLoader !"); + delete rl; + return 2; + } + + AliTPC *TPC=(AliTPC*)rl->GetAliRun()->GetDetector("TPC"); Int_t ver = TPC->IsVersion(); - cerr<<"TPC version "<Get("75x40_100x60"); - if (!digp) { cerr<<"TPC parameters have not been found !\n"; exit(6); } - TPC->SetParam(digp); + cout<<"TPC version "<LoadRecPoints(); + else if (ver==2) tpcl->LoadDigits(); + else { + ::Error("GoodTracksTPC","Wrong TPC version !"); + delete rl; + return 3; + } + + + AliCDBManager *man=AliCDBManager::Instance(); + man->SetDefaultStorage("local://$ALICE_ROOT/OCDB"); + man->SetRun(0); + AliTPCParamSR *digp= + (AliTPCParamSR*)(AliTPCcalibDB::Instance()->GetParameters()); + if (!digp) { + ::Error("AliTPCComparison.C","TPC parameters have not been found !"); + delete rl; + return 4; + } Int_t nrow_up=digp->GetNRowUp(); Int_t nrows=digp->GetNRowLow()+nrow_up; Int_t zero=digp->GetZeroSup(); - Int_t gap=Int_t(0.125*nrows); + Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap); Int_t good_number=Int_t(0.4*nrows); + Int_t nev=rl->GetNumberOfEvents(); + ::Info("GoodTracksTPC","Number of events : %d\n",nev); + + + sprintf(fname,"%s/GoodTracksTPC.root",dir); + TFile *file=TFile::Open(fname,"recreate"); + + TClonesArray dummy("AliTrackReference",1000), *refs=&dummy; + TTree tpcTree("tpcTree","Tree with info about the reconstructable TPC tracks"); + tpcTree.Branch("TPC",&refs); + + //******** Loop over generated events + for (Int_t e=0; eClear(); - //MI change - Int_t nt=0; //reset counter - char treeName[100]; //declare event identifier - - for (Int_t event=0;eventGetEvent(event); - Int_t *good=new Int_t[np]; - for (Int_t ii=0; iiGetEvent(e); file->cd(); + + Int_t np = rl->GetHeader()->GetNtrack(); + cout<<"Event "<IsOpen()){cerr<<"Can't open AliTPCclusters.root !\n";exit(5);} - AliTPCClustersArray *ca=new AliTPCClustersArray; - ca->Setup(digp); - ca->SetClusterType("AliTPCcluster"); - //ca->ConnectTree("Segment Tree"); - ca->ConnectTree("TreeC_TPC_0"); - Int_t nrows=Int_t(ca->GetTree()->GetEntries()); - for (Int_t n=0; nLoadEntry(n); - Int_t sec,row; - digp->AdjustSectorRow(s->GetID(),sec,row); - AliTPCClustersRow &clrow = *ca->GetRow(sec,row); - Int_t ncl=clrow.GetArray()->GetEntriesFast(); - while (ncl--) { - AliTPCcluster *c=(AliTPCcluster*)clrow[ncl]; - Int_t lab=c->GetLabel(0); - if (lab<0) continue; //noise cluster - lab=TMath::Abs(lab); - if (sec>=digp->GetNInnerSector()) - if (row==nrow_up-1 ) good[lab]|=0x1000; - if (sec>=digp->GetNInnerSector()) - if (row==nrow_up-1-gap) good[lab]|=0x800; - good[lab]++; - } - ca->ClearRow(sec,row); - } - cf->Close(); - } - break; + { + AliTPCClustersArray *pca=new AliTPCClustersArray, &ca=*pca; + ca.Setup(digp); + ca.SetClusterType("AliTPCcluster"); + ca.ConnectTree(tpcl->TreeR()); + Int_t nrows=Int_t(ca.GetTree()->GetEntries()); + for (Int_t n=0; nAdjustSectorRow(s->GetID(),sec,row); + AliTPCClustersRow &clrow = *ca.GetRow(sec,row); + Int_t ncl=clrow.GetArray()->GetEntriesFast(); + while (ncl--) { + AliTPCcluster *c=(AliTPCcluster*)clrow[ncl]; + Int_t lab=c->GetLabel(0); + if (lab<0) continue; //noise cluster + lab=TMath::Abs(lab); + + if (sec>=digp->GetNInnerSector()) + if (row==nrow_up-1) good[lab]|=0x4000; + if (sec>=digp->GetNInnerSector()) + if (row==nrow_up-1-gap) good[lab]|=0x1000; + + if (sec>=digp->GetNInnerSector()) + if (row==nrow_up-1-shift) good[lab]|=0x2000; + if (sec>=digp->GetNInnerSector()) + if (row==nrow_up-1-gap-shift) good[lab]|=0x800; + + good[lab]++; + } + ca.ClearRow(sec,row); + } + delete pca; + } + break; case 2: - - sprintf(treeName,"TreeD_75x40_100x60_%d",event); - TD=(TTree*)gDirectory->Get(treeName); - TD->GetBranch("Segment")->SetAddress(&digits); - count = new Int_t[np]; - Int_t i; - for (i=0; iGetEntries(); - for (i=0; iGetEvent(i)) continue; - Int_t sec,row; - digp->AdjustSectorRow(digits->GetID(),sec,row); - cerr<First(); - do { //Many thanks to J.Chudoba who noticed this - Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn(); - Short_t dig = digits->GetDigit(it,ip); - Int_t idx0=digits->GetTrackID(it,ip,0); - Int_t idx1=digits->GetTrackID(it,ip,1); - Int_t idx2=digits->GetTrackID(it,ip,2); - if (idx0>=0 && dig>=zero) count[idx0]+=1; - if (idx1>=0 && dig>=zero) count[idx1]+=1; - if (idx2>=0 && dig>=zero) count[idx2]+=1; + { + TTree *TD=tpcl->TreeD(); + + AliSimDigits da, *digits=&da; + TD->GetBranch("Segment")->SetAddress(&digits); + + Int_t *count = new Int_t[np]; + Int_t i; + for (i=0; iGetEntries(); + for (i=0; iGetEvent(i)) continue; + Int_t sec,row; + digp->AdjustSectorRow(digits->GetID(),sec,row); + //cerr<First(); + do { //Many thanks to J.Chudoba who noticed this + Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn(); + Short_t dig = digits->GetDigit(it,ip); + Int_t idx0=digits->GetTrackID(it,ip,0); + Int_t idx1=digits->GetTrackID(it,ip,1); + Int_t idx2=digits->GetTrackID(it,ip,2); + if (idx0>=0 && dig>=zero && idx0=0 && dig>=zero && idx1=0 && dig>=zero && idx2Next()); - for (Int_t j=0; j1) { - if (sec>=digp->GetNInnerSector()) - if (row==nrow_up-1 ) good[j]|=0x1000; - if (sec>=digp->GetNInnerSector()) - if (row==nrow_up-1-gap) good[j]|=0x800; - good[j]++; - } - count[j]=0; - } - } - delete[] count; - delete TD; //Thanks to Mariana Bondila - break; - default: - cerr<<"Invalid TPC version !\n"; - file->Close(); - exit(7); + for (Int_t j=0; j1) { + if (sec>=digp->GetNInnerSector()) + if (row==nrow_up-1 ) good[j]|=0x4000; + if (sec>=digp->GetNInnerSector()) + if (row==nrow_up-1-gap) good[j]|=0x1000; + + if (sec>=digp->GetNInnerSector()) + if (row==nrow_up-1-shift) good[j]|=0x2000; + if (sec>=digp->GetNInnerSector()) + if (row==nrow_up-1-gap-shift) good[j]|=0x800; + good[j]++; + } + count[j]=0; + } + } + delete[] count; + } + break; } - - TTree *TH=gAlice->TreeH(); - Int_t npart=(Int_t)TH->GetEntries(); - - while (npart--) { - AliTPChit *hit0=0; - - TPC->ResetHits(); - TH->GetEvent(npart); - AliTPChit * hit = (AliTPChit*) TPC->FirstHit(-1); - while (hit){ - if (hit->fQ==0.) break; - hit = (AliTPChit*) TPC->NextHit(); - } - if (hit) { - hit0 = new AliTPChit(*hit); //Make copy of hit - hit = hit0; - } - else continue; - AliTPChit *hit1=(AliTPChit*)TPC->NextHit(); - if (hit1==0) continue; - if (hit1->fQ != 0.) continue; - Int_t i=hit->Track(); - TParticle *p = (TParticle*)gAlice->Particle(i); - - if (p->GetFirstMother()>=0) continue; //secondary particle - if (good[i] < 0x1000+0x800+2+good_number) continue; - if (p->Pt()<0.100) continue; - if (TMath::Abs(p->Pz()/p->Pt())>0.999) continue; - - gt[nt].lab=i; - gt[nt].code=p->GetPdgCode(); - //**** px py pz - in global coordinate system, x y z - in local ! - gt[nt].px=hit->X(); gt[nt].py=hit->Y(); gt[nt].pz=hit->Z(); - Float_t cs,sn; digp->AdjustCosSin(hit1->fSector,cs,sn); - gt[nt].fEventN=event; //MI change - gt[nt].x = hit1->X()*cs + hit1->Y()*sn; - gt[nt].y =-hit1->X()*sn + hit1->Y()*cs; - gt[nt].z = hit1->Z(); - nt++; - if (hit0) delete hit0; - cerr<Stack(); + for (i=0; iParticle(i); + if (p == 0x0) { + cerr<<"Can not get particle "<Pt()<0.100) continue; + if (TMath::Abs(p->Pz()/p->Pt())>0.999) continue; + + Double_t vx=p->Vx(),vy=p->Vy(),vz=p->Vz(); + if (TMath::Sqrt(vx*vx+vy*vy)>3.5) continue; + if (TMath::Abs(vz) > 50.) continue; + + AliTrackReference *ref=new((*refs)[nt]) AliTrackReference(); + + ref->SetLabel(i); + Int_t pdg=p->GetPdgCode(); + ref->SetLength(pdg); //This will the particle's PDG ! + ref->SetMomentum(0.,0.,0.); + ref->SetPosition(0.,0.,0.); + + nt++; } - delete[] good; - } - delete gAlice; gAlice=0; - file->Close(); - gBenchmark->Stop("AliTPCComparison"); - gBenchmark->Show("AliTPCComparison"); - return nt; -} + //**** check if there is also information at the entrance of the TPC + TTree *TR=rl->TreeTR(); + TBranch *branch=TR->GetBranch("TrackReferences"); + if (branch==0) { + ::Error("GoodTracksTPC","No track references !"); + delete rl; + return 5; + } + TClonesArray tpcdummy("AliTrackReference",1000), *tpcRefs=&tpcdummy; + branch->SetAddress(&tpcRefs); + Int_t nr=(Int_t)TR->GetEntries(); + for (Int_t r=0; rGetEvent(r); + + Int_t nref = tpcRefs->GetEntriesFast(); + if (!nref) continue; + AliTrackReference *tpcRef= 0x0; + for (Int_t iref=0; irefUncheckedAt(iref); + if (tpcRef->DetectorId() == AliTrackReference::kTPC) break; + tpcRef = 0x0; + } + + if (!tpcRef) continue; + Int_t j; + AliTrackReference *ref=0; + for (j=0; jUncheckedAt(j); + if (ref->Label()==tpcRef->Label()) break; + ref=0; + } + if (ref==0) continue; + ref->SetMomentum(tpcRef->Px(),tpcRef->Py(),tpcRef->Pz()); + ref->SetPosition(tpcRef->LocalX(),tpcRef->LocalY(),tpcRef->Z()); + + tpcRefs->Clear(); + } + + tpcTree.Fill(); + + delete[] good; + } //****** end of the loop over generated events + + + tpcTree.Write(); + file->Close(); + + delete rl; + return 0; +}