-/****************************************************************************
- * Very important, delicate and rather obscure macro. *
- * *
- * Creates list of "trackable" tracks, *
- * sorts tracks for matching with the ITS, *
- * calculates efficiency, resolutions etc. *
- * *
- * Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch *
- * with several nice improvements by: M.Ivanov, GSI, m.ivanov@gsi.de *
- ****************************************************************************/
+/// \file AliTPCComparison.C
+/// \brief Very important, delicate and rather obscure macro
+///
+/// Creates list of "trackable" tracks, sorts tracks for matching with the ITS,
+/// calculates efficiency, resolutions etc.
+///
+/// There is a possibility to run this macro over several events.
+///
+/// \author I.Belikov, CERN, Jouri.Belikov@cern.ch, M.Ivanov, GSI, m.ivanov@gsi.de
#if !defined(__CINT__) || defined(__MAKECINT__)
#include<fstream.h>
#include<TTree.h>
#include <AliMagF.h>
#include <AliRun.h>
-#include <AliTPCtracker.h>
#include <alles.h>
#endif
/***********************************************************************/
TFile *inkin = TFile::Open("rfio:galice.root");
- // if(gAlice)delete gAlice; COMMENTED BECAUSE OF A BUG (IN COMPILED MODE)
+// \file AliTPCComparison2.C
+// if(gAlice)delete gAlice; COMMENTED BECAUSE OF A BUG (IN COMPILED MODE)
+
gAlice = (AliRun*)inkin->Get("gAlice");
cout<<"AliRun object found on file "<<gAlice<<endl;
- Float_t fifac=gAlice->Field()->Factor();
- AliKalmanTrack::SetConvConst(100/0.299792458/0.2/fifac);
- cout<<fifac*0.2<<" T Magnetic field is used\n";
+ AliKalmanTrack::SetFieldMap(gAlice->Field());
inkin->Close();
/*
delete gAlice; COMMENTED BECAUSE OF A BUG IN COMPILED MODE
if(kOLD){
cf=TFile::Open("AliTPCclusters.root");
if (!cf->IsOpen()) {cerr<<"Can't open AliTPCclusters.root !\n"; return 1;}
- digp= (AliTPCParam*)cf->Get("75x40_100x60");
+ digp= (AliTPCParam*)cf->Get("75x40_100x60_150x60");
if (!digp) { cerr<<"TPC parameters have not been found !\n"; return 2; }
}
- ///////////
- AliTPCtracker *tracker =0;
+
TObjArray tarray(MAX);
AliTPCtrack *iotrack=0;
Int_t nentr= 0;
cout<<"================================================"<<endl;
cout<<"Processing event "<<event<<endl;
cout<<"================================================"<<endl;
- if(kOLD){
- cf->cd();
- tracker = new AliTPCtracker(digp,event);
- tracker->LoadInnerSectors();
- tracker->LoadOuterSectors();
- }
char tname[100];
if (eventn==-1) {
iotrack=new AliTPCtrack;
tbranch->SetAddress(&iotrack);
tracktree->GetEvent(i);
- if(kOLD)tracker->CookLabel(iotrack,0.1);
tarray.AddLast(iotrack);
}
eventptr[event+1] = nentr; //store end of the event
- delete tracker;
delete tracktree;
}
tf->Close();
gt[ngood].x >>gt[ngood].y >>gt[ngood].z) {
ngood++;
cerr<<ngood<<"\r";
+ //cout<<ngood<<"\r";
if (ngood==MAX) {
cerr<<"Too many good tracks !\n";
break;
} else {
cerr<<"Marking good tracks (this will take a while)...\n";
ngood=good_tracks(gt,45000,firstev,eventn);
+ printf("Goood %d\n", ngood);
ofstream out("good_tracks_tpc");
- ofstream out2("good_tracks_tpc_par");
-
if (out) {
cout<<"File good_tracks_tpc opened\n";
for (Int_t ngd=0; ngd<ngood; ngd++) {
- Float_t pt = TMath::Sqrt(gt[ngd].px*gt[ngd].px+gt[ngd].py*gt[ngd].py);
- Float_t angle = gt[ngd].pz/pt;
out<<gt[ngd].fEventN<<' '<<gt[ngd].lab<<' '<<gt[ngd].code<<' '<<
gt[ngd].px<<' '<<gt[ngd].py<<' '<<gt[ngd].pz<<' '<<
gt[ngd].x <<' '<<gt[ngd].y <<' '<<gt[ngd].z <<endl;
- out2<<gt[ngd].fEventN<<"\t"<<gt[ngd].lab<<"\t"<<gt[ngd].code<<"\t"<<
- pt<<"\t"<<angle<<"\t"<<endl;
}
} else cerr<<"Can not open file (good_tracks_tpc) !\n";
+ out<<flush;
out.close();
+
+ ofstream out2("good_tracks_tpc_par");
+
+ if (out2) {
+ //cout<<"File good_tracks_tpc opened\n";
+ for (Int_t ngd=0; ngd<ngood; ngd++) {
+ Float_t pt = TMath::Sqrt(gt[ngd].px*gt[ngd].px+gt[ngd].py*gt[ngd].py);
+ Float_t angle = 0;
+ if (TMath::Abs(pt)>0.01) angle = gt[ngd].pz/pt;
+ out2<<gt[ngd].fEventN<<"\t"<<gt[ngd].lab<<"\t"<<gt[ngd].code<<"\t"<<
+ pt<<"\t"<<angle<<"\t"<<endl;
+ }
+ } else cerr<<"Can not open file (good_tracks_tpc) !\n";
+ out2<<flush;
out2.close();
+
}
cerr<<"Number of good tracks : "<<ngood<<endl;
+ cout<<"Number of good tracks : "<<ngood<<endl;
if(ngood==0)return 5;
- TH1F *hp=new TH1F("hp","PHI resolution",50,-200.,200.); hp->SetFillColor(4);
- TH1F *hl=new TH1F("hl","LAMBDA resolution",50,-200,200);hl->SetFillColor(4);
+ TH1F *hp=new TH1F("hp","PHI resolution",50,-20.,20.); hp->SetFillColor(4);
+ TH1F *hl=new TH1F("hl","LAMBDA resolution",50,-20,20);hl->SetFillColor(4);
TH1F *hpt=new TH1F("hpt","Relative Pt resolution",30,-10.,10.);
hpt->SetFillColor(2);
TH1F *hmpt=new TH1F("hmpt","Relative Pt resolution (pt>4GeV/c)",30,-60,60);
//
Double_t xk=gt[ngood].x;
- printf("Track =%p\n",track);
+ if (!track) continue;
+ // printf("Track =%p\n",track);
track->PropagateTo(xk);
if (lab==tlab) hfound->Fill(ptg);
if (TMath::Abs(gt[ngood].code)==11 && ptg>4.) {
hmpt->Fill((pt_1 - 1/ptg)/(1/ptg)*100.);
- } else {
+ }
+ //else
+ {
Float_t phig=TMath::ATan2(gt[ngood].py,gt[ngood].px);
hp->Fill((phi - phig)*1000.);
Stat_t ng=hgood->GetEntries(), nf=hfound->GetEntries();
if (ng!=0) cerr<<"\n\nIntegral efficiency is about "<<nf/ng*100.<<" %\n";
+ if (ng!=0) cout<<"\n\nIntegral efficiency is about "<<nf/ng*100.<<" %\n";
cout<<"Total number of found tracks ="<<nentr<<endl;
cout<<"Total number of \"good\" tracks ="
<<mingood<<" (selected for comparison: "<<ng<<')'<<endl<<endl;
Int_t good_tracks(GoodTrackTPC *gt, Int_t max, Int_t firstev, Int_t eventn) {
- //eventn - number of events in file
+ /// eventn - number of events in file
- TFile *file=TFile::Open("rfio:galice.root");
+ TFile *file=TFile::Open("galice.root");
if (!file->IsOpen()) {cerr<<"Can't open galice.root !\n"; exit(4);}
// delete gAlice; gAlice = 0;
if (!(gAlice=(AliRun*)file->Get("gAlice"))) {
cerr<<"gAlice have not been found on galice.root !\n";
exit(5);
}
-
+
+ TFile *fdigit = TFile::Open("digits.root");
+ file->cd();
AliTPC *TPC=(AliTPC*)gAlice->GetDetector("TPC");
Int_t ver = TPC->IsVersion();
cerr<<"TPC version "<<ver<<" has been found !\n";
AliTPCParam *digp=(AliTPCParam*)file->Get("75x40_100x60");
+ if(digp){
+ cerr<<"2 pad-lenght geom hits with 3 pad-length geom digits...\n";
+ delete digp;
+ digp = new AliTPCParamSR();
+ }
+ else
+ {
+ digp =(AliTPCParam*)gDirectory->Get("75x40_100x60_150x60");
+ }
+
if (!digp) { cerr<<"TPC parameters have not been found !\n"; exit(6); }
TPC->SetParam(digp);
break;
case 2:
{
- sprintf(treeName,"TreeD_75x40_100x60_%d",event);
- TD=(TTree*)gDirectory->Get(treeName);
+ sprintf(treeName,"TreeD_75x40_100x60_150x60_%d",event);
+ TD=(TTree*)fdigit->Get(treeName); // To be revised according to
+ // NewIO schema M.Kowalski
TD->GetBranch("Segment")->SetAddress(&digits);
count = new Int_t[np];
Int_t i;
Int_t sec,row;
digp->AdjustSectorRow(digits->GetID(),sec,row);
digits->First();
+ digits->ExpandTrackBuffer();
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);
+ // Short_t dig = digits->GetDigit(it,ip);
+ Short_t dig = digits->CurrentDigit();
+
+ Int_t idx0=digits->GetTrackIDFast(it,ip,0)-2;
+ Int_t idx1=digits->GetTrackIDFast(it,ip,1)-2;
+ Int_t idx2=digits->GetTrackIDFast(it,ip,2)-2;
if (idx0>=0 && dig>=zero && idx0<np ) count[idx0]+=1; //background events - higher track ID's
if (idx1>=0 && dig>=zero && idx1<np ) count[idx1]+=1;
if (idx2>=0 && dig>=zero && idx2<np ) count[idx2]+=1;
}
/** select tracks which are "good" enough **/
+ //printf("\t %d \n",np);
for (Int_t i=0; i<np; i++) {
if ((good[i]&0x5000) != 0x5000)
if ((good[i]&0x2800) != 0x2800) continue;
//if (!(pp->TestBit(kPrimaryCharged))) continue; //only one decay is allowed
}
- //if(!(p->TestBit(kPrimaryCharged)))continue; // only primaries
+ if(!(p->TestBit(kPrimaryCharged)))continue; // only primaries
+ // printf("1");
gt[nt].fEventN=event;
gt[nt].lab=i;
gt[nt].code=p->GetPdgCode();
gt[nt].px=0.; gt[nt].py=0.; gt[nt].pz=0.;
- gt[nt].x=0.; gt[nt].z=0.; gt[nt].z=0.;
+ gt[nt].x=0.; gt[nt].y=0.; gt[nt].z=0.;
nt++;
if (nt==max) {cerr<<"Too many good tracks !\n"; break;}
cerr<<np-i<<" \r";
Int_t j, lab=phit->Track();
for (j=0; j<nt; j++) {if (gt[j].fEventN==event && gt[j].lab==lab) break;}
if (j==nt) continue;
-
+ //printf("1-");
// (px,py,pz) - in global coordinate system, (x,y,z) - in local !
gt[j].px=px; gt[j].py=py; gt[j].pz=pz;
Float_t cs,sn; digp->AdjustCosSin(phit->fSector,cs,sn);
}
cerr<<np-i<<" \r";
}
+ //printf("\n%d\n",nt);
cout<<endl;
delete[] good;
} // /// loop on events
for(Int_t iprim = 0; iprim<nparticles; iprim++){ //loop on tracks
part = gAlice->Particle(iprim);
+ Double_t ptot=TMath::Sqrt(part->Px()*part->Px()+part->Py()*part->Py()+part->Pz()*part->Pz());
+ if (ptot<0.01) continue;
char *xxx=strstr(part->GetName(),"XXX");
if(xxx)continue;
if(part->T()>timecut)continue;
- Double_t ptot=TMath::Sqrt(part->Px()*part->Px()+part->Py()*part->Py()+part->Pz()*part->Pz());
+ //Double_t ptot=TMath::Sqrt(part->Px()*part->Px()+part->Py()*part->Py()+part->Pz()*part->Pz());
if(ptot==(TMath::Abs(part->Pz())))continue; // no beam particles
Bool_t prmch = kTRUE; // candidate primary track
return nprch1;
}
-
-
-
-
-
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff