--- /dev/null
+/****************************************************************************
+ * 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 <TMath.h>
+ #include <TError.h>
+ #include <Riostream.h>
+ #include <TH1F.h>
+ #include <TH2F.h>
+ #include <TTree.h>
+ #include <TParticle.h>
+ #include <TCanvas.h>
+ #include <TLine.h>
+ #include <TText.h>
+ #include <TBenchmark.h>
+ #include <TStyle.h>
+ #include <TFile.h>
+ #include <TROOT.h>
+
+ #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.4<p<0.5 GeV/c",50,0.,200.);
+
+ TH2F *hep=(TH2F*)gROOT->FindObject("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<<endl<<endl<<"********* Processing event number: "<<e<<"*******\n";
+
+ Int_t nentr=event->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; k<ngood; k++) {
+ AliTrackReference *ref=(AliTrackReference*)refs->UncheckedAt(k);
+ Int_t lab=ref->Label(), tlab=-1;
+ Float_t ptg=TMath::Sqrt(ref->Px()*ref->Px() + ref->Py()*ref->Py());
+
+ if (ptg<ptcutl) continue;
+ if (ptg>ptcuth) continue;
+
+ allselected++;
+
+ hgood->Fill(ptg);
+
+ AliESDtrack *esd=0;
+ Int_t cnt=0;
+ for (Int_t i=0; i<nentr; i++) {
+ AliESDtrack *t=event->GetTrack(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; k<nnotf; k++){
+ cout<<notf[k]<<"\t";
+ if ((k%9)==8) cout<<"\n";
+ }
+ cout<<"\n\nList of fake tracks :\n";
+ for (k=0; k<nfake; k++){
+ cout<<fake[k]<<"\t";
+ if ((k%9)==8) cout<<"\n";
+ }
+ cout<<"\n\nList of multiple found tracks :\n";
+ for (k=0; k<nmult; k++) {
+ cout<<"id. "<<mult[k]
+ <<" found - "<<numb[k]<<"times\n";
+ }
+ cout<<endl;
+
+ cout<<"Number of found tracks : "<<nentr<<endl;
+ cout<<"Number of \"good\" tracks : "<<ngood<<endl;
+
+ refs->Clear();
+ } //***** 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 "<<nf/ng*100.<<" %\n";
+ cout<<"Total number selected of \"good\" tracks ="<<allselected<<endl<<endl;
+ cout<<"Total number of found tracks ="<<allfound<<endl;
+ cout<<"Total number of \"good\" tracks ="<<allgood<<endl;
+ cout<<endl;
+
+ gStyle->SetOptStat(111110);
+ gStyle->SetOptFit(1);
+
+ TCanvas *c1=new TCanvas("c1","",0,0,700,850);
+
+ Int_t minc=33;
+
+ TPad *p1=new TPad("p1","",0,0.3,.5,.6); p1->Draw();
+ p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10);
+ hp->SetFillColor(4); hp->SetXTitle("(mrad)");
+ if (hp->GetEntries()<minc) hp->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()<minc) hl->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()<minc) hpt->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()<minc) hmpt->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()<minc) he->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; e<nev; e++) {
+ Int_t k;
+
+ rl->GetEvent(e); itsFile->cd();
+
+ Int_t np = rl->GetHeader()->GetNtrack();
+ cout<<"Event "<<e<<" Number of particles: "<<np<<endl;
+
+ //******** Fill the "good" masks
+ Int_t *good=new Int_t[np]; for (k=0; k<np; k++) good[k]=0;
+
+ 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; layer<nLayers; layer++) {
+ TClonesArray *ptr =
+ new(clusters+layer) TClonesArray("AliITSUClusterPix",1000);
+ Char_t bname[33];
+ sprintf(bname,"ITSRecPoints%d\0",layer);
+ branch[layer]=cTree->GetBranch(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; k<entr; k++) {
+ cTree->GetEvent(k);
+ for (Int_t lay=0; lay<nLayers; lay++) {
+ Int_t ncl=clusters[lay].GetEntriesFast(); if (ncl==0) continue;
+ while (ncl--) {
+ AliITSUClusterPix *pnt=
+ (AliITSUClusterPix*)clusters[lay].UncheckedAt(ncl);
+ Int_t l0=pnt->GetLabel(0);
+ if (l0>=np) {
+// cerr<<"Wrong label: "<<l0<<endl;
+ continue;
+ }
+ Int_t l1=pnt->GetLabel(1);
+ if (l1>=np) {
+// cerr<<"Wrong label: "<<l1<<endl;
+ continue;
+ }
+ Int_t l2=pnt->GetLabel(2);
+ if (l2>=np) {
+// cerr<<"Wrong label: "<<l2<<endl;
+ continue;
+ }
+ Int_t mask=1<<lay;
+ if (l0>=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; k<nk; k++) {
+ AliTrackReference *tpcRef=(AliTrackReference *)tpcRefs->UncheckedAt(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 "<<lab<<endl;
+ continue;
+ }
+
+ AliTrackReference *ref=new((*itsRefs)[nt]) AliTrackReference(*tpcRef);
+ ref->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;
+}
+
+
--- /dev/null
+/****************************************************************************
+ * 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 <TMath.h>
+ #include <TError.h>
+ #include <Riostream.h>
+ #include <TH1F.h>
+ #include <TH2F.h>
+ #include <TTree.h>
+ #include <TParticle.h>
+ #include <TCanvas.h>
+ #include <TLine.h>
+ #include <TText.h>
+ #include <TBenchmark.h>
+ #include <TStyle.h>
+ #include <TFile.h>
+ #include <TROOT.h>
+
+ #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; i<nb; i++) {
+ Float_t c=h->GetBinContent(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; i<nb; i++) {
+ Float_t c=h->GetBinContent(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<<endl<<endl<<"********* Processing event number: "<<e<<"*******\n";
+
+ Int_t nentr=event->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; k<ngood; k++) {
+ AliTrackReference *ref=(AliTrackReference*)refs->UncheckedAt(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 (ptg<ptcutl) continue;
+ if (ptg>ptcuth) continue;
+
+ allselected++;
+
+ hgood->Fill(ptg);
+
+ AliESDtrack *esd=0;
+ Int_t cnt=0;
+ for (Int_t i=0; i<nentr; i++) {
+ AliESDtrack *t=event->GetTrack(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; k<nnotf; k++){
+ cout<<notf[k]<<"\t";
+ if ((k%9)==8) cout<<"\n";
+ }
+ cout<<"\n\nList of fake tracks :\n";
+ for (k=0; k<nfake; k++){
+ cout<<fake[k]<<"\t";
+ if ((k%9)==8) cout<<"\n";
+ }
+ cout<<"\n\nList of multiple found tracks :\n";
+ for (k=0; k<nmult; k++) {
+ cout<<"id. "<<mult[k]
+ <<" found - "<<numb[k]<<"times\n";
+ }
+ cout<<endl;
+
+ cout<<"Number of found tracks : "<<nentr<<endl;
+ cout<<"Number of \"good\" tracks : "<<ngood<<endl;
+
+ refs->Clear();
+ } //***** 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 "<<nf/ng*100.<<" %\n";
+ cout<<"Total number selected of \"good\" tracks ="<<allselected<<endl<<endl;
+ cout<<"Total number of found tracks ="<<allfound<<endl;
+ cout<<"Total number of \"good\" tracks ="<<allgood<<endl;
+ cout<<endl;
+
+ gStyle->SetOptStat(111110);
+ gStyle->SetOptFit(1);
+
+ TCanvas *c1=new TCanvas("c1","",0,0,700,850);
+
+ Int_t minc=33;
+
+ TPad *p1=new TPad("p1","",0,0.3,.5,.6); p1->Draw();
+ p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10);
+ hp->SetFillColor(4); hp->SetXTitle("(mrad)");
+ if (hp->GetEntries()<minc) hp->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()<minc) hl->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; e<nev; e++) {
+ Int_t k;
+
+ rl->GetEvent(e); itsFile->cd();
+
+ Int_t np = rl->GetHeader()->GetNtrack();
+ cout<<"Event "<<e<<" Number of particles: "<<np<<endl;
+
+ AliGenEventHeader *h=rl->GetHeader()->GenEventHeader();
+ TArrayF vtx(3);
+ h->PrimaryVertex(vtx);
+
+ Bool_t skip=kFALSE;
+ if (TMath::Abs(vtx[2]) > 10.) {
+ cout<<"Skipping an event with Zv="<<vtx[2]<<endl;
+ skip=kTRUE;
+ }
+
+ //******** Fill the "good" masks
+ Int_t *good=new Int_t[np]; for (k=0; k<np; k++) good[k]=0;
+
+ TTree *cTree=itsl->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; layer<nLayers; layer++) {
+ TClonesArray *ptr =
+ new(clusters+layer) TClonesArray("AliITSUClusterPix",1000);
+ Char_t bname[33];
+ sprintf(bname,"ITSRecPoints%d\0",layer);
+ branch[layer]=cTree->GetBranch(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; k<entr; k++) {
+ cTree->GetEvent(k);
+ for (Int_t lay=0; lay<nLayers; lay++) {
+ Int_t ncl=clusters[lay].GetEntriesFast(); if (ncl==0) continue;
+ while (ncl--) {
+ AliITSUClusterPix *pnt=
+ (AliITSUClusterPix*)clusters[lay].UncheckedAt(ncl);
+ Int_t l0=pnt->GetLabel(0);
+ if (l0>=np) {
+// cerr<<"Wrong label: "<<l0<<endl;
+ continue;
+ }
+ Int_t l1=pnt->GetLabel(1);
+ if (l1>=np) {
+// cerr<<"Wrong label: "<<l1<<endl;
+ continue;
+ }
+ Int_t l2=pnt->GetLabel(2);
+ if (l2>=np) {
+// cerr<<"Wrong label: "<<l2<<endl;
+ continue;
+ }
+ Int_t mask=1<<lay;
+ if (l0>=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; k<np; k++) {
+ if (good[k] != 0x7F) continue;
+
+ if (skip) continue; //No good primary vertex for this event
+
+ TParticle *p = (TParticle*)stack->Particle(k);
+ if (p == 0x0) {
+ cerr<<"Can not get particle "<<k<<endl;
+ continue;
+ }
+
+ if (p->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;
+}
+
+
--- /dev/null
+/****************************************************************************
+ * 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 <TMath.h>
+ #include <TError.h>
+ #include <Riostream.h>
+ #include <TH1F.h>
+ #include <TH2F.h>
+ #include <TTree.h>
+ #include <TParticle.h>
+ #include <TCanvas.h>
+ #include <TLine.h>
+ #include <TText.h>
+ #include <TBenchmark.h>
+ #include <TStyle.h>
+ #include <TFile.h>
+ #include <TROOT.h>
+ #include <TGeoGlobalMagField.h>
+
+ #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.4<p<0.5 GeV/c",50,0.,100.);
+
+ TH2F *hep=(TH2F*)gROOT->FindObject("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<<endl<<endl<<"********* Processing event number: "<<e<<"*******\n";
+
+ Int_t nentr=event->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; k<ngood; k++) {
+ AliTrackReference *ref=(AliTrackReference*)refs->UncheckedAt(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 (ptg<ptcutl) continue;
+ if (ptg>ptcuth) continue;
+
+ allselected++;
+
+ hgood->Fill(ptg);
+
+ AliESDtrack *track=0;
+ for (i=0; i<nentr; i++) {
+ track=event->GetTrack(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; mi<nentr; mi++) {
+ mitrack=event->GetTrack(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<<track_notfound[i]<<"\t";
+ if ((i%5)==4) cout<<"\n";
+ }
+ cout<<"\nList of fake tracks :\n";
+ for ( i = 0; i< itrack_fake; i++){
+ cout<<track_fake[i]<<"\t";
+ if ((i%5)==4) cout<<"\n";
+ }
+ cout<<"\nList of multiple found tracks :\n";
+ for ( i=0; i<itrack_multifound; i++) {
+ cout<<"id. "<<track_multifound[i]
+ <<" found - "<<track_multifound_n[i]<<"times\n";
+ }
+
+ cout<<"Number of found tracks ="<<nentr<<endl;
+ cout<<"Number of \"good\" tracks ="<<ngood<<endl;
+
+ refs->Clear();
+ }// ***** End of the loop over events
+
+ delete event;
+ delete esdTree;
+ ef->Close();
+
+ delete tpcTree;
+ refFile->Close();
+
+ Stat_t ng=hgood->GetEntries(), nf=hfound->GetEntries();
+ if (ng!=0) cout<<"\n\nIntegral efficiency is about "<<nf/ng*100.<<" %\n";
+ cout<<"Total number selected of \"good\" tracks ="<<allselected<<endl<<endl;
+ cout<<"Total number of found tracks ="<<allfound<<endl;
+ cout<<"Total number of \"good\" tracks ="<<allgood<<endl;
+ cout<<endl;
+
+ gStyle->SetOptStat(111110);
+ gStyle->SetOptFit(1);
+
+ TCanvas *c1=new TCanvas("c1","",0,0,700,850);
+
+ Int_t minc=33;
+
+ TPad *p1=new TPad("p1","",0,0.3,.5,.6); p1->Draw();
+ p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10);
+ hp->SetFillColor(4); hp->SetXTitle("(mrad)");
+ if (hp->GetEntries()<minc) hp->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()<minc) hl->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()<minc) hpt->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()<minc) hmpt->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()<minc) he->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 "<<ver<<" has been found !\n";
+ if (ver==1) tpcl->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; e<nev; e++) {
+ Int_t nt=0;
+ refs->Clear();
+
+ Int_t i;
+
+ rl->GetEvent(e); file->cd();
+
+ Int_t np = rl->GetHeader()->GetNtrack();
+ cout<<"Event "<<e<<" Number of particles: "<<np<<endl;
+
+ //******** Fill the "good" masks
+ Int_t *good=new Int_t[np]; for (i=0; i<np; i++) good[i]=0;
+ switch (ver) {
+ case 1:
+ /*
+ {
+ 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; n<nrows; n++) {
+ AliSegmentID *s=ca.LoadEntry(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]|=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; i<np; i++) count[i]=0;
+
+ Int_t sectors_by_rows=(Int_t)TD->GetEntries();
+ for (i=0; i<sectors_by_rows; i++) {
+ if (!TD->GetEvent(i)) continue;
+ Int_t sec,row;
+ digp->AdjustSectorRow(digits->GetID(),sec,row);
+ //cerr<<sec<<' '<<row<<'\r';
+ digits->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<np) count[idx0]+=1;
+ if (idx1>=0 && dig>=zero && idx1<np) count[idx1]+=1;
+ if (idx2>=0 && dig>=zero && idx2<np) count[idx2]+=1;
+ } while (digits->Next());
+ for (Int_t j=0; j<np; j++) {
+ if (count[j]>1) {
+ 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; i<np; i++) {
+ if ((good[i]&0x5000) != 0x5000)
+ if ((good[i]&0x2800) != 0x2800) continue;
+ if ((good[i]&0x7FF ) < good_number) continue;
+ TParticle *p = (TParticle*)stack->Particle(i);
+ if (p == 0x0) {
+ cerr<<"Can not get particle "<<i<<endl;
+ continue;
+ }
+ if (p->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; r<nr; r++) {
+ //cerr<<r<<' '<<nr<<'\r';
+ TR->GetEvent(r);
+
+ Int_t nref = tpcRefs->GetEntriesFast();
+ if (!nref) continue;
+ AliTrackReference *tpcRef= 0x0;
+ for (Int_t iref=0; iref<nref; ++iref) {
+ tpcRef = (AliTrackReference*)tpcRefs->UncheckedAt(iref);
+ if (tpcRef->DetectorId() == AliTrackReference::kTPC) break;
+ tpcRef = 0x0;
+ }
+
+ if (!tpcRef) continue;
+
+ Int_t j;
+ AliTrackReference *ref=0;
+ for (j=0; j<nt; j++) {
+ ref=(AliTrackReference *)refs->UncheckedAt(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;
+}
--- /dev/null
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+ #include <TMath.h>
+ #include <TError.h>
+ #include <Riostream.h>
+ #include <TH1F.h>
+ #include <TH2F.h>
+ #include <TTree.h>
+ #include <TParticle.h>
+ #include <TCanvas.h>
+ #include <TLine.h>
+ #include <TText.h>
+ #include <TBenchmark.h>
+ #include <TStyle.h>
+ #include <TFile.h>
+ #include <TROOT.h>
+ #include <TNtuple.h>
+ #include <TEllipse.h>
+ #include <TGeoManager.h>
+
+ #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; e<nev; e++) {
+ Int_t k;
+
+ rl->GetEvent(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; layer<nLayers; layer++) {
+ TClonesArray *ptr =
+ new(clusters+layer) TClonesArray("AliITSUClusterPix",1000);
+ Char_t bname[33];
+ sprintf(bname,"ITSRecPoints%d\0",layer);
+ branch[layer]=cTree->GetBranch(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; k<entr; k++) {
+ cTree->GetEvent(k);
+ for (Int_t lay=0; lay<nLayers; lay++) {
+ Int_t ncl=clusters[lay].GetEntriesFast(); if (ncl==0) continue;
+ while (ncl--) {
+ AliITSUClusterPix *pnt=
+ (AliITSUClusterPix*)clusters[lay].UncheckedAt(ncl);
+ Float_t g[3];
+ pnt->GetGlobalXYZ(g);
+ //pnt->GetLocalXYZ(g);
+ //cout<<g[0]<<' '<<g[1]<<' '<<g[2]<<endl;
+ nt->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;
+}
+
+
--- /dev/null
+#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;iEnt<hitTree->GetEntries();iEnt++){//entries loop degli hits
+ hitTree->GetEntry(iEnt);
+ int nh = hitList->GetEntries();
+ for(Int_t iHit=0; iHit<nh;iHit++){
+ AliITSUHit *pHit = (AliITSUHit*)hitList->At(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;ilr<nlr;ilr++) {
+ AliITSURecoLayer* lr = its->GetLayerActive(ilr);
+ TClonesArray* clr = lr->GetClusters();
+ int nClu = clr->GetEntries();
+ printf("Layer %d : %d clusters\n",ilr,nClu);
+ //
+ for (int icl=0;icl<nClu;icl++) {
+ AliITSUClusterPix *cl = (AliITSUClusterPix*)clr->At(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;ics<nclSn;ics++) {
+ AliITSUClusterPix* clusT = (AliITSUClusterPix*)lr->GetCluster(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;il<nLab;il++) {
+ TClonesArray* htArr = (TClonesArray*)arrMCTracks.At(labels[il]);
+ if (!htArr) {printf("did not find MChits for label %d ",labels[il]); cl->Print(); 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;i<nf;i++) {
+ TPaveStats *fnc = (TPaveStats*) lst->At(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;
+}
+
--- /dev/null
+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 ? evStart : nevTot-1;
+ if (evStart<0) evStart = 0;
+
+ Int_t lastEv = nev<0 ? nevTot : evStart+nev;
+ if (lastEv > 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; iEnt<hitTree->GetEntries(); iEnt++){
+ hitTree->GetEntry(iEnt);
+
+ for(Int_t iHit=0; iHit<hitList->GetEntries();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();
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff