--- /dev/null
+/****************************************************************************
+ * This macro is supposed to do reconstruction in the barrel ALICE trackers *
+ * (Make sure you have TPC digits and ITS hits before using this macro !!!) *
+ * It does the following steps (April 12, 2001): *
+ * 1) TPC cluster finding *
+ * 2) TPC track finding *
+ * 3) ITS cluster finding V2 (via fast points !) *
+ * 4) ITS track finding V2 *
+ * 5) ITS back track propagation V2 *
+ * 6) TPC back track propagation *
+ * (Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch) *
+ ****************************************************************************/
+
+#ifndef __CINT__
+ #include "alles.h"
+ #include "AliMagF.h"
+ #include "AliTPCtracker.h"
+
+ #include "AliITS.h"
+ #include "AliITSgeom.h"
+ #include "AliITSRecPoint.h"
+ #include "AliITSclusterV2.h"
+ #include "AliITSsimulationFastPoints.h"
+ #include "AliITStrackerV2.h"
+ #include<iostream.h>
+#endif
+
+Int_t TPCFindClusters(const Char_t *inname, const Char_t *outname, Int_t n);
+Int_t TPCFindTracks(const Char_t *inname, const Char_t *outname, Int_t n);
+Int_t TPCSortTracks(const Char_t *inname, const Char_t *inname2, const Char_t *outname, Int_t n);
+Int_t TPCPropagateBack(const Char_t *inname, const Char_t *outname);
+
+Int_t ITSFindClusters(const Char_t *inname, const Char_t *outname, Int_t n);
+Int_t ITSFindTracks(const Char_t *inname, const Char_t *inname2, const Char_t *outname, Int_t n);
+Int_t ITSPropagateBack(const Char_t *inname, const Char_t *outname);
+
+
+Int_t AliBarrelReconstructionV2(Int_t n=1) {
+ //const Char_t *TPCdigName="rfio:galice.root";
+ const Char_t *TPCdigName="galice.root";
+ const Char_t *TPCclsName="AliTPCclusters.root";
+ const Char_t *TPCtrkName="AliTPCtracks.root";
+ const Char_t *TPCtrkNameS="AliTPCtracksSorted.root";
+
+
+ //const Char_t *ITSdigName="rfio:galice.root";
+ const Char_t *ITSdigName="galice.root";
+ const Char_t *ITSclsName="AliITSclustersV2.root";
+ const Char_t *ITStrkName="AliITStracksV2.root";
+
+ AliKalmanTrack::SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor());
+
+
+// ********** Find TPC clusters *********** //
+ if (TPCFindClusters(TPCdigName,TPCclsName, n)) {
+ cerr<<"Failed to get TPC clusters !\n";
+ return 1;
+ }
+
+// ********** Find TPC tracks *********** //
+ if (TPCFindTracks(TPCclsName,TPCtrkName,n)) {
+ cerr<<"Failed to get TPC tracks !\n";
+ return 1;
+ }
+
+// ********** Sort and label TPC tracks *********** //
+ if (TPCSortTracks(TPCclsName,TPCtrkName,TPCtrkNameS,n)) {
+ cerr<<"Failed to sort TPC tracks !\n";
+ return 1;
+ }
+
+// ********** Find ITS clusters *********** //
+ if (ITSFindClusters(ITSdigName,ITSclsName,n)) {
+ cerr<<"Failed to get ITS clusters !\n";
+ return 1;
+ }
+
+// ********** Find ITS tracks *********** //
+ {
+ //TFile *clsFile=TFile::Open(ITSclsName);
+ if (ITSFindTracks(TPCtrkNameS,ITSclsName,ITStrkName,n)) {
+ cerr<<"Failed to get ITS tracks !\n";
+ return 1;
+ }
+ //clsFile->Close();
+ }
+ return 1;
+// ********** Back propagation of the ITS tracks *********** //
+ {TFile *clsFile=TFile::Open(ITSclsName);
+ if (ITSPropagateBack(ITStrkName,TPCtrkNameS)) {
+ cerr<<"Failed to propagate back the ITS tracks !\n";
+ return 1;
+ }
+ clsFile->Close();}
+
+
+// ********** Back propagation of the TPC tracks *********** //
+ {TFile *clsFile=TFile::Open(TPCclsName);
+ if (TPCPropagateBack(TPCtrkName,TPCtrkName)) {
+ cerr<<"Failed to propagate back the TPC tracks !\n";
+ return 1;
+ }
+ clsFile->Close();}
+
+ return 0;
+}
+
+
+Int_t TPCFindClusters(const Char_t *inname, const Char_t *outname, Int_t n) {
+ Int_t rc=0;
+ const Char_t *name="TPCFindClusters";
+ cerr<<'\n'<<name<<"...\n";
+ gBenchmark->Start(name);
+
+ TFile *out=TFile::Open(outname,"recreate");
+ TFile *in =TFile::Open(inname);
+
+ AliTPCParam *param=(AliTPCParam *)in->Get("75x40_100x60");
+ if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
+ AliTPCv2 tpc;
+ tpc.SetParam(param);
+
+ //tpc.Digits2Clusters(out); //MI change
+ cout<<"TPCFindClusters: nev ="<<n<<endl;
+ for (Int_t i=0;i<n;i++){
+ printf("Processing event %d\n",i);
+ tpc.Digits2Clusters(out,i);
+ // AliTPCclusterer::Digits2Clusters(dig, out, i);
+ }
+ in->Close();
+ out->Close();
+ gBenchmark->Stop(name);
+ gBenchmark->Show(name);
+
+ return rc;
+}
+
+Int_t TPCFindTracks(const Char_t *inname, const Char_t *outname, Int_t n) {
+ Int_t rc=0;
+ const Char_t *name="TPCFindTracks";
+ cerr<<'\n'<<name<<"...\n";
+ gBenchmark->Start(name);
+ TFile *out=TFile::Open(outname,"recreate");
+ TFile *in =TFile::Open(inname);
+ AliTPCParam *param=(AliTPCParam *)in->Get("75x40_100x60");
+ if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
+
+ //AliTPCtracker *tracker=new AliTPCtracker(param);
+ //rc=tracker->Clusters2Tracks(0,out);
+ //delete tracker;
+ cout<<"TPCFindTracks: nev ="<<n<<endl;
+
+ for (Int_t i=0;i<n;i++){
+ printf("Processing event %d\n",i);
+ AliTPCtracker *tracker = new AliTPCtracker(param,i);
+ //Int_t rc=
+ tracker->Clusters2Tracks(0,out);
+ delete tracker;
+ }
+
+ in->Close();
+ out->Close();
+ gBenchmark->Stop(name);
+ gBenchmark->Show(name);
+
+ return rc;
+}
+Int_t TPCSortTracks(const Char_t *inname, const Char_t *inname2, const Char_t * outname, Int_t eventn){
+ Int_t rc=0;
+ const Char_t *name="TPCSortTracks";
+ cerr<<'\n'<<name<<"...\n";
+ gBenchmark->Start(name);
+
+ TFile *out =TFile::Open(outname,"recreate");
+ TFile *in1 =TFile::Open(inname);
+ TFile *in2 =TFile::Open(inname2);
+
+
+ AliTPCParam *param=(AliTPCParam *)in1->Get("75x40_100x60");
+ if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
+
+ cout<<"TPCSortedtracks: nev ="<<eventn<<endl;
+
+
+ // loop over events
+ for (Int_t event=0;event<eventn; event++){
+ TObjArray tarray(10000);
+ AliTPCtrack *iotrack=0;
+ Int_t i;
+
+
+ in2->cd();
+ char tname[100];
+ sprintf(tname,"TreeT_TPC_%d",event);
+
+ TTree *tracktree=(TTree*)in2->Get(tname);
+ TBranch *tbranch=tracktree->GetBranch("tracks");
+ Int_t nentr=(Int_t)tracktree->GetEntries();
+ for (i=0; i<nentr; i++) {
+ iotrack=new AliTPCtrack;
+ tbranch->SetAddress(&iotrack);
+ tracktree->GetEvent(i);
+ tarray.AddLast(iotrack);
+ }
+ tarray.Sort();
+ // in2->Close();
+
+ //assign thacks GEANT labels
+ in1->cd();
+ AliTPCtracker *tracker = new AliTPCtracker(param,event);
+ tracker->LoadInnerSectors();
+ tracker->LoadOuterSectors();
+ for (i=0; i<nentr; i++) {
+ iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
+ tracker->CookLabel(iotrack,0.1);
+ }
+ delete tracker;
+ //in->Close();
+ //end of GEANT label assignment
+
+ tracktree=new TTree(tname,"Tree with TPC tracks");
+ tracktree->Branch("tracks","AliTPCtrack",&iotrack,32000,0);
+ for (i=0; i<nentr; i++) {
+ iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
+ tracktree->Fill();
+ }
+ out->cd();
+ tracktree->Write();
+ }
+
+ out->Close();
+ in2->Close();
+ in1->Close();
+
+ gBenchmark->Stop(name);
+ gBenchmark->Show(name);
+
+ return rc;
+}
+
+Int_t ITSFindClusters(const Char_t *inname, const Char_t *outname, Int_t n) {
+ Int_t rc=0;
+ const Char_t *name="ITSFindClusters";
+ cerr<<'\n'<<name<<"...\n";
+ gBenchmark->Start(name);
+ TFile *out=TFile::Open(outname,"recreate");
+ TFile *in =TFile::Open(inname,"update");
+
+ if (!(gAlice=(AliRun*)in->Get("gAlice"))) {
+ cerr<<"Can't get gAlice !\n";
+ return 1;
+ }
+
+ AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
+ if (!ITS) { cerr<<"Can't get the ITS !\n"; return 1;}
+ AliITSgeom *geom=ITS->GetITSgeom();
+ out->cd();
+ geom->Write();
+
+ cout<<"ITSFindClusters: nev ="<<n<<endl;
+ Int_t ev=0;
+ for (ev = 0; ev<n; ev++){
+ in->cd(); // !!!! MI directory must point to galice. - othervise problem with Tree -connection
+ gAlice->GetEvent(ev);
+ //gAlice->TreeR()->Reset(); //reset reconstructed tree
+
+
+ TTree *pTree=gAlice->TreeR();
+ if (!pTree){
+ gAlice->MakeTree("R");
+ pTree = gAlice->TreeR();
+ }
+ TBranch *branch=pTree->GetBranch("ITSRecPoints");
+ /*
+ if (!branch) {
+ //if not reconstructed ITS branch do reconstruction
+ ITS->MakeBranch("R",0);
+
+ //////////////// Taken from ITSHitsToFastPoints.C ///////////////////////
+ AliITSsimulationFastPoints *sim = new AliITSsimulationFastPoints();
+ for (Int_t i=0;i<3;i++) { ITS->SetSimulationModel(i,sim); }
+ Int_t nsignal=25;
+ Int_t size=-1;
+ Int_t bgr_ev=Int_t(ev/nsignal);
+ ITS->HitsToFastRecPoints(ev,bgr_ev,size," ","All"," ");
+ //////////////////////////////////////////////////////////////////////////
+ gAlice->GetEvent(ev); //MI comment - in HitsToFast... they reset treeR to 0
+ //they overwrite full reconstructed event ???? ... so lets connect TreeR one more
+ //time
+ }
+ */
+
+
+ out->cd();
+ TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
+ char cname[100];
+ sprintf(cname,"TreeC_ITS_%d",ev);
+
+ TTree *cTree=new TTree(cname,"ITS clusters");
+ cTree->Branch("Clusters",&clusters);
+
+ pTree=gAlice->TreeR();
+ if (!pTree) { cerr<<"Can't get TreeR !\n"; return 1; }
+ branch=pTree->GetBranch("ITSRecPoints");
+ if (!branch) { cerr<<"Can't get ITSRecPoints branch !\n"; return 1;}
+ TClonesArray *points=new TClonesArray("AliITSRecPoint",10000);
+ branch->SetAddress(&points);
+
+ TClonesArray &cl=*clusters;
+ Int_t nclusters=0;
+ Int_t nentr=(Int_t)pTree->GetEntries();
+ Int_t lab[6];
+ Float_t lp[5];
+ AliITSgeom *geom=ITS->GetITSgeom();
+
+ cout<<"!!!! nentr ="<<nentr<<endl;
+ for (Int_t i=0; i<nentr; i++) {
+ if (!pTree->GetEvent(i)) {cTree->Fill(); continue;}
+ Int_t lay,lad,det; geom->GetModuleId(i,lay,lad,det);
+ Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift);
+ Double_t rot[9]; geom->GetRotMatrix(lay,lad,det,rot);
+ Double_t yshift = x*rot[0] + y*rot[1];
+ Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1);
+ Int_t ncl=points->GetEntriesFast();
+ nclusters+=ncl;
+
+ Float_t kmip=1.;
+ if(lay<5&&lay>2){kmip=280.;}; // b.b.
+ if(lay<7&&lay>4){kmip=38.;};
+
+
+ for (Int_t j=0; j<ncl; j++) {
+ AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j);
+ lp[0]=-p->GetX()-yshift; if (lay==1) lp[0]=-lp[0];
+ lp[1]=p->GetZ()+zshift;
+ lp[2]=p->GetSigmaX2();
+ lp[3]=p->GetSigmaZ2();
+ lp[4]=p->GetQ(); lp[4]/=kmip; // b.b.
+ if(ncl==11)cout<<"mod,cl,Q ="<<i<<","<<j<<","<<lp[4]<<endl;
+ lab[0]=p->GetLabel(0);lab[1]=p->GetLabel(1);lab[2]=p->GetLabel(2);
+ lab[3]=ndet;
+
+ Int_t label=lab[0];
+ if(label>=0) { // b.b.
+ TParticle *part=(TParticle*)gAlice->Particle(label);
+ label=-3;
+ while (part->P() < 0.005) {
+ Int_t m=part->GetFirstMother();
+ if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
+ label=m;
+ part=(TParticle*)gAlice->Particle(label);
+ }
+ }
+ if (lab[1]<0) lab[1]=label;
+ else if (lab[2]<0) lab[2]=label;
+ else cerr<<"No empty labels !\n";
+
+ new(cl[j]) AliITSclusterV2(lab,lp);
+ }
+ cTree->Fill(); clusters->Delete();
+ points->Delete();
+ }
+ cTree->Write();
+ cerr<<"Number of clusters: "<<nclusters<<endl;
+ delete cTree; delete clusters; delete points;
+
+ }
+
+ delete gAlice; gAlice=0;
+ in->Close();
+ out->Close();
+ gBenchmark->Stop(name);
+ gBenchmark->Show(name);
+
+ return rc;
+}
+
+Int_t ITSFindTracks(const Char_t * inname, const Char_t *inname2, const Char_t *outname, Int_t n) {
+ Int_t rc=0;
+ const Char_t *name="ITSFindTracks";
+ cerr<<'\n'<<name<<"...\n";
+ gBenchmark->Start(name);
+
+
+ TFile *out=TFile::Open(outname,"recreate");
+ TFile *in =TFile::Open(inname);
+ TFile *in2 =TFile::Open(inname2);
+
+ AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
+ if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
+
+ cout<<"ITSFindtracks: nev ="<<n<<endl;
+
+ for (Int_t i=0;i<n;i++){
+ AliITStrackerV2 tracker(geom,i);
+ rc=tracker.Clusters2Tracks(in,out);
+ }
+
+ in->Close();
+ in2->Close();
+ out->Close();
+
+ gBenchmark->Stop(name);
+ gBenchmark->Show(name);
+
+ return rc;
+}
+
+Int_t ITSPropagateBack(const Char_t *inname, const Char_t *outname) {
+
+ Int_t rc=0;
+ /*
+ const Char_t *name="ITSPropagateBack";
+ cerr<<'\n'<<name<<"...\n";
+ gBenchmark->Start(name);
+
+ AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
+ if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
+ AliITStrackerV2 tracker(geom);
+
+ TFile *out=TFile::Open(outname,"update");
+ TFile *in =TFile::Open(inname);
+ rc=tracker.PropagateBack(in,out);
+ in->Close();
+ out->Close();
+
+ gBenchmark->Stop(name);
+ gBenchmark->Show(name);
+ */
+ return rc;
+}
+
+Int_t TPCPropagateBack(const Char_t *inname, const Char_t *outname) {
+ Int_t rc=0;
+ const Char_t *name="TPCPropagateBack";
+ cerr<<'\n'<<name<<"...\n";
+ gBenchmark->Start(name);
+
+ AliTPCParam *param=(AliTPCParam *)gFile->Get("75x40_100x60");
+ if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
+ AliTPCtracker *tracker=new AliTPCtracker(param);
+
+ TFile *out=TFile::Open(outname,"update");
+ TFile *in =TFile::Open(inname);
+ rc=tracker->PropagateBack(in,out);
+ delete tracker;
+ in->Close();
+ out->Close();
+
+ gBenchmark->Stop(name);
+ gBenchmark->Show(name);
+
+ return rc;
+}
+
+
+
+
+
#include "AliITSPid.h"
#include "TMath.h"
+#include "AliITSIOTrack.h"
#include <iostream.h>
ClassImp(AliITSPid)
if(fcorr<=0.1)fcorr=0.1;
return qtot/fcorr;
}
-
//-----------------------------------------------------------
#include "vector.h"
#include <algorithm>
for(Int_t i=0;i<nl;i++){q(i+1)=vf[i];} this->SetVec(track,q);
return (q(6));
}
+//........................................................
+Float_t AliITSPid::qtrm(Float_t qarr[6],Int_t narr)
+{
+ Float_t q[6],qm;
+ Int_t nl,ml;
+ narr>6?ml=6:ml=narr;
+ nl=0; for(Int_t i=0;i<ml;i++){if(qarr[i]>0.){q[i]=qarr[i];nl++;}}
+ if(nl==0)return 0.;
+ vector<float> vf(nl);
+ switch(nl){
+ case 1:qm=q[0]; break;
+ case 2:qm=(q[0]+q[1])/2.; break;
+ default:
+ for(Int_t i=0;i<nl;i++){vf[i]=q[i];}
+ sort(vf.begin(),vf.end());
+ qm= (vf[0]+vf[1])/2.;
+ break;
+ }
+ return qm;
+}
//-----------------------------------------------------------
+//inline
Int_t AliITSPid::wpik(Int_t nc,Float_t q)
{
return pion();
return pion();
}
//-----------------------------------------------------------
+//inline
Int_t AliITSPid::wpikp(Int_t nc,Float_t q)
{
return pion();
return 0;
}
//-----------------------------------------------------------
+Int_t AliITSPid::GetPcode(AliTPCtrack *track)
+{
+ Double_t xk,par[5]; track->GetExternalParameters(xk,par);
+ Float_t phi=TMath::ASin(par[2]) + track->GetAlpha();
+ if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
+ if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
+ Float_t lam=TMath::ATan(par[3]);
+ Float_t pt_1=TMath::Abs(par[4]);
+ Float_t mom=1./(pt_1*TMath::Cos(lam));
+ Float_t dedx=track->GetdEdx();
+ Int_t pcode=GetPcode(dedx/40.,mom);
+ cout<<"TPCtrack dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
+ return pcode?pcode:211;
+ }
+//------------------------------------------------------------
+/*
+Int_t AliITSPid::GetPcode(AliITSIOTrack *track)
+{
+ Double_t px,py,pz;
+ px=track->GetPx();
+ py=track->GetPy();
+ pz=track->GetPz();
+ Float_t mom=TMath::Sqrt(px*px+py*py+pz*pz);
+//???????????????????
+ // Float_t dedx=1.0;
+ Float_t dedx=track->GetdEdx();
+//???????????????????
+ Int_t pcode=GetPcode(dedx,mom);
+ cout<<"ITSV1 dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
+return pcode?pcode:211;
+}
+*/
+//-----------------------------------------------------------
+Int_t AliITSPid::GetPcode(AliITStrackV2 *track)
+{
+ if(track==0)return 0;
+ //track->Propagate(track->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848);
+ track->PropagateTo(3.,0.0028,65.19);
+
+ track->PropagateToVertex();
+ Double_t xk,par[5]; track->GetExternalParameters(xk,par);
+ Float_t lam=TMath::ATan(par[3]);
+ Float_t pt_1=TMath::Abs(par[4]);
+ Float_t mom=0.;
+ if( (pt_1*TMath::Cos(lam))!=0. ){ mom=1./(pt_1*TMath::Cos(lam)); }else{mom=0.;};
+ Float_t dedx=track->GetdEdx();
+ //cout<<"lam,pt_1,mom,dedx="<<lam<<","<<pt_1<<","<<mom<<","<<dedx<<endl;
+ Int_t pcode=GetPcode(dedx,mom);
+ //cout<<"ITS V2 dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
+return pcode?pcode:211;
+}
+//-----------------------------------------------------------
Int_t AliITSPid::GetPcode(Float_t q,Float_t pm)
{
fWpi=fWk=fWp=-1.; fPcode=0;
#include <TObject.h>
#include <TClonesArray.h>
#include <TVector.h>
+#include "../TPC/AliTPCtrack.h"
+#include "AliITSIOTrack.h"
+#include "AliITStrackV2.h"
#include <assert.h>
//#include <stl.h>
//___________________________________________________________________________
TVector* GetVec(Int_t track);
Int_t GetPcode(TClonesArray*,Float_t);
Int_t GetPcode(Float_t,Float_t);
+ Int_t GetPcode(AliTPCtrack*); // For PID TPC
+// Int_t GetPcode(AliITSIOTrack*); // For PID ITS tracking V1
+ Int_t GetPcode(AliITStrackV2*); // For PID ITS tracking V2
void SetCut(Int_t,Float_t,Float_t,Float_t,
Float_t,Float_t,Float_t,Float_t);
Float_t GetWpi(){return fWpi;}
Float_t qcorr(Float_t);
int qcomp(Float_t* qa,Float_t* qb){return qa[0]>qb[0]?1:0;}
Float_t qtrm(Int_t track);
+ Float_t qtrm(Float_t qarr[6],Int_t narr);
Int_t wpik(Int_t,Float_t);
Int_t wpikp(Int_t,Float_t);
Int_t pion(){return fWpi=1.,fPcode=211;}
--- /dev/null
+#!/bin/sh
+# This script processes the following steps for n events
+# (use the parameter n):
+# - the TPC and ITS slow simulation (hits+digits+clusters),
+# - the TPC+ITS tracking,
+# - the TPC and ITS PID
+# (Dubnna version)
+#
+# delete eventual old files from the last run
+rm -f *.root
+#
+# run the hit generation
+aliroot -q -b "$ALICE_ROOT/macros/grun.C($1)"
+#
+# digitize TPC
+aliroot -q -b "$ALICE_ROOT/TPC/AliTPCHits2Digits.C($1)"
+#
+# digitize ITS
+aliroot -q -b "$ALICE_ROOT/ITS/AliITSHits2DigitsDefault.C"
+#
+# create reconstructed points for the ITS
+aliroot -q -b "$ALICE_ROOT/ITS/ITSDigitsToClusters.C(0,$[$1-1])"
+#
+# do the TPC+ITS tracking
+aliroot -q -b "$ALICE_ROOT/ITS/AliBarrelReconstructionV2.C($1)"
+#
+# Do the PID procedure for the ITS. The output file AliITStrackV2Pid.root is
+# created with the information for each track:
+# the ITS ADC trancated mean signal, the particle reconstructed momentum,
+# the probable weights for the different particle kinds (electron, pion, kaon,
+# proton). These weights are under study now.
+#
+aliroot -q -b "$ALICE_ROOT/ITS/save_pidV2.C(0,$[$1-1])"
+#
+# This last line checks the ITS PID only, i.e. creates and draws
+#the dEdx-momentum plot (comment if it's not necessary)
+aliroot "$ALICE_ROOT/ITS/scan_pidV2.C(0,$[$1-1])"
+
+
#!/bin/sh
-# delete eventual old files from the last run
-./DeleteOldV2Files
-# run the hit generation
-aliroot -q -b "$ALICE_ROOT/macros/grun.C"
-# digitize TPC and prepare TPC tracks for matching with the ITS
-aliroot -q -b "$ALICE_ROOT/TPC/AliTPCtest.C"
-# digitize ITS
-aliroot -q -b "$ALICE_ROOT/ITS/AliITSHitsToDigitsDefault.C"
-# create reconstructed points for the ITS
-aliroot -q -b "$ALICE_ROOT/ITS/ITSDigitsToClusters.C"
-#prepare slow or fast recpoints for the tracking
-aliroot -q -b "$ALICE_ROOT/ITS/"AliITSFindClustersV2.C\(\'$2\'\)
-# ITS tracking
-aliroot -q -b "$ALICE_ROOT/ITS/SetConvConst.C" "$ALICE_ROOT/ITS/AliITSFindTracksV2.C"
-# do the PID
-aliroot -q -b "$ALICE_ROOT/ITS/SetConvConst.C" "$ALICE_ROOT/ITS/save_pid.C"
-aliroot "$ALICE_ROOT/ITS/ITSPIDtest.C"
-#
-# After all of the above the PID menu will appear. To check the PID package
-# put the buttons:
-#
-#1. "fill tab_tr" and "dEdX-P plot" (or dEdX-P pions, ... kaons, ... elect,
-# ... prot) and the dEdX versus momentum taken from the tracking for all
-# particle (pions, kaons, electrons, protons) will be drown.
-#
-#2. "dEdX.C" and the same other ones as for the point 1. In this case the
-# same plots but fot the generated particle momentum will be drown.
-#
-# The output file AliITStrackV2Pid.root will be written with the information
-# for each track:
+# This script processes the ITS PID for n events
+# (use the parameter n):
+# - the AliTPCtest.C and AliITStest.C should be processed befor
+# (Dubnna version)
+#
+# Do the PID procedure for the ITS. The output file AliITStrackV2Pid.root is
+# created with the information for each track:
# the ITS ADC trancated mean signal, the particle reconstructed momentum,
# the probable weights for the different particle kinds (electron, pion, kaon,
# proton). These weights are under study now.
-
-
-
-
-
+#
+rm -f AliITStrackingV2Pid.root
+#
+aliroot -q -b "$ALICE_ROOT/ITS/save_pidV2.C(0,$[$1-1])"
+#
+# This last line checks the ITS PID only, i.e. creates and draws
+#the dEdx-momentum plot (comment if it's not necessary)
+aliroot "$ALICE_ROOT/ITS/scan_pidV2.C(0,$[$1-1])"
//qplotP->GetXaxis()->SetTitleSize(0.05);
//qplotP->GetYaxis()->SetTitleSize(0.05);
+ //gPad->SetFillColor(kWhite); // b.b.
gStyle->SetOptStat(0);
qplotP->SetXTitle("(Mev/c)");
qplotP->SetYTitle("(mips)");
cerr<<"Track "<<lab<<" was not found !\n";
continue;
}
- track->Propagate(track->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848);
+ //track->Propagate(track->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848);
+ track->PropagateTo(3.,0.0028,65.19);
+
track->PropagateToVertex();
if (lab==tlab) hfound->Fill(ptg);
--- /dev/null
+////////////////////////////////////////////////////////////
+// Test macro for AliITSPid class and tracking version V2 //
+////////////////////////////////////////////////////////////
+void
+save_pidV2(Int_t evNumber1=0,Int_t evNumber2=0) {
+ const char *filename="AliITStracksV2.root";
+
+ if (gClassTable->GetID("AliRun") < 0) {
+ gROOT->LoadMacro("loadlibs.C"); loadlibs();
+ } else { delete gAlice; gAlice=0;
+ }
+
+ TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
+ if (!file) file = new TFile(filename);
+
+TFile *fpid = new TFile("AliITStracksV2Pid.root","recreate");
+ Float_t factor=0;
+ if(gAlice!=0)factor=gAlice->Field()->Factor();
+ if(factor==0.)factor=1.;
+ AliKalmanTrack::SetConvConst(100/0.299792458/0.2/factor);
+
+//
+// Loop over events
+//
+ for (int nev=0; nev<= evNumber2; nev++) {
+ char tname[30];
+ sprintf(tname,"TreeT_ITS_%d;1",nev);
+ TTree *tracktree=(TTree*)file->Get(tname);
+ TBranch *tbranch=tracktree->GetBranch("tracks");
+
+ Int_t nentr=tracktree->GetEntries();
+ cout<<"Found "<<nentr<<" ITS tracks in event No "<<nev<<endl;
+
+ char tpidname[30];
+ sprintf(tpidname,"TreeT%d",nev);
+ AliITStrackV2Pid pidtmp;
+ TTree itspidTree(tpidname,"Tree with PID");
+ AliITStrackV2Pid *outpid=&pidtmp;
+ itspidTree.Branch("pids","AliITStrackV2Pid",&outpid,32000,1);
+ AliITSPid pid;
+
+ AliITStrackV2 *iotrack=0;
+ for (Int_t i=0; i<nentr; i++) {
+ AliITStrackV2 *iotrack=new AliITStrackV2;
+ tbranch->SetAddress(&iotrack);
+ tracktree->GetEvent(i);
+ Int_t pcode=pid->GetPcode(iotrack);
+ Float_t signal=iotrack->GetdEdx();
+ //iotrack->Propagate(iotrack->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848);
+ iotrack->PropagateTo(3.,0.0028,65.19);
+
+ iotrack->PropagateToVertex();
+ Double_t xk,par[5]; iotrack->GetExternalParameters(xk,par);
+ Float_t lam=TMath::ATan(par[3]);
+ Float_t pt_1=TMath::Abs(par[4]);
+ Float_t mom=0.;
+ if( (pt_1*TMath::Cos(lam))!=0. ){ mom=1./(pt_1*TMath::Cos(lam)); }else{mom=0.;};
+ pidtmp.fPcode=pcode;
+ pidtmp.fSignal=signal;
+ pidtmp.fMom=mom;
+ itspidTree.Fill();
+ delete iotrack;
+ }// End for i (tracks)
+ fpid->Write();
+ }// End for nev (events)
+ fpid->Close();
+ file->Close();
+ cout<<"File AliITStracksV2Pid.root written"<<endl;
+ return;
+///////////////////////////////////////////////////////
+}
+
+
+
+
+
+
+
+
--- /dev/null
+///////////////////////////////////////////////////////////
+// Test macro for AliITStracksV1Pid.root file //
+// JINR Dubna Jan 2002 //
+///////////////////////////////////////////////////////////
+void
+scan_pidV2(Int_t evNumber1=0,Int_t evNumber2=0) {
+ //................. Prepare histogramms ................
+ TH2F *qplot = new TH2F("Qtrm","Qtrm vs Pmom",100,0,1300,100,0,13);
+ TH2F *qplotP= new TH2F("Qtrm","Qtrm vs Pmom",100,0,1300,100,0,13);
+ TH2F *qplotKa= new TH2F("Qtrm","Qtrm vs Pmom",100,0,1300,100,0,13);
+ TH2F *qplotPi= new TH2F("Qtrm","Qtrm vs Pmom",100,0,1300,100,0,13);
+ TH2F *qplotE= new TH2F("Qtrm","Qtrm vs Pmom",100,0,1300,100,0,13);
+ qplotP.SetMarkerStyle(8); qplotP.SetMarkerColor(kBlack); qplotP.SetMarkerSize(.3);
+ qplotKa.SetMarkerStyle(8); qplotKa.SetMarkerColor(kRed); qplotKa.SetMarkerSize(.3);
+ qplotPi.SetMarkerStyle(8); qplotPi.SetMarkerColor(kBlue); qplotPi.SetMarkerSize(.3);
+ qplotE.SetMarkerStyle(8); qplotE.SetMarkerColor(kGreen); qplotE.SetMarkerSize(.3);
+ //......................................................
+ TH1F *signal_mip = new TH1F("signal_mip","Signal (mips) for track",100,0.,15.);
+
+TFile *fpid = new TFile("AliITStracksV2Pid.root","read");
+fpid->ls();
+//
+// Loop over events
+//
+for (int nev=0; nev<= evNumber2; nev++) {
+ char tpidname[30];
+ sprintf(tpidname,"TreeT%d",nev);
+ TTree *tracktree=(TTree*)fpid->Get(tpidname);
+ TBranch *tbranch=tracktree->GetBranch("pids");
+
+ Int_t nentr=tracktree->GetEntries();
+ cout<<"Found PID for "<<nentr<<" ITS V1 tracks on "<<tpidname<<endl;
+
+ AliITStrackV2Pid *iopid=0;
+for(Int_t ii=0;ii<nentr;ii++)
+ {
+ AliITStrackV2Pid *iopid=new AliITStrackV2Pid;
+ tbranch->SetAddress(&iopid);
+ tracktree->GetEvent(ii);
+
+ signal_mip->Fill(iopid->fSignal);
+
+ if(iopid->fPcode ==2212)qplotP.Fill(1000*iopid->fMom,iopid->fSignal);
+ if(iopid->fPcode == 321)qplotKa.Fill(1000*iopid->fMom,iopid->fSignal );
+ if(iopid->fPcode == 211)qplotPi.Fill(1000*iopid->fMom,iopid->fSignal );
+ if(iopid->fPcode == 11)qplotE.Fill(1000*iopid->fMom,iopid->fSignal );
+
+ cout<<"PID pcode,fsignal,fmom= "
+ <<iopid->fPcode<<","<<iopid->fSignal<<","<<iopid->fMom<<endl;
+ delete iopid;
+ }// Enf for ii (tracks)
+ }// End for nev (events)
+ fpid->Close();
+ //...................... Draw histogramms .................
+ TCanvas *c1 = new TCanvas("PID_test","Scan PID ",200,10,900,700);
+ c1->Divide(2,1);
+ //.........................................................
+ c1->cd(1); gPad->SetFillColor(33);
+ signal_mip->Draw();
+
+ c1->cd(2); //gPad->SetFillColor(33);
+ qplot->Draw();
+ qplotP.Draw("same"); qplotKa.Draw("same"); qplotPi.Draw("same"); qplotE.Draw("same");
+ AliITSPid *pid =new AliITSPid(100);
+ c1->Range(0,0,1300,10);
+ gStyle->SetLineColor(kRed);
+ gStyle->SetLineWidth(2);
+ TLine *lj[3],*lk[3];
+ for(Int_t j=0;j<3;j++){
+ Float_t x1,x2,y1,y2,xx1,xx2,yy1,yy2;
+ x1=pid->cut[j+1][0]; x2=pid->cut[j+2][0];
+ y1=y2=pid->cut[j+2][2];
+ lj[j]=new TLine(1000*x1,y1,1000*x2,y2); lj[j]->Draw();
+ if(j==0){yy1=10.;}else{yy1=lj[j-1]->GetY1();}
+ yy2=lj[j]->GetY1();
+ xx1=xx2=x1;
+ lk[j]=new TLine(1000*xx1,yy1,1000*xx2,yy2); lk[j]->Draw();
+ }
+ //Draw pions-kaons cuts.
+ TLine *mj[7],*mk[7];
+ for(Int_t j=0;j<7;j++){
+ Float_t x1,x2,y1,y2,xx1,xx2,yy1,yy2;
+ x1=pid->cut[j+2][0]; x2=pid->cut[j+3][0];
+ y1=y2=pid->cut[j+3][5];
+ mj[j]=new TLine(1000*x1,y1,1000*x2,y2); mj[j]->Draw();
+ if(j==0){yy1=10.;}else{yy1=mj[j-1]->GetY1();}
+ yy2=mj[j]->GetY1();
+ xx1=xx2=x1;
+ mk[j]=new TLine(1000*xx1,yy1,1000*xx2,yy2); mk[j]->Draw();
+ }
+ cout<<"End of file AliITStracksV1Pid.root "<<endl;
+ return;
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff