From ff2225f6c66905c5f0f9160da67a229be816826b Mon Sep 17 00:00:00 2001 From: belikov Date: Mon, 14 Jul 2014 11:18:25 +0200 Subject: [PATCH] Starting a collection of QA/Comparison macros --- ITS/UPGRADE/macros/QA/AliITSUComparison.C | 540 +++++++++++++++ .../macros/QA/AliITSUComparisonCooked.C | 555 ++++++++++++++++ ITS/UPGRADE/macros/QA/AliTPCUComparison.C | 626 ++++++++++++++++++ ITS/UPGRADE/macros/QA/check_radii.C | 141 ++++ ITS/UPGRADE/macros/QA/compClusHits.C | 541 +++++++++++++++ ITS/UPGRADE/macros/QA/plotHits.C | 188 ++++++ 6 files changed, 2591 insertions(+) create mode 100644 ITS/UPGRADE/macros/QA/AliITSUComparison.C create mode 100644 ITS/UPGRADE/macros/QA/AliITSUComparisonCooked.C create mode 100644 ITS/UPGRADE/macros/QA/AliTPCUComparison.C create mode 100644 ITS/UPGRADE/macros/QA/check_radii.C create mode 100644 ITS/UPGRADE/macros/QA/compClusHits.C create mode 100644 ITS/UPGRADE/macros/QA/plotHits.C diff --git a/ITS/UPGRADE/macros/QA/AliITSUComparison.C b/ITS/UPGRADE/macros/QA/AliITSUComparison.C new file mode 100644 index 00000000000..2a734df48c5 --- /dev/null +++ b/ITS/UPGRADE/macros/QA/AliITSUComparison.C @@ -0,0 +1,540 @@ +/**************************************************************************** + * Legacy comparison macro, adapted to the upgraded ITS. * + * * + * Creates list of "trackable" tracks, * + * calculates efficiency, resolutions etc. * + * The ESD tracks must be in an appropriate state (kITSin or kITSrefit) * + * * + * Before running, load the ITSU libraries: * + * gSystem->Load("libITSUpgradeBase");gSystem->Load("libITSUpgradeRec"); * + * * + * Origin: I.Belikov, IPHC, Iouri.Belikov@iphc.cnrs.fr * + ****************************************************************************/ + +#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 "UPGRADE/AliITSUClusterPix.h" + #include "UPGRADE/AliITSULoader.h" +#endif + +Int_t GoodTracksITS(const Char_t *dir="."); + +extern AliRun *gAlice; +extern TBenchmark *gBenchmark; +extern TROOT *gROOT; + +static Int_t allgood=0; +static Int_t allselected=0; +static Int_t allfound=0; + +Int_t AliITSUComparison +(Float_t ptcutl=0., Float_t ptcuth=2., const Char_t *dir=".") { + gBenchmark->Start("AliITSUComparison"); + + ::Info("AliITSUComparison.C","Doing comparison..."); + + + 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","Transverse impact parameter",30,-777,777); + hmpt->SetFillColor(6); + + TH1F *hz=(TH1F*)gROOT->FindObject("hz"); + if (!hz) hz=new TH1F("hz","Longitudinal impact parameter",30,-777,777); + + + Int_t nb=100; + TH1F *hgood=(TH1F*)gROOT->FindObject("hgood"); + if (!hgood) hgood=new TH1F("hgood","Good tracks",nb,ptcutl,ptcuth); + + TH1F *hfound=(TH1F*)gROOT->FindObject("hfound"); + if (!hfound) hfound=new TH1F("hfound","Found tracks",nb,ptcutl,ptcuth); + + TH1F *hfake=(TH1F*)gROOT->FindObject("hfake"); + if (!hfake) hfake=new TH1F("hfake","Fake tracks",nb,ptcutl,ptcuth); + + TH1F *hg=(TH1F*)gROOT->FindObject("hg"); + if (!hg) hg=new TH1F("hg","Efficiency for good tracks",nb,ptcutl,ptcuth); + hg->SetLineColor(4); hg->SetLineWidth(2); + + TH1F *hf=(TH1F*)gROOT->FindObject("hf"); + if (!hf) hf=new TH1F("hf","Efficiency for fake tracks",nb,ptcutl,ptcuth); + hf->SetFillColor(1); hf->SetFillStyle(3013); hf->SetLineWidth(2); + + + TH1F *he=(TH1F*)gROOT->FindObject("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); + + + Char_t fname[100]; + sprintf(fname,"%s/GoodTracksITS.root",dir); + + TFile *refFile=TFile::Open(fname,"old"); + if (!refFile || !refFile->IsOpen()) { + ::Info("AliITSUComparison.C","Marking good tracks (will take a while)..."); + if (GoodTracksITS(dir)) { + ::Error("AliITSUComparison.C","Can't generate the reference file !"); + return 1; + } + } + refFile=TFile::Open(fname,"old"); + if (!refFile || !refFile->IsOpen()) { + ::Error("AliITSUComparison.C","Can't open the reference file !"); + return 1; + } + + TTree *itsTree=(TTree*)refFile->Get("itsTree"); + if (!itsTree) { + ::Error("AliITSUComparison.C","Can't get the reference tree !"); + return 2; + } + TBranch *branch=itsTree->GetBranch("ITS"); + if (!branch) { + ::Error("AliITSUComparison.C","Can't get the ITS branch !"); + return 3; + } + 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/AliESDits.root",dir); + ef=TFile::Open(fname); + if ((!ef)||(!ef->IsOpen())) { + ::Error("AliITSUComparison.C","Can't open AliESDits.root !"); + return 4; + } + } + AliESDEvent* event = new AliESDEvent(); + TTree* esdTree = (TTree*) ef->Get("esdTree"); + if (!esdTree) { + ::Error("AliITSComparison.C", "no ESD tree found"); + return 6; + } + event->ReadFromTree(esdTree); + + + //******* Loop over events ********* + Int_t e=0; + while (esdTree->GetEvent(e)) { + cout<GetNumberOfTracks(); + allfound+=nentr; + + if (itsTree->GetEvent(e++)==0) { + cerr<<"No reconstructable tracks !\n"; + continue; + } + + Int_t ngood=refs->GetEntriesFast(); + allgood+=ngood; + + const Int_t MAX=15000; + Int_t notf[MAX], nnotf=0; + Int_t fake[MAX], nfake=0; + Int_t mult[MAX], numb[MAX], nmult=0; + Int_t k; + for (k=0; kUncheckedAt(k); + Int_t lab=ref->Label(), tlab=-1; + Float_t ptg=TMath::Sqrt(ref->Px()*ref->Px() + ref->Py()*ref->Py()); + + if (ptgptcuth) continue; + + allselected++; + + hgood->Fill(ptg); + + AliESDtrack *esd=0; + Int_t cnt=0; + for (Int_t i=0; iGetTrack(i); + UInt_t status=t->GetStatus(); + + if ((status&AliESDtrack::kITSrefit)==0) continue; + + Int_t lbl=t->GetLabel(); + if (lab==TMath::Abs(lbl)) { + if (cnt==0) {esd=t; tlab=lbl;} + if (lbl> 0) {esd=t; tlab=lbl;} + cnt++; + } + } + if (cnt==0) { + notf[nnotf++]=lab; + continue; + } else if (cnt>1){ + mult[nmult]=lab; + numb[nmult]=cnt; nmult++; + } + + if (lab==tlab) hfound->Fill(ptg); + else { + fake[nfake++]=lab; + hfake->Fill(ptg); + } + + Double_t alpha=esd->GetAlpha(),xv,par[5]; + esd->GetExternalParameters(xv,par); + Float_t phi=TMath::ASin(par[2]) + alpha; + 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]); + + 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.); + + Float_t d,z; esd->GetImpactParameters(d,z); + hmpt->Fill(10000*d); + hz->Fill(10000*z); + + hpt->Fill((pt_1 - 1/ptg)/(1/ptg)*100.); + + Float_t mom=(TMath::Abs(pt_1)>1e-5) ? 1./(pt_1*TMath::Cos(lam)) : 1e+5; + Float_t dedx=esd->GetITSsignal(); + hep->Fill(mom,dedx,1.); + + Int_t pdg=(Int_t)ref->GetLength(); //this is particle's PDG ! + + if (TMath::Abs(pdg)==211) //pions + if (mom>0.4 && mom<0.5) he->Fill(dedx,1.); + + } + + cout<<"\nList of Not found tracks :\n"; + for (k=0; kClear(); + } //***** End of the loop over events + + delete event; + delete esdTree; + ef->Close(); + + delete itsTree; + refFile->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)"); + 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)"); + 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("(%)"); + 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("(micron)"); + if (hmpt->GetEntries()Draw(); else hmpt->Fit("gaus"); + hz->Draw("same"); c1->cd(); + + + TPad *p5=new TPad("p5","",0,0.6,1,1); p5->Draw(); p5->cd(); + p5->SetFillColor(41); p5->SetFrameFillColor(10); + hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2(); + hg->Divide(hfound,hgood,1,1.,"b"); + hf->Divide(hfake,hgood,1,1.,"b"); + hg->SetMaximum(1.4); + hg->SetYTitle("Tracking efficiency"); + hg->SetXTitle("Pt (GeV/c)"); + hg->Draw(); + + TLine *line1 = new TLine(ptcutl,1.0,ptcuth,1.0); line1->SetLineStyle(4); + line1->Draw("same"); + TLine *line2 = new TLine(ptcutl,0.9,ptcuth,0.9); line2->SetLineStyle(4); + line2->Draw("same"); + + hf->SetFillColor(1); + hf->SetFillStyle(3013); + hf->SetLineColor(2); + hf->SetLineWidth(2); + hf->Draw("histsame"); + TText *text = new TText(0.461176,0.248448,"Fake tracks"); + text->SetTextSize(0.05); + text->Draw(); + text = new TText(0.453919,1.11408,"Good tracks"); + text->SetTextSize(0.05); + text->Draw(); + + 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"); + 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("AliITSUComparison.root","RECREATE"); + c1->Write(); + c2->Write(); + fc.Close(); + + gBenchmark->Stop("AliITSUComparison"); + gBenchmark->Show("AliITSUComparison"); + + return 0; +} + + + +Int_t GoodTracksITS(const Char_t *dir) { + if (gAlice) { + delete AliRunLoader::Instance(); + delete gAlice;//if everything was OK here it is already NULL + gAlice = 0x0; + } + + Char_t fname[100]; + sprintf(fname,"%s/galice.root",dir); + + AliRunLoader *rl = AliRunLoader::Open(fname,"COMPARISON"); + if (!rl) { + ::Error("GoodTracksITS","Can't start session !"); + return 1; + } + + rl->LoadgAlice(); + rl->LoadHeader(); + rl->LoadKinematics(); + + AliITSULoader* itsl = (AliITSULoader*)rl->GetLoader("ITSLoader"); + if (itsl == 0x0) { + ::Error("GoodTracksITS","Can not find the ITSLoader"); + delete rl; + return 4; + } + itsl->LoadRecPoints(); + + + Int_t nev=rl->GetNumberOfEvents(); + ::Info("GoodTracksITS","Number of events : %d\n",nev); + + sprintf(fname,"%s/GoodTracksTPC.root",dir); + TFile *tpcFile=TFile::Open(fname); + if ((!tpcFile)||(!tpcFile->IsOpen())) { + ::Error("GoodTracksITS","Can't open the GoodTracksTPC.root !"); + delete rl; + return 5; + } + TClonesArray dum("AliTrackReference",1000), *tpcRefs=&dum; + TTree *tpcTree=(TTree*)tpcFile->Get("tpcTree"); + if (!tpcTree) { + ::Error("GoodTracksITS","Can't get the TPC reference tree !"); + delete rl; + return 6; + } + TBranch *tpcBranch=tpcTree->GetBranch("TPC"); + if (!tpcBranch) { + ::Error("GoodTracksITS","Can't get the TPC reference branch !"); + delete rl; + return 7; + } + tpcBranch->SetAddress(&tpcRefs); + + sprintf(fname,"%s/GoodTracksITS.root",dir); + TFile *itsFile=TFile::Open(fname,"recreate"); + TClonesArray dummy2("AliTrackReference",1000), *itsRefs=&dummy2; + TTree itsTree("itsTree","Tree with info about the reconstructable ITS tracks"); + itsTree.Branch("ITS",&itsRefs); + + //******** Loop over generated events + for (Int_t e=0; eGetEvent(e); itsFile->cd(); + + Int_t np = rl->GetHeader()->GetNtrack(); + cout<<"Event "<TreeR(); + if (!cTree) { + ::Error("GoodTracksITS","Can't get the cluster tree !"); + delete rl; + return 8; + } + + const Int_t nLayers=7; + TBranch *branch[nLayers]; + TClonesArray clusters[nLayers]; + for (Int_t layer=0; layerGetBranch(bname); + if (!branch[layer]) { + ::Error("GoodTracksITS","Can't get the clusters branch !"); + delete rl; + return 9; + } + branch[layer]->SetAddress(&ptr); + } + + Int_t entr=(Int_t)cTree->GetEntries(); + for (k=0; kGetEvent(k); + for (Int_t lay=0; layGetLabel(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; + if (l2>=0) good[l2]|=mask; + } + clusters[lay].Clear(); + } + } + + + + //****** select tracks which are "good" enough + AliStack* stack = rl->Stack(); + + tpcTree->GetEvent(e); + Int_t nk=tpcRefs->GetEntriesFast(); + Int_t nt=0; + for (k=0; kUncheckedAt(k); + Int_t lab=tpcRef->Label(); + if (good[lab] != 0x7F) continue; + TParticle *p = (TParticle*)stack->Particle(lab); + if (p == 0x0) { + cerr<<"Can not get particle "<SetMomentum(p->Px(),p->Py(),p->Pz()); + ref->SetPosition(p->Vx(),p->Vy(),p->Vz()); + nt++; + } + tpcRefs->Clear(); + + itsTree.Fill(); + itsRefs->Clear(); + + delete[] good; + + } //*** end of the loop over generated events + + itsTree.Write(); + itsFile->Close(); + + delete tpcTree; + tpcFile->Close(); + + delete rl; + return 0; +} + + diff --git a/ITS/UPGRADE/macros/QA/AliITSUComparisonCooked.C b/ITS/UPGRADE/macros/QA/AliITSUComparisonCooked.C new file mode 100644 index 00000000000..3aec0d468cf --- /dev/null +++ b/ITS/UPGRADE/macros/QA/AliITSUComparisonCooked.C @@ -0,0 +1,555 @@ +/**************************************************************************** + * A standalone comparison macro for the upgraded ITS. * + * * + * Creates list of "trackable" tracks, * + * calculates efficiency, resolutions etc. * + * The ESD tracks must be in an appropriate state: kITSrefit * + * * + * The efficiency and resolutions are calculated * + * wrt "primary-like" pions within the acceptance tan(lambda)~[-1,1] * + * and for tracks having all 7 clusters correctly assigned. * + * * + * Before running, load the ITSU libraries: * + * gSystem->Load("libITSUpgradeBase");gSystem->Load("libITSUpgradeRec"); * + * * + * Origin: I.Belikov, IPHC, Iouri.Belikov@iphc.cnrs.fr * + ****************************************************************************/ + +#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 "AliGenEventHeader.h" + #include "AliTrackReference.h" + #include "AliRunLoader.h" + #include "AliRun.h" + #include "AliESDEvent.h" + #include "AliESDtrack.h" + + #include "UPGRADE/AliITSUClusterPix.h" + #include "UPGRADE/AliITSULoader.h" +#endif + +Int_t GoodTracksCooked(const Char_t *dir="."); + +extern AliRun *gAlice; +extern TBenchmark *gBenchmark; +extern TROOT *gROOT; + +static Int_t allgood=0; +static Int_t allselected=0; +static Int_t allfound=0; + +void root(TH1 *h) { + Int_t nb=h->GetNbinsX(); + for (Int_t i=0; iGetBinContent(i); + c=TMath::Sqrt(c); + h->SetBinContent(i,c); + + Float_t e=h->GetBinError(i); + if (c!=0) e = 0.5*e/c; + h->SetBinError(i,e); + } +} +void divide(TH1 *h) { + Int_t nb=h->GetNbinsX(); + for (Int_t i=0; iGetBinContent(i); + c *= h->GetBinCenter(i); + h->SetBinContent(i,c); + + Float_t e=h->GetBinError(i); + e *= h->GetBinCenter(i); + h->SetBinError(i,e); + } +} + +Int_t AliITSUComparisonCooked +(Float_t ptcutl=0.01, Float_t ptcuth=10., const Char_t *dir=".") { + gBenchmark->Start("AliITSUComparisonCooked"); + + ::Info("AliITSUComparisonCooked.C","Doing comparison..."); + + + 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 *hz=(TH1F*)gROOT->FindObject("hz"); + //if (!hz) hz=new TH1F("hz","Longitudinal impact parameter",30,-777.,777.); + + + Int_t nb=100; + Float_t xbins[nb+1]; + Double_t a=TMath::Log(ptcuth/ptcutl)/nb; + for (Int_t i=0; i<=nb; i++) xbins[i] = ptcutl*TMath::Exp(i*a); + + TH1F *hgood=(TH1F*)gROOT->FindObject("hgood"); + if (!hgood) hgood=new TH1F("hgood","Good tracks",nb,xbins); + + TH1F *hfound=(TH1F*)gROOT->FindObject("hfound"); + if (!hfound) hfound=new TH1F("hfound","Found tracks",nb,xbins); + + TH1F *hfake=(TH1F*)gROOT->FindObject("hfake"); + if (!hfake) hfake=new TH1F("hfake","Fake tracks",nb,xbins); + + TH1F *hg=(TH1F*)gROOT->FindObject("hg"); + if (!hg) hg=new TH1F("hg","Efficiency for good tracks",nb,xbins); + hg->SetLineColor(4); hg->SetLineWidth(2); + + TH1F *hf=(TH1F*)gROOT->FindObject("hf"); + if (!hf) hf=new TH1F("hf","Efficiency for fake tracks",nb,xbins); + hf->SetFillColor(1); hf->SetFillStyle(3013); hf->SetLineWidth(2); + + TH1F *hpt=(TH1F*)gROOT->FindObject("hpt"); + if (!hpt) hpt=new TH1F("hpt","Relative Pt resolution",nb,xbins); + hpt->Sumw2(); + TH1F *hd=(TH1F*)gROOT->FindObject("hd"); + if (!hd) + hd=new TH1F("hd","Transverse impact parameter",nb,xbins); + hd->Sumw2(); + + + Char_t fname[100]; + sprintf(fname,"%s/GoodTracksCooked.root",dir); + + TFile *refFile=TFile::Open(fname,"old"); + if (!refFile || !refFile->IsOpen()) { + ::Info("AliITSUComparisonCooked.C","Marking good tracks (will take a while)..."); + if (GoodTracksCooked(dir)) { + ::Error("AliITSUComparisonCooked.C","Can't generate the reference file !"); + return 1; + } + } + refFile=TFile::Open(fname,"old"); + if (!refFile || !refFile->IsOpen()) { + ::Error("AliITSUComparisonCooked.C","Can't open the reference file !"); + return 1; + } + + TTree *itsTree=(TTree*)refFile->Get("itsTree"); + if (!itsTree) { + ::Error("AliITSUComparisonCooked.C","Can't get the reference tree !"); + return 2; + } + TBranch *branch=itsTree->GetBranch("ITS"); + if (!branch) { + ::Error("AliITSUComparisonCooked.C","Can't get the ITS branch !"); + return 3; + } + 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/AliESDits.root",dir); + ef=TFile::Open(fname); + if ((!ef)||(!ef->IsOpen())) { + ::Error("AliITSUComparisonCooked.C","Can't open AliESDits.root !"); + return 4; + } + } + AliESDEvent* event = new AliESDEvent(); + TTree* esdTree = (TTree*) ef->Get("esdTree"); + if (!esdTree) { + ::Error("AliITSComparison.C", "no ESD tree found"); + return 6; + } + event->ReadFromTree(esdTree); + + + //******* Loop over events ********* + Int_t e=0; + while (esdTree->GetEvent(e)) { + cout<GetNumberOfTracks(); + allfound+=nentr; + + if (itsTree->GetEvent(e++)==0) { + cerr<<"No reconstructable tracks !\n"; + continue; + } + + Int_t ngood=refs->GetEntriesFast(); + allgood+=ngood; + + const Int_t MAX=15000; + Int_t notf[MAX], nnotf=0; + Int_t fake[MAX], nfake=0; + Int_t mult[MAX], numb[MAX], nmult=0; + Int_t k; + for (k=0; kUncheckedAt(k); + Int_t lab=ref->Label(), tlab=-1; + Float_t ptg=TMath::Sqrt(ref->Px()*ref->Px() + ref->Py()*ref->Py()); + + Int_t pdg=(Int_t)ref->GetLength(); //this is particle's PDG ! + if (TMath::Abs(pdg)!=211) continue; //select pions only + + if (ptgptcuth) continue; + + allselected++; + + hgood->Fill(ptg); + + AliESDtrack *esd=0; + Int_t cnt=0; + for (Int_t i=0; iGetTrack(i); + UInt_t status=t->GetStatus(); + + if ((status&AliESDtrack::kITSrefit)==0) continue; + if (t->GetITSclusters(0)<4) continue; + + Int_t lbl=t->GetLabel(); + if (lab==TMath::Abs(lbl)) { + if (cnt==0) {esd=t; tlab=lbl;} + if (lbl> 0) {esd=t; tlab=lbl;} + cnt++; + } + } + if (cnt==0) { + notf[nnotf++]=lab; + continue; + } else if (cnt>1){ + mult[nmult]=lab; + numb[nmult]=cnt; nmult++; + } + + if (lab==tlab) hfound->Fill(ptg); + else { + fake[nfake++]=lab; + hfake->Fill(ptg); + } + + if (esd->GetLabel()<0) continue; //resolutions for good tracks only + + Double_t alpha=esd->GetAlpha(),xv,par[5]; + esd->GetExternalParameters(xv,par); + Float_t phi=TMath::ASin(par[2]) + alpha; + 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]); + + 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.); + + Float_t d,z; esd->GetImpactParameters(d,z); + d*=10000; //microns + Float_t w=d*d; + hd->Fill(ptg, w); + //z*=10000; //microns + //hz->Fill(z); + + w=(pt_1 - 1/ptg)*100 * (pt_1-1/ptg)*100; + hpt->Fill(ptg, w); + + } + + cout<<"\nList of Not found tracks :\n"; + for (k=0; kClear(); + } //***** End of the loop over events + + delete event; + delete esdTree; + ef->Close(); + + delete itsTree; + refFile->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)"); + 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)"); + if (hl->GetEntries()Draw(); else hl->Fit("gaus"); c1->cd(); + + TPad *p3=new TPad("p3","",0,0,0.5,0.3); + p3->SetLogx(); p3->SetGridx(); p3->SetGridy(); + p3->Draw(); + p3->cd(); p3->SetFillColor(42); p3->SetFrameFillColor(10); + hpt->SetXTitle("p_{T} (GeV/c)"); + hpt->SetYTitle("(%)"); + TH1F *hh=new TH1F(*hfound); + //hh->Add(hfake); + hpt->Divide(hh); + root(hpt); + divide(hpt); + hpt->Draw(); c1->cd(); + + TPad *p4=new TPad("p4","",0.5,0,1,0.3); + p4->SetLogx(); p4->SetGridx(); p4->SetGridy(); + p4->Draw(); + p4->cd(); p4->SetFillColor(42); p4->SetFrameFillColor(10); + hd->SetXTitle("p_{T} (GeV/c)"); + hd->SetYTitle("(micron)"); + hd->Divide(hh); + root(hd); + hd->Draw(); c1->cd(); + + //hz->Draw("same"); c1->cd(); + + + TPad *p5=new TPad("p5","",0,0.6,1,1); + p5->SetLogx(); p5->SetGridx(); p5->SetGridy(); + p5->Draw(); p5->cd(); + p5->SetFillColor(41); p5->SetFrameFillColor(10); + hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2(); + hg->Divide(hfound,hgood,1,1.,"b"); + hf->Divide(hfake,hgood,1,1.,"b"); + hg->SetYTitle("Tracking efficiency (%)"); + hg->SetXTitle("p_{T} (GeV/c)"); + hg->Scale(100); + hg->Draw(); + + hf->SetFillColor(1); + hf->SetFillStyle(3013); + hf->SetLineColor(2); + hf->SetLineWidth(2); + hf->Scale(100); + hf->Draw("histsame"); + TText *text = new TText(0.4, 20., "Fake tracks"); + text->SetTextSize(0.05); + text->Draw(); + text = new TText(0.4, 80., "Good tracks"); + text->SetTextSize(0.05); + text->Draw(); + + TFile fc("AliITSUComparisonCooked.root","RECREATE"); + c1->Write(); + fc.Close(); + + gBenchmark->Stop("AliITSUComparisonCooked"); + gBenchmark->Show("AliITSUComparisonCooked"); + + return 0; +} + + + +Int_t GoodTracksCooked(const Char_t *dir) { + if (gAlice) { + delete AliRunLoader::Instance(); + delete gAlice;//if everything was OK here it is already NULL + gAlice = 0x0; + } + + Char_t fname[100]; + sprintf(fname,"%s/galice.root",dir); + + AliRunLoader *rl = AliRunLoader::Open(fname,"COMPARISON"); + if (!rl) { + ::Error("GoodTracksCooked","Can't start session !"); + return 1; + } + + rl->LoadgAlice(); + rl->LoadHeader(); + rl->LoadKinematics(); + + AliITSULoader* itsl = (AliITSULoader*)rl->GetLoader("ITSLoader"); + if (itsl == 0x0) { + ::Error("GoodTracksCooked","Can not find the ITSLoader"); + delete rl; + return 4; + } + itsl->LoadRecPoints(); + + + Int_t nev=rl->GetNumberOfEvents(); + ::Info("GoodTracksCooked","Number of events : %d\n",nev); + + sprintf(fname,"%s/GoodTracksCooked.root",dir); + TFile *itsFile=TFile::Open(fname,"recreate"); + TClonesArray dummy2("AliTrackReference",1000), *itsRefs=&dummy2; + TTree itsTree("itsTree","Tree with info about the reconstructable ITS tracks"); + itsTree.Branch("ITS",&itsRefs); + + //******** Loop over generated events + for (Int_t e=0; eGetEvent(e); itsFile->cd(); + + Int_t np = rl->GetHeader()->GetNtrack(); + cout<<"Event "<GetHeader()->GenEventHeader(); + TArrayF vtx(3); + h->PrimaryVertex(vtx); + + Bool_t skip=kFALSE; + if (TMath::Abs(vtx[2]) > 10.) { + cout<<"Skipping an event with Zv="<TreeR(); + if (!cTree) { + ::Error("GoodTracksCooked","Can't get the cluster tree !"); + delete rl; + return 8; + } + + const Int_t nLayers=7; + TBranch *branch[nLayers]; + TClonesArray clusters[nLayers]; + for (Int_t layer=0; layerGetBranch(bname); + if (!branch[layer]) { + ::Error("GoodTracksCooked","Can't get the clusters branch !"); + delete rl; + return 9; + } + branch[layer]->SetAddress(&ptr); + } + + Int_t entr=(Int_t)cTree->GetEntries(); + for (k=0; kGetEvent(k); + for (Int_t lay=0; layGetLabel(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; + if (l2>=0) good[l2]|=mask; + } + clusters[lay].Clear(); + } + } + + + + //****** select tracks which are "good" enough + AliStack* stack = rl->Stack(); + Int_t nt=0; + for (k=0; kParticle(k); + if (p == 0x0) { + cerr<<"Can not get particle "<Pt() <= 0.) continue; + if (TMath::Abs(p->Pz()/p->Pt())>0.999) continue; + + Double_t dx=p->Vx()-vtx[0], dy=p->Vy()-vtx[1], dz=p->Vz()-vtx[2]; + if (TMath::Sqrt(dx*dx+dy*dy)>0.0001) continue; //Primary-like + if (TMath::Abs(dz) > 0.0001) continue; + + AliTrackReference *ref=new((*itsRefs)[nt]) AliTrackReference(); + ref->SetLabel(k); + Int_t pdg=p->GetPdgCode(); + ref->SetLength(pdg); //This will the particle's PDG ! + ref->SetMomentum(p->Px(),p->Py(),p->Pz()); + ref->SetPosition(p->Vx(),p->Vy(),p->Vz()); + nt++; + } + + itsTree.Fill(); + itsRefs->Clear(); + + delete[] good; + + } //*** end of the loop over generated events + + itsTree.Write(); + itsFile->Close(); + + delete rl; + return 0; +} + + diff --git a/ITS/UPGRADE/macros/QA/AliTPCUComparison.C b/ITS/UPGRADE/macros/QA/AliTPCUComparison.C new file mode 100644 index 00000000000..077c6d1975e --- /dev/null +++ b/ITS/UPGRADE/macros/QA/AliTPCUComparison.C @@ -0,0 +1,626 @@ +/**************************************************************************** + * Legacy comparison macro, adapted to the upgraded ITS. * + * * + * Creates list of "trackable" tracks, * + * calculates efficiency, resolutions etc. * + * There is a possibility to run this macro over several events. * + * * + * Origin: I.Belikov, IPHC, Iouri.Belikov@iphc.cnrs.fr * + * 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 + + #include "AliStack.h" + #include "AliHeader.h" + #include "AliTrackReference.h" + #include "AliRunLoader.h" + #include "AliRun.h" + #include "AliMagF.h" + #include "AliESDEvent.h" + #include "AliESDtrack.h" + + #include "Base/AliSimDigits.h" + #include "Base/AliTPCParamSR.h" + #include "Base/AliTPCLoader.h" + #include "Base/AliTPCcalibDB.h" + #include "Sim/AliTPC.h" + #include "Rec/AliTPCClustersRow.h" + + #include "AliCDBManager.h" +#endif + +Int_t GoodTracksTPC(const Char_t *dir="."); + +extern TBenchmark *gBenchmark; +extern TROOT *gROOT; + +static Int_t allgood=0; +static Int_t allselected=0; +static Int_t allfound=0; + +Int_t AliTPCUComparison +(Float_t ptcutl=0., Float_t ptcuth=2., const Char_t *dir=".") { + gBenchmark->Start("AliTPCUComparison"); + + ::Info("AliTPCUComparison.C","Doing comparison..."); + + + + 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); + + + Int_t nb=100; + TH1F *hgood=(TH1F*)gROOT->FindObject("hgood"); + if (!hgood) hgood=new TH1F("hgood","Good tracks",nb,ptcutl,ptcuth); + + TH1F *hfound=(TH1F*)gROOT->FindObject("hfound"); + if (!hfound) hfound=new TH1F("hfound","Found tracks",nb,ptcutl,ptcuth); + + TH1F *hfake=(TH1F*)gROOT->FindObject("hfake"); + if (!hfake) hfake=new TH1F("hfake","Fake tracks",nb,ptcutl,ptcuth); + + TH1F *hg=(TH1F*)gROOT->FindObject("hg"); + if (!hg) hg=new TH1F("hg","Efficiency for good tracks",nb,ptcutl,ptcuth); + hg->SetLineColor(4); hg->SetLineWidth(2); + + TH1F *hf=(TH1F*)gROOT->FindObject("hf"); + if (!hf) hf=new TH1F("hf","Efficiency for fake tracks",nb,ptcutl,ptcuth); + hf->SetFillColor(1); hf->SetFillStyle(3013); hf->SetLineWidth(2); + + + TH1F *he=(TH1F*)gROOT->FindObject("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); + + + Char_t fname[100]; + sprintf(fname,"%s/GoodTracksTPC.root",dir); + + TFile *refFile=TFile::Open(fname,"old"); + if (!refFile || !refFile->IsOpen()) { + ::Info("AliTPCUComparison.C","Marking good tracks (will take a while)..."); + if (GoodTracksTPC(dir)) { + ::Error("AliTPCUComparison.C","Can't generate the reference file !"); + return 1; + } + } + refFile=TFile::Open(fname,"old"); + if (!refFile || !refFile->IsOpen()) { + ::Error("AliTPCUComparison.C","Can't open the reference file !"); + return 2; + } + + TTree *tpcTree=(TTree*)refFile->Get("tpcTree"); + if (!tpcTree) { + ::Error("AliTPCUComparison.C","Can't get the reference tree !"); + return 3; + } + TBranch *branch=tpcTree->GetBranch("TPC"); + if (!branch) { + ::Error("AliTPCUComparison.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("AliTPCUComparison.C","Can't open AliESDtpc.root !"); + return 5; + } + } + AliESDEvent* event = new AliESDEvent(); + TTree* esdTree = (TTree*) ef->Get("esdTree"); + if (!esdTree) { + ::Error("AliTPCUComparison.C", "no ESD tree found"); + return 6; + } + event->ReadFromTree(esdTree); + + + //******* Loop over events ********* + Int_t e=0; + while (esdTree->GetEvent(e)) { + cout<GetNumberOfTracks(); + allfound+=nentr; + + if (tpcTree->GetEvent(e++)==0) { + cerr<<"No reconstructable tracks !\n"; + continue; + } + + Int_t ngood=refs->GetEntriesFast(); + allgood+=ngood; + + 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; + + Int_t i; + for (Int_t k=0; kUncheckedAt(k); + Int_t lab=ref->Label(), tlab=-1; + Float_t ptg=TMath::Sqrt(ref->Px()*ref->Px() + ref->Py()*ref->Py()); + + if (ptg<1e-33) continue; // for those not crossing 0 pad row + + if (ptgptcuth) continue; + + allselected++; + + hgood->Fill(ptg); + + 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.); + } + } + + cout<<"\nList of Not found tracks :\n"; + for ( i = 0; i< itrack_notfound; i++){ + cout<Draw(); + p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10); + 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)"); + 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("(%)"); + 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("(%)"); + 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); + hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2(); + hg->Divide(hfound,hgood,1,1.,"b"); + hf->Divide(hfake,hgood,1,1.,"b"); + hg->SetMaximum(1.4); + hg->SetYTitle("Tracking efficiency"); + hg->SetXTitle("Pt (GeV/c)"); + hg->Draw(); + + TLine *line1 = new TLine(ptcutl,1.0,ptcuth,1.0); line1->SetLineStyle(4); + line1->Draw("same"); + TLine *line2 = new TLine(ptcutl,0.9,ptcuth,0.9); line2->SetLineStyle(4); + line2->Draw("same"); + + hf->SetFillColor(1); + hf->SetFillStyle(3013); + hf->SetLineColor(2); + hf->SetLineWidth(2); + hf->Draw("histsame"); + TText *text = new TText(0.461176,0.248448,"Fake tracks"); + text->SetTextSize(0.05); + text->Draw(); + text = new TText(0.453919,1.11408,"Good tracks"); + text->SetTextSize(0.05); + text->Draw(); + + + + 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"); + 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("AliTPCUComparison.root","RECREATE"); + c1->Write(); + c2->Write(); + fc.Close(); + + gBenchmark->Stop("AliTPCUComparison"); + gBenchmark->Show("AliTPCUComparison"); + + return 0; +} + + +Int_t GoodTracksTPC(const Char_t *dir) { + Char_t fname[100]; + sprintf(fname,"%s/galice.root",dir); + + AliRunLoader *rl = AliRunLoader::Open(fname,"COMPARISON"); + if (!rl) { + ::Error("GoodTracksTPC","Can't start session !"); + return 1; + } + + rl->LoadgAlice(); + rl->LoadHeader(); + rl->LoadKinematics(); + rl->LoadTrackRefs(); + + 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(); + cout<<"TPC version "<LoadRecPoints(); + else if (ver==2) tpcl->LoadDigits(); + else { + ::Error("GoodTracksTPC","Wrong TPC version !"); + delete rl; + return 3; + } + + TGeoGlobalMagField::Instance()->SetField( + new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG)); + + AliCDBManager *man=AliCDBManager::Instance(); + man->SetDefaultStorage("local://$ALICE_ROOT/OCDB"); + man->SetRun(0); + AliTPCParamSR *digp= + (AliTPCParamSR*)(AliTPCcalibDB::Instance()->GetParameters()); + if (!digp) { + ::Error("AliTPCUComparison.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), 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(); + + Int_t i; + + rl->GetEvent(e); file->cd(); + + Int_t np = rl->GetHeader()->GetNtrack(); + cout<<"Event "<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: + { + 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]|=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; + } + + + //****** select tracks which are "good" enough + AliStack* stack = rl->Stack(); + for (i=0; iParticle(i); + if (p == 0x0) { + cerr<<"Can not get particle "<Pt() <= 0.) 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++; + } + + //**** 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; +} diff --git a/ITS/UPGRADE/macros/QA/check_radii.C b/ITS/UPGRADE/macros/QA/check_radii.C new file mode 100644 index 00000000000..bcf3de995ad --- /dev/null +++ b/ITS/UPGRADE/macros/QA/check_radii.C @@ -0,0 +1,141 @@ + +#if !defined(__CINT__) || defined(__MAKECINT__) + #include + #include + #include + #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 "UPGRADE/AliITSUClusterPix.h" + #include "UPGRADE/AliITSULoader.h" + #include "UPGRADE/AliITSUGeomTGeo.h" +#endif + + +Int_t check_radii(const Char_t *dir=".") { + TFile *f=TFile::Open("xyz.root","recreate"); + TNtuple *nt=new TNtuple("nt","my ntuple","x:y:z"); + + if (gAlice) { + delete AliRunLoader::Instance(); + delete gAlice;//if everything was OK here it is already NULL + gAlice = 0x0; + } + + Char_t fname[100]; + sprintf(fname,"%s/galice.root",dir); + + AliRunLoader *rl = AliRunLoader::Open(fname,"COMPARISON"); + if (!rl) { + ::Error("GoodTracksITS","Can't start session !"); + return 1; + } + + rl->LoadgAlice(); + rl->LoadHeader(); + rl->LoadKinematics(); + + AliITSULoader* itsl = (AliITSULoader*)rl->GetLoader("ITSLoader"); + if (itsl == 0x0) { + ::Error("GoodTracksITS","Can not find the ITSLoader"); + delete rl; + return 4; + } + itsl->LoadRecPoints(); + + TGeoManager::Import("geometry.root"); + AliITSUGeomTGeo* gm = new AliITSUGeomTGeo(kTRUE,kTRUE); + AliITSUClusterPix::SetGeom(gm); + + Int_t nev=rl->GetNumberOfEvents(); + ::Info("GoodTracksITS","Number of events : %d\n",nev); + + //******** Loop over generated events + for (Int_t e=0; eGetEvent(e); + + TTree *cTree=itsl->TreeR(); + if (!cTree) { + ::Error("GoodTracksITS","Can't get the cluster tree !"); + delete rl; + return 8; + } + + const Int_t nLayers=7; + TBranch *branch[nLayers]; + TClonesArray clusters[nLayers]; + for (Int_t layer=0; layerGetBranch(bname); + if (!branch[layer]) { + ::Error("GoodTracksITS","Can't get the clusters branch !"); + delete rl; + return 9; + } + branch[layer]->SetAddress(&ptr); + } + + Int_t entr=(Int_t)cTree->GetEntries(); + for (k=0; kGetEvent(k); + for (Int_t lay=0; layGetGlobalXYZ(g); + //pnt->GetLocalXYZ(g); + //cout<Fill(g); + } + clusters[lay].Clear(); + } + } + + + } //*** end of the loop over generated events + + nt->Draw("y:x"); + Float_t rmin=19.44; + TEllipse *ellipse = new TEllipse(0,0,rmin,rmin,0,360,0); + ellipse->SetFillStyle(0); + ellipse->SetLineColor(2); + ellipse->Draw(); + + Float_t rmax=19.77; + ellipse = new TEllipse(0,0,rmax,rmax,0,360,0); + ellipse->SetFillStyle(0); + ellipse->SetLineColor(4); + ellipse->Draw(); + + delete rl; + return 0; +} + + diff --git a/ITS/UPGRADE/macros/QA/compClusHits.C b/ITS/UPGRADE/macros/QA/compClusHits.C new file mode 100644 index 00000000000..dfba17898a0 --- /dev/null +++ b/ITS/UPGRADE/macros/QA/compClusHits.C @@ -0,0 +1,541 @@ +#if !defined(__CINT__) || defined(__MAKECINT__) +#include "TObjArray.h" +#include "TFile.h" +#include "TTree.h" +#include "TH1F.h" +#include "../ITS/UPGRADE/AliITSUClusterPix.h" +#include "../ITS/UPGRADE/AliITSURecoLayer.h" +#include "../ITS/UPGRADE/AliITSURecoDet.h" +#include "../ITS/UPGRADE/AliITSUHit.h" +#include "../ITS/UPGRADE/AliITSUGeomTGeo.h" +#include "AliITSsegmentation.h" +#include "AliGeomManager.h" +#include "AliStack.h" +#include "AliLoader.h" +#include "AliCDBManager.h" + +#include "TROOT.h" +#include "TStyle.h" +#include "TGeoMatrix.h" +#include "TParticle.h" +#include "TCanvas.h" +#include "TPaveStats.h" +#include "TClonesArray.h" + +#endif + +TObjArray histoArr; +enum {kNPixAll=0,kNPixSPL=1,kDR=0,kDTXodd,kDTXeven,kDTZ, kDTXoddSPL,kDTXevenSPL,kDTZSPL}; + +TPaveStats* GetStPad(TH1* hst); +TPaveStats* SetStPadPos(TH1* hst,float x1,float x2,float y1,float y2, Int_t stl=-1,Int_t col=-1); +TCanvas* DrawNP(int np, TObjArray* harr=0, TCanvas* cnv=0); +TH1* GetHistoClSize(int npix,int id,TObjArray* harr=0); +void DrawReport(const char* psname, TObjArray* harr=0); + + +typedef struct { + Int_t evID; + Int_t volID; + Int_t lrID; + Int_t clID; + Int_t nPix; + Int_t nX; + Int_t nZ; + Int_t q; + Float_t pt; + Float_t eta; + Float_t phi; + Float_t xyz[3]; + Float_t dX; + Float_t dY; + Float_t dZ; + Bool_t split; + Bool_t prim; + Int_t pdg; + Int_t ntr; +} clSumm; + +void compClusHits(int nev=-1) +{ + const int kSplit=0x1<<22; + const int kSplCheck=0x1<<23; + // + gSystem->Load("libITSUpgradeBase"); + gSystem->Load("libITSUpgradeSim"); + gSystem->Load("libITSUpgradeRec"); + gROOT->SetStyle("Plain"); + + AliCDBManager* man = AliCDBManager::Instance(); + man->SetDefaultStorage("local://$ALICE_ROOT/OCDB"); + man->SetSpecificStorage("GRP/GRP/Data", + Form("local://%s",gSystem->pwd())); + man->SetSpecificStorage("ITS/Align/Data", + Form("local://%s",gSystem->pwd())); + man->SetSpecificStorage("ITS/Calib/RecoParam", + Form("local://%s",gSystem->pwd())); + man->SetRun(0); + + + gAlice=NULL; + AliRunLoader* runLoader = AliRunLoader::Open("galice.root"); + runLoader->LoadgAlice(); + + gAlice = runLoader->GetAliRun(); + + runLoader->LoadHeader(); + runLoader->LoadKinematics(); + runLoader->LoadRecPoints(); + runLoader->LoadSDigits(); + runLoader->LoadHits(); + + AliLoader *dl = runLoader->GetDetectorLoader("ITS"); + + AliGeomManager::LoadGeometry("geometry.root"); + TObjArray algITS; + AliGeomManager::LoadAlignObjsFromCDBSingleDet("ITS",algITS); + AliGeomManager::ApplyAlignObjsToGeom(algITS); + // + AliITSUGeomTGeo* gm = new AliITSUGeomTGeo(kTRUE); + AliITSUClusterPix::SetGeom(gm); + // + AliITSURecoDet *its = new AliITSURecoDet(gm, "ITSinterface"); + its->CreateClusterArrays(); + // + Double_t xg1,yg1,zg1=0.,xg0,yg0,zg0=0.,tg0; + // + TTree * cluTree = 0x0; + TTree *hitTree = 0x0; + TClonesArray *hitList=new TClonesArray("AliITSUHit"); + // + TObjArray arrMCTracks; // array of hit arrays for each particle + // + Float_t xyzClGloF[3]; + Double_t xyzClGlo[3],xyzClTr[3]; + Int_t labels[3]; + int nLab = 0; + int nlr=its->GetNLayersActive(); + int ntotev = (Int_t)runLoader->GetNumberOfEvents(); + printf("N Events : %i \n",ntotev); + if (nev>0) ntotev = TMath::Min(nev,ntotev); + // + // output tree + TFile* flOut = TFile::Open("clInfo.root","recreate"); + TTree* trOut = new TTree("clitsu","clitsu"); + clSumm cSum; + trOut->Branch("evID", &cSum.evID ,"evID/I"); + trOut->Branch("volID",&cSum.volID,"volID/I"); + trOut->Branch("lrID", &cSum.lrID ,"lrID/I"); + trOut->Branch("clID", &cSum.clID ,"clID/I"); + trOut->Branch("nPix", &cSum.nPix ,"nPix/I"); + trOut->Branch("nX" , &cSum.nX ,"nX/I"); + trOut->Branch("nZ" , &cSum.nZ ,"nZ/I"); + trOut->Branch("q" , &cSum.q ,"q/I"); + trOut->Branch("pt" , &cSum.pt ,"pt/F"); + trOut->Branch("eta" ,&cSum.eta ,"eta/F"); + trOut->Branch("phi" , &cSum.phi ,"phi/F"); + trOut->Branch("xyz", cSum.xyz, "xyz[3]/F"); + trOut->Branch("dX" , &cSum.dX ,"dX/F"); + trOut->Branch("dY" , &cSum.dY ,"dY/F"); + trOut->Branch("dZ" , &cSum.dZ ,"dZ/F"); + trOut->Branch("split",&cSum.split,"split/O"); + trOut->Branch("prim", &cSum.prim, "prim/O"); + trOut->Branch("pdg", &cSum.pdg, "pdg/I"); + trOut->Branch("ntr", &cSum.ntr, "ntr/I"); + // + for (Int_t iEvent = 0; iEvent < ntotev; iEvent++) { + printf("\n Event %i \n",iEvent); + runLoader->GetEvent(iEvent); + AliStack *stack = runLoader->Stack(); + cluTree=dl->TreeR(); + hitTree=dl->TreeH(); + hitTree->SetBranchAddress("ITS",&hitList); + // + // read clusters + for (int ilr=nlr;ilr--;) { + TBranch* br = cluTree->GetBranch(Form("ITSRecPoints%d",ilr)); + if (!br) {printf("Did not find cluster branch for lr %d\n",ilr); exit(1);} + br->SetAddress(its->GetLayerActive(ilr)->GetClustersAddress()); + } + cluTree->GetEntry(0); + its->ProcessClusters(); + // + // read hits + for(Int_t iEnt=0;iEntGetEntries();iEnt++){//entries loop degli hits + hitTree->GetEntry(iEnt); + int nh = hitList->GetEntries(); + for(Int_t iHit=0; iHitAt(iHit); + int mcID = pHit->GetTrack(); + TClonesArray* harr = arrMCTracks.GetEntriesFast()>mcID ? (TClonesArray*)arrMCTracks.At(mcID) : 0; + if (!harr) { + harr = new TClonesArray("AliITSUHit"); // 1st encounter of the MC track + arrMCTracks.AddAtAndExpand(harr,mcID); + } + // + new ( (*harr)[harr->GetEntriesFast()] ) AliITSUHit(*pHit); + } + } + // + // compare clusters and hits + // + printf(" tree entries: %lld\n",cluTree->GetEntries()); + // + for (int ilr=0;ilrGetLayerActive(ilr); + TClonesArray* clr = lr->GetClusters(); + int nClu = clr->GetEntries(); + printf("Layer %d : %d clusters\n",ilr,nClu); + // + for (int icl=0;iclAt(icl); + int modID = cl->GetVolumeId(); + + //------------ check if this is a split cluster + int sInL = modID - gm->GetFirstChipIndex(ilr); + if (!cl->TestBit(kSplCheck)) { + cl->SetBit(kSplCheck); + // check if there is no other cluster with same label on this module + AliITSURecoSens* sens = lr->GetSensor(sInL); + int nclSn = sens->GetNClusters(); + int offs = sens->GetFirstClusterId(); + // printf("To check for %d (mod:%d) N=%d from %d\n",icl,modID,nclSn,offs); + for (int ics=0;icsGetCluster(offs+ics); // access to clusters + if (clusT==cl) continue; + for (int ilb0=0;ilb0<3;ilb0++) { + int lb0 = cl->GetLabel(ilb0); if (lb0<=-1) break; + for (int ilb1=0;ilb1<3;ilb1++) { + int lb1 = clusT->GetLabel(ilb1); if (lb1<=-1) break; + if (lb1==lb0) { + cl->SetBit(kSplit); + clusT->SetBit(kSplit); + /* + printf("Discard clusters of module %d:\n",modID); + cl->Print(); + clusT->Print(); + */ + break; + } + } + } + } + } + //------------ + const AliITSsegmentation* segm = gm->GetSegmentation(ilr); + // + cl->GetGlobalXYZ(xyzClGloF); + int clsize = cl->GetNPix(); + for (int i=3;i--;) xyzClGlo[i] = xyzClGloF[i]; + const TGeoHMatrix* mat = gm->GetMatrixSens(modID); + if (!mat) {printf("failed to get matrix for module %d\n",cl->GetVolumeId());} + mat->MasterToLocal(xyzClGlo,xyzClTr); + // + int col,row; + segm->LocalToDet(xyzClTr[0],xyzClTr[2],row,col); // effective col/row + nLab = 0; + for (int il=0;il<3;il++) { + if (cl->GetLabel(il)>=0) labels[nLab++] = cl->GetLabel(il); + else break; + } + // find hit info + for (int il=0;ilPrint(); continue;} + // + int nh = htArr->GetEntriesFast(); + AliITSUHit *pHit=0; + for (int ih=nh;ih--;) { + AliITSUHit* tHit = (AliITSUHit*)htArr->At(ih); + if (tHit->GetChip()!=modID) continue; + pHit = tHit; + break; + } + if (!pHit) { + printf("did not find MChit for label %d on module %d ",il,modID); + cl->Print(); + htArr->Print(); + continue; + } + // + pHit->GetPositionG(xg1,yg1,zg1); + pHit->GetPositionG0(xg0,yg0,zg0,tg0); + // + double txyzH[3],gxyzH[3] = { (xg1+xg0)/2, (yg1+yg0)/2, (zg1+zg0)/2 }; + mat->MasterToLocal(gxyzH,txyzH); + double rcl = TMath::Sqrt(xyzClTr[0]*xyzClTr[0]+xyzClTr[1]*xyzClTr[1]); + double rht = TMath::Sqrt(txyzH[0]*txyzH[0]+txyzH[1]*txyzH[1]); + // + GetHistoClSize(clsize,kDR,&histoArr)->Fill((rht-rcl)*1e4); + if (cl->TestBit(kSplit)) { + if (col%2) GetHistoClSize(clsize,kDTXoddSPL,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4); + else GetHistoClSize(clsize,kDTXevenSPL,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4); + GetHistoClSize(clsize,kDTZSPL,&histoArr)->Fill((txyzH[2]-xyzClTr[2])*1e4); + GetHistoClSize(0,kNPixSPL,&histoArr)->Fill(clsize); + } + if (col%2) GetHistoClSize(clsize,kDTXodd,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4); + else GetHistoClSize(clsize,kDTXeven,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4); + GetHistoClSize(clsize,kDTZ,&histoArr)->Fill((txyzH[2]-xyzClTr[2])*1e4); + GetHistoClSize(0,kNPixAll,&histoArr)->Fill(clsize); + // + cSum.evID = iEvent; + cSum.volID = cl->GetVolumeId(); + cSum.lrID = ilr; + cSum.clID = icl; + cSum.nPix = cl->GetNPix(); + cSum.nX = cl->GetNx(); + cSum.nZ = cl->GetNz(); + cSum.q = cl->GetQ(); + cSum.split = cl->TestBit(kSplit); + cSum.dX = (txyzH[0]-xyzClTr[0])*1e4; + cSum.dY = (txyzH[1]-xyzClTr[1])*1e4; + cSum.dZ = (txyzH[2]-xyzClTr[2])*1e4; + int label = cl->GetLabel(0); + TParticle* part = 0; + if (label>=0 && (part=stack->Particle(label)) ) { + cSum.pdg = part->GetPdgCode(); + cSum.eta = part->Eta(); + cSum.pt = part->Pt(); + cSum.phi = part->Phi(); + cSum.prim = stack->IsPhysicalPrimary(label); + } + cSum.ntr = 0; + for (int ilb=0;ilb<3;ilb++) if (cl->GetLabel(ilb)>=0) cSum.ntr++; + for (int i=0;i<3;i++) cSum.xyz[i] = xyzClGloF[i]; + // + trOut->Fill(); + /* + if (clsize==5) { + printf("\nL%d(%c) Mod%d, Cl:%d | %+5.1f %+5.1f (%d/%d)|H:%e %e %e | C:%e %e %e\n",ilr,cl->TestBit(kSplit) ? 'S':'N', + modID,icl,(txyzH[0]-xyzClTr[0])*1e4,(txyzH[2]-xyzClTr[2])*1e4, row,col, + gxyzH[0],gxyzH[1],gxyzH[2],xyzClGlo[0],xyzClGlo[1],xyzClGlo[2]); + cl->Print(); + pHit->Print(); + // + double a0,b0,c0,a1,b1,c1,e0; + pHit->GetPositionL0(a0,b0,c0,e0); + pHit->GetPositionL(a1,b1,c1); + float cloc[3]; + cl->GetLocalXYZ(cloc); + printf("LocH: %e %e %e | %e %e %e\n",a0,b0,c0,a1,b1,c1); + printf("LocC: %e %e %e | %e %e %e\n",cloc[0],cloc[1],cloc[2],xyzClTr[0],xyzClTr[1],xyzClTr[2]); + } + */ + // + } + } + } + + // layerClus.Clear(); + // + arrMCTracks.Delete(); + }//event loop + // + flOut->cd(); + trOut->Write(); + delete trOut; + flOut->Close(); + flOut->Delete(); + DrawReport("clinfo.ps",&histoArr); + // +} + +void DrawReport(const char* psname, TObjArray* harr) +{ + gStyle->SetOptFit(1); + if (!harr) harr = &histoArr; + TCanvas* cnv = new TCanvas("cl","cl",900,600); + // + TString psnm1 = psname; + if (psnm1.IsNull()) psnm1 = "clusters.ps"; + TString psnm0 = psnm1.Data(); + psnm0 += "["; + TString psnm2 = psnm1.Data(); + psnm2 += "]"; + cnv->Print(psnm0.Data()); + // + TH1* clall = GetHistoClSize(0,kNPixAll,harr); + clall->SetLineColor(kRed); + clall->Draw(); + TH1* clSpl = GetHistoClSize(0,kNPixSPL,harr); + clSpl->SetLineColor(kBlue); + clSpl->Draw("sames"); + gPad->Modified(); + gPad->Update(); + SetStPadPos(clall,0.75,0.97,0.8,1.,-1,clall->GetLineColor()); + SetStPadPos(clSpl,0.75,0.97,0.6,0.8,-1,clSpl->GetLineColor()); + gPad->Modified(); + gPad->Update(); + gPad->SetLogy(1); + // + cnv->cd(); + cnv->Print(psnm1.Data()); + // + // plot cluster sized from 1 to 10 + for (int i=1;i<=10;i++) { + if (clall->GetBinContent(clall->FindBin(i))<100) continue; + DrawNP(i,harr,cnv); + cnv->Print(psnm1.Data()); + } + cnv->Print(psnm2.Data()); +} + +//------------------------------------------ +TH1* GetHistoClSize(int npix,int id,TObjArray* harr) +{ + // book histos + TH1* h = 0; + if (!harr) harr = &histoArr; + // + if (npix<1) { + if (harr->GetEntriesFast()>=id && (h=(TH1*)harr->At(id))) return h; + h = new TH1F("npixAll","npixAll",150,0.5,54.5); + h->SetDirectory(0); + h->SetLineColor(kRed); + harr->AddAtAndExpand(h, kNPixAll); + // + h = new TH1F("npixSpl","npixSpl",150,0.5,54.5); + h->SetLineColor(kBlue); + h->SetDirectory(0); + harr->AddAtAndExpand(h, kNPixSPL); + // + h = (TH1*)harr->At(id); + if (!h) {printf("Unknown histo id=%d\n",id); exit(1);} + return h; + } + // + int idh = npix*10+id; + if (harr->GetEntriesFast()>=idh && (h=(TH1*)harr->At(idh))) return h; + // + const int nbin=100; + const double kdiff=80; + // need to create set of histos + // + h = new TH1F(Form("dxy_npix%d",npix),Form("dr_npix%d",npix),nbin,-kdiff,kdiff); + h->SetDirectory(0); + harr->AddAtAndExpand(h, npix*10 + kDR); + // + h = new TH1F(Form("dtxODD_npix%d",npix),Form("dtxODD_npix%d",npix),nbin,-kdiff,kdiff); + h->SetDirectory(0); + h->SetLineColor(kRed); + harr->AddAtAndExpand(h, npix*10 + kDTXodd); + h = new TH1F(Form("dtxEVN_npix%d",npix),Form("dtxEVN_npix%d",npix),nbin,-kdiff,kdiff); + h->SetDirectory(0); + h->SetLineColor(kBlue); + harr->AddAtAndExpand(h, npix*10 + kDTXeven); + // + h = new TH1F(Form("dtz_npix%d",npix),Form("dtz_npix%d",npix),nbin,-kdiff,kdiff); + h->SetLineColor(kGreen); + h->SetDirectory(0); + harr->AddAtAndExpand(h, npix*10 + kDTZ); + // + // + h = new TH1F(Form("SPL_dtxODD_npix%d",npix),Form("SPL_dtxODD_npix%d",npix),nbin,-kdiff,kdiff); + h->SetLineColor(kMagenta); + h->SetFillColor(kMagenta); + h->SetFillStyle(3001); + h->SetLineStyle(2); + h->SetDirectory(0); + + harr->AddAtAndExpand(h, npix*10 + kDTXoddSPL); + h = new TH1F(Form("SPL_dtxEVN_npix%d",npix),Form("SPL_dtxEVN_npix%d",npix),nbin,-kdiff,kdiff); + h->SetLineColor(kCyan); + h->SetFillColor(kCyan); + h->SetFillStyle(3006); + h->SetLineStyle(2); + h->SetDirectory(0); + harr->AddAtAndExpand(h, npix*10 + kDTXevenSPL); + // + h = new TH1F(Form("SPL_dtz_npix%d",npix),Form("SPLdtz_npix%d",npix),nbin,-kdiff,kdiff); + harr->AddAtAndExpand(h, npix*10 + kDTZSPL); + h->SetDirectory(0); + // + h->SetLineColor(kGreen+2); + h->SetFillColor(kGreen+2); + h->SetLineStyle(2); + h->SetFillStyle(3001); + h = (TH1*)harr->At(idh); + if (!h) {printf("Unknown histo id=%d\n",idh); exit(1);} + return h; +} + + + +TCanvas* DrawNP(int np, TObjArray* harr, TCanvas* cnv) +{ + if (!harr) harr = &histoArr; + if (!cnv) cnv = new TCanvas(Form("cnv%d",np),Form("cnv%d",np),900,700); + cnv->Clear(); + cnv->Divide(2,1); + cnv->cd(1); + // + TH1* dxodd = (TH1*)harr->At(np*10+kDTXodd); + TH1* dxevn = (TH1*)harr->At(np*10+kDTXeven); + TH1* dxoddS =(TH1*)harr->At(np*10+kDTXoddSPL); + TH1* dxevnS =(TH1*)harr->At(np*10+kDTXevenSPL); + double max = TMath::Max(dxodd->GetMaximum(),dxevn->GetMaximum()); + dxodd->SetMaximum(1.1*max); + dxodd->GetXaxis()->SetTitle("#DeltaX, #mum"); + dxodd->SetTitle(Form("#DeltaX for clSize=%d",np)); + dxodd->Fit("gaus","",""); + dxevn->Fit("gaus","","sames"); + // + dxoddS->Draw("sames"); + dxevnS->Draw("sames"); + // + gPad->Modified(); + gPad->Update(); + SetStPadPos(dxodd,0.75,0.97,0.8,1., -1,dxodd->GetLineColor()); + SetStPadPos(dxevn,0.75,0.97,0.6,0.8, -1,dxevn->GetLineColor()); + SetStPadPos(dxoddS,0.75,0.97,0.4,0.6, -1,dxoddS->GetLineColor()); + SetStPadPos(dxevnS,0.75,0.97,0.2,0.4, -1,dxevnS->GetLineColor()); + // + cnv->cd(2); + TH1* dz = (TH1*)harr->At(np*10+kDTZ); + dz->SetTitle(Form("#DeltaZ for clSize=%d",np)); + dz->GetXaxis()->SetTitle("#DeltaZ, #mum"); + dz->Fit("gaus"); + TH1* dzS = (TH1*)harr->At(np*10+kDTZSPL); + dz->Draw("sames"); + gPad->Modified(); + gPad->Update(); + SetStPadPos(dz,0.75,0.97,0.8,1., -1, dz->GetLineColor()); + SetStPadPos(dzS,0.75,0.97,0.5,0.7, -1, dzS->GetLineColor()); + gPad->Modified(); + gPad->Update(); + // + cnv->cd(); + return cnv; +} + + +TPaveStats* GetStPad(TH1* hst) +{ + TList *lst = hst->GetListOfFunctions(); + if (!lst) return 0; + int nf = lst->GetSize(); + for (int i=0;iAt(i); + if (fnc->InheritsFrom("TPaveStats")) return fnc; + } + return 0; + // +} + + +TPaveStats* SetStPadPos(TH1* hst,float x1,float x2,float y1,float y2, Int_t stl, Int_t col) +{ + TPaveStats* pad = GetStPad(hst); + if (!pad) return 0; + pad->SetX1NDC( x1 ); + pad->SetX2NDC( x2 ); + pad->SetY1NDC( y1 ); + pad->SetY2NDC( y2 ); + if (stl>=0) pad->SetFillStyle(stl); + if (col>=0) pad->SetTextColor(col); + pad->SetFillColor(0); + // + gPad->Modified(); + return pad; +} + diff --git a/ITS/UPGRADE/macros/QA/plotHits.C b/ITS/UPGRADE/macros/QA/plotHits.C new file mode 100644 index 00000000000..7c25df3243c --- /dev/null +++ b/ITS/UPGRADE/macros/QA/plotHits.C @@ -0,0 +1,188 @@ +void plotHits(Int_t ilay=3, Int_t ilrs=1, Int_t nev=-1, Int_t evStart=0) +{ + // Adapted from readHit.C macro - M.S. 10 jul 14 + + enum {kLSLayer, kLSStave , kLSHalfStave, kLSModule, + kLSChip , kLSSensor, kLSNumPar }; + char lrsNames[kLSNumPar][10]; + Int_t i = -1; + strcpy(lrsNames[++i],"Layer"); + strcpy(lrsNames[++i],"Stave"); + strcpy(lrsNames[++i],"HalfStave"); + strcpy(lrsNames[++i],"Module"); + strcpy(lrsNames[++i],"Chip"); + strcpy(lrsNames[++i],"Sensor"); + + Int_t idWrapVol[7] = {0, 0, 0, 1, 1, 2, 2}; + + Bool_t argerr = kFALSE; + + if (ilay < 0 || ilay > 6 || + ilrs < 0 || ilrs >= kLSNumPar ) argerr = kTRUE; + + if (argerr) { + printf("Wrong parameters! ilay = %d - ilrs = %d\n\n",ilay,ilrs); + printf("Usage:\n.x plotHit(ilay, ilrs, nev, evStart)\n"); + printf(" ilay : layer number (0 < ilay < 7, default 0)\n"); + printf(" ilrs : local reference system (0 layer, 1 stave, 2 half stave, 3 module, 4 chip, 5 sensor, default 1)\n"); + printf(" nev : number of events (default -1 = all)\n"); + printf(" evStart : first event (default 0)\n"); + return; + } + + printf("--> Plotting hits at %s level for Layer %d <--\n", + lrsNames[ilrs],ilay); + + gROOT->SetStyle("Plain"); + + // Load libraries, geometry, pointers, whatever + gSystem->Load("libITSUpgradeBase"); + gSystem->Load("libITSUpgradeSim"); + + gAlice=NULL; + AliRunLoader* runLoader = AliRunLoader::Open("galice.root"); + runLoader->LoadgAlice(); + + gAlice = runLoader->GetAliRun(); + + runLoader->LoadHeader(); + runLoader->LoadKinematics(); + runLoader->LoadSDigits(); + runLoader->LoadHits(); + + AliGeomManager::LoadGeometry("geometry.root"); + AliITSUGeomTGeo* gm = new AliITSUGeomTGeo(kTRUE); + + AliLoader *dl = runLoader->GetDetectorLoader("ITS"); + + // Create histos + char volName[15]; + sprintf(volName,"ITSU%s%d",lrsNames[ilrs],ilay); + TGeoBBox *vShape = (TGeoBBox*)(gGeoManager->GetVolume(volName)->GetShape()); + Double_t xHist = 1.1*vShape->GetDX(); // 10% greater + Double_t yHist = 1.1*vShape->GetDY(); + Double_t zHist = 1.1*vShape->GetDZ(); + + // OB stave are strongly asymmetric in Y: dimension increased for better view + if (ilay > 2 && ilrs == kLSStave) yHist*=1.20; + + char histTitle[50]; + + sprintf(histTitle," X - Y %s Local coordinates",lrsNames[ilrs]); + TH2F *xy = new TH2F("xy",histTitle,100,-xHist,xHist,100,-yHist,yHist); + xy->SetXTitle("X (cm)"); + xy->SetYTitle("Y (cm)"); + xy->SetMarkerStyle(7); + + sprintf(histTitle," X - Z %s Local coordinates",lrsNames[ilrs]); + TH2F *xz = new TH2F("xz",histTitle,100,-xHist,xHist,100,-zHist,zHist); + xz->SetXTitle("X (cm)"); + xz->SetYTitle("Z (cm)"); + xz->SetMarkerStyle(7); + + sprintf(histTitle," Z - Y %s Local coordinates",lrsNames[ilrs]); + TH2F *zy = new TH2F("zy",histTitle,100,-zHist,zHist,100,-yHist,yHist); + zy->SetXTitle("Z (cm)"); + zy->SetYTitle("Y (cm)"); + zy->SetMarkerStyle(7); + + // Get number of events which to loop on + Int_t nevTot = (Int_t)runLoader->GetNumberOfEvents(); + printf("Total events : %i\n\n",nevTot); + evStart = evStart nevTot) lastEv = nevTot; + + // Loop on hit structure + TTree *hitTree = 0x0; + TClonesArray *hitList=new TClonesArray("AliITSUHit"); + + TString path; + TGeoMatrix *mat = 0x0; + char nodName[15], parName[15]; + + for (Int_t iEvent = evStart; iEvent < lastEv; iEvent++) { + + printf("Event\t%d\n",iEvent); + + runLoader->GetEvent(iEvent); + + hitTree = dl->TreeH(); + hitTree->SetBranchAddress("ITS",&hitList); + + for(Int_t iEnt=0; iEntGetEntries(); iEnt++){ + hitTree->GetEntry(iEnt); + + for(Int_t iHit=0; iHitGetEntries();iHit++){ + AliITSUHit *pHit = (AliITSUHit*)hitList->At(iHit); + + Int_t id = pHit->GetChip(); + Int_t lr = gm->GetLayer(id); + Int_t st = gm->GetStave(id); + Int_t hs = gm->GetHalfStave(id); + Int_t md = gm->GetModule(id); + Int_t cp = gm->GetChipIdInModule(id); + + if (lr == ilay) { + Double_t xg, yg, zg=0.; + Double_t loc[3], mas[3]; + + pHit->GetPositionG(xg, yg, zg); + mas[0] = xg; + mas[1] = yg; + mas[2] = zg; + path.Clear(); + path = Form("/ALIC_1/ITSV_2/"); + path += Form("ITSUWrapVol%d_1/ITSULayer%d_1",idWrapVol[ilay],ilay); + if (ilrs > kLSLayer) + path += Form("/ITSUStave%d_%d",ilay,st); + if (ilrs > kLSStave) + path += Form("/ITSUHalfStave%d_%d",ilay,hs); + if (ilrs > kLSHalfStave) + path += Form("/ITSUModule%d_%d",ilay,md); + if (ilrs > kLSModule) + path += Form("/ITSUChip%d_%d",ilay,cp); + if (ilrs > kLSChip) + path += Form("/ITSUSensor%d_1",ilay); + + // Get the matrix for the current volume + gGeoManager->PushPath(); + if (!gGeoManager->cd(path.Data())) { + gGeoManager->PopPath(); + printf("Error in cd-ing to %s!\n",path.Data()); + printf("Chip id %5d | Lr:%2d Stave:%3d, HStave:%3d Module:%3d Chip:%3d\n", + id,lr,st,hs,md,cp); + + return 1; + } + mat = gGeoManager->GetCurrentMatrix(); + gGeoManager->PopPath(); + + // Convert global coordinate to local reference system + mat->MasterToLocal(mas, loc); + xy->Fill(loc[0],loc[1]); + xz->Fill(loc[0],loc[2]); + zy->Fill(loc[2],loc[1]); + + } // lr == ilay + }//loop hit + }//entryloopHitList + + }//event loop + + // Finally plot the histos + TCanvas *xyCanv = new TCanvas("xyCanv","Hit X-Y positions",500,500); + xyCanv->cd(); + xy->Draw(); + + TCanvas *xzCanv = new TCanvas("xzCanv","Hit X-Z positions",500,500); + xzCanv->cd(); + xz->Draw(); + + TCanvas *zyCanv = new TCanvas("zyCanv","Hit Z-Y positions",500,500); + zyCanv->cd(); + zy->Draw(); + +} -- 2.43.0