#include <TParticle.h>
#include <TTree.h>
#include <TClonesArray.h>
-#include <TNtupleD.h>
#include <AliRun.h>
#include <AliSimDigits.h>
#include <AliTPCClustersRow.h>
#include <AliTPCParam.h>
#include <AliComplexCluster.h>
+#include <AliStack.h>
-#if GCCVERSION == 3
+#if __GNUC__ == 3
#include <fstream>
#include <iosfwd>
#else
#include "AliL3TrackArray.h"
#include "AliL3Evaluate.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
fMinSlice=0;
fMaxSlice=0;
fGoodTracks=0;
+ fNtupleRes=0;
+ fDigitsTree=0;
+ fPtRes=0;
+ fNGoodTracksPt=0;
+ fNFoundTracksPt=0;
+ fNFakeTracksPt=0;
+ fTrackEffPt=0;
+ fFakeTrackEffPt=0;
+ fNGoodTracksEta=0;
+ fNFoundTracksEta=0;
+ fNFakeTracksEta=0;
+ fTrackEffEta=0;
+ fFakeTrackEffEta=0;
+ fMcIndex=0;
+ fMcId=0;
+ fNtuppel=0;
+ fStandardComparison=kTRUE;
}
-AliL3Evaluate::AliL3Evaluate(Char_t *datapath,Int_t min_clusters,Int_t minhits,Double_t minpt,Int_t *slice)
+AliL3Evaluate::AliL3Evaluate(Char_t *datapath,Int_t min_clusters,Int_t minhits,Double_t minpt,Double_t maxpt,Int_t *slice)
{
if(slice)
fMinPointsOnTrack = min_clusters;
fMinHitsFromParticle = minhits;
fMinGoodPt = minpt;
+ fMaxGoodPt = maxpt;
memset(fClusters,0,36*6*sizeof(AliL3SpacePointData*));
fTracks=0;
fGoodTracks=0;
+ fNtupleRes=0;
+ fDigitsTree=0;
+ fPtRes=0;
+ fNGoodTracksPt=0;
+ fNFoundTracksPt=0;
+ fNFakeTracksPt=0;
+ fTrackEffPt=0;
+ fFakeTrackEffPt=0;
+ fNGoodTracksEta=0;
+ fNFoundTracksEta=0;
+ fNFakeTracksEta=0;
+ fTrackEffEta=0;
+ fFakeTrackEffEta=0;
+ fMcIndex=0;
+ fMcId=0;
+ fNtuppel=0;
+ fStandardComparison=kTRUE;
}
void AliL3Evaluate::LoadData(Int_t event,Bool_t sp)
break;
}
}
-
+
sprintf(fname,"%s/tracks_%d.raw",fPath,event);
AliL3FileHandler *tfile = new AliL3FileHandler();
if(!tfile->SetBinaryInput(fname)){
tfile->Binary2TrackArray(fTracks);
tfile->CloseBinaryInput();
delete tfile;
+ fTracks->QSort();
}
if(fMcIndex) delete [] fMcIndex;
if(fMcId) delete [] fMcId;
if(fNtuppel) delete fNtuppel;
+ if(fNtupleRes) delete fNtupleRes;
+ for(Int_t s=0; s<=35; s++)
+ for(Int_t p=0; p<6; p++)
+ if(fClusters[s][p]) delete fClusters[s][p];
}
+void AliL3Evaluate::AssignPIDs()
+{
+ fTracks->QSort();
+ LOG(AliL3Log::kDebug,"AliL3Evaluate::AssignPIDs","Track Loop")
+ <<"Assigning pid to the found tracks...."<<ENDLOG;
+ for(Int_t i=0; i<fTracks->GetNTracks(); i++)
+ {
+ AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i);
+ if(!track) continue;
+ if(track->GetNumberOfPoints() < fMinPointsOnTrack)
+ track->SetPID(0);
+ else {
+ Float_t pid = GetTrackPID(track);
+ track->SetPID(pid);
+ }
+ }
+}
void AliL3Evaluate::AssignIDs()
{
//cout<<"Found "<<fGoodFound<<" good tracks "<<endl;
}
+Float_t AliL3Evaluate::GetTrackPID(AliL3Track *track)
+{
+ track->CalculateHelix();
+ // Track dEdx
+ Int_t nc=track->GetNHits();
+ UInt_t *hits = track->GetHitNumbers();
+ Float_t sampleDEdx[159];
+ for (Int_t iHit = 0; iHit < nc; iHit++) {
+ UInt_t hitID = hits[iHit];
+ Int_t iSector = (hitID>>25) & 0x7f;
+ Int_t patch = (hitID>>22) & 0x7;
+ UInt_t position = hitID&0x3fffff;
+ AliL3SpacePointData *points = fClusters[iSector][patch];
+ if(!points) continue;
+ if(position>=fNcl[iSector][patch])
+ {
+ LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
+ <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
+ continue;
+ }
+ UChar_t padrow = points[position].fPadRow;
+ Float_t pWidth = AliL3Transform::GetPadPitchWidthLow();
+ if (padrow>63)
+ pWidth = AliL3Transform::GetPadPitchWidthUp();
+ Float_t corr=1.; if (padrow>63) corr=0.67;
+ sampleDEdx[iHit] = points[position].fCharge/pWidth*corr;
+ Double_t crossingangle = track->GetCrossingAngle(padrow,iSector);
+ Double_t s = sin(crossingangle);
+ Double_t t = track->GetTgl();
+ sampleDEdx[iHit] *= sqrt((1-s*s)/(1+t*t));
+ }
+
+ /* Cook dEdx */
+ Int_t i;
+ Int_t swap;//stupid sorting
+ do {
+ swap=0;
+ for (i=0; i<nc-1; i++) {
+ if (sampleDEdx[i]<=sampleDEdx[i+1]) continue;
+ Float_t tmp=sampleDEdx[i];
+ sampleDEdx[i]=sampleDEdx[i+1]; sampleDEdx[i+1]=tmp;
+ swap++;
+ }
+ } while (swap);
+
+ Double_t low=0.05; Double_t up=0.7;
+ Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
+ Float_t trackDEdx=0;
+ for (i=nl; i<=nu; i++) trackDEdx += sampleDEdx[i];
+ trackDEdx /= (nu-nl+1);
+
+ // cout<<" PID: "<<nc<<" "<<nl<<" "<<nu<<" "<<trackDEdx<<" "<<track->GetPt()<<endl;
+ return trackDEdx;
+}
struct S {Int_t lab; Int_t max;};
Int_t AliL3Evaluate::GetMCTrackLabel(AliL3Track *track){
fGoodGen=0;
fGoodTracks = new GoodTrack[MaxTracks];
- while (in>>fGoodTracks[fGoodGen].label>>fGoodTracks[fGoodGen].code>>
- fGoodTracks[fGoodGen].px>>fGoodTracks[fGoodGen].py>>fGoodTracks[fGoodGen].pz>>
- fGoodTracks[fGoodGen].x>>fGoodTracks[fGoodGen].y >>fGoodTracks[fGoodGen].z>>fGoodTracks[fGoodGen].nhits>>fGoodTracks[fGoodGen].sector)
-
- {
- fGoodGen++;
- if (fGoodGen==MaxTracks)
- {
- cerr<<"AliL3Evaluate::GetGoodParticles : Too many good tracks !\n";
- break;
- }
- }
+ if(fStandardComparison){
+ while (in>>fGoodTracks[fGoodGen].label>>fGoodTracks[fGoodGen].code>>
+ fGoodTracks[fGoodGen].px>>fGoodTracks[fGoodGen].py>>fGoodTracks[fGoodGen].pz>>
+ fGoodTracks[fGoodGen].x>>fGoodTracks[fGoodGen].y>>fGoodTracks[fGoodGen].z)
+ {
+ fGoodTracks[fGoodGen].nhits=-1;
+ fGoodTracks[fGoodGen].sector=-1;
+ fGoodGen++;
+ if (fGoodGen==MaxTracks)
+ {
+ cerr<<"AliL3Evaluate::GetGoodParticles : Too many good tracks !\n";
+ break;
+ }
+ }
+ } else {
+ while (in>>fGoodTracks[fGoodGen].label>>fGoodTracks[fGoodGen].code>>
+ fGoodTracks[fGoodGen].px>>fGoodTracks[fGoodGen].py>>fGoodTracks[fGoodGen].pz>>
+ fGoodTracks[fGoodGen].x>>fGoodTracks[fGoodGen].y >>fGoodTracks[fGoodGen].z>>fGoodTracks[fGoodGen].nhits>>fGoodTracks[fGoodGen].sector)
+ {
+ fGoodGen++;
+ if (fGoodGen==MaxTracks)
+ {
+ cerr<<"AliL3Evaluate::GetGoodParticles : Too many good tracks !\n";
+ break;
+ }
+ }
+ }
}
+//has to be modified for fakes.
void AliL3Evaluate::FillEffHistos()
{
for(Int_t i=0; i<fGoodGen; i++)
{
//cout<<"Checking particle "<<i<<endl;
- if(fGoodTracks[i].nhits < fMinHitsFromParticle) continue;
+ if(!fStandardComparison)
+ if(fGoodTracks[i].nhits < fMinHitsFromParticle) continue;
Float_t ptg = TMath::Sqrt(fGoodTracks[i].px*fGoodTracks[i].px + fGoodTracks[i].py*fGoodTracks[i].py);
- if(ptg < fMinGoodPt) continue;
+ if(ptg < fMinGoodPt || ptg > fMaxGoodPt) continue;
Float_t pzg=fGoodTracks[i].pz;
Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi();
+
+ //If we are only considering tracks on one side of the TPC:
+ if(fMaxSlice <= 17)
+ if(dipangle < 0)
+ continue;
fNGoodTracksPt->Fill(ptg);
fNGoodTracksEta->Fill(dipangle);
}
//if(!found)
- // cout<<"Track "<<fGoodTracks[i].label<<" was not found"<<endl;
+ //cout<<"Track "<<fGoodTracks[i].label<<" was not found"<<endl;
}
}
cout<<endl<<"Note: Doing NAIVE evaluation "<<endl;
for(Int_t i=0; i<fGoodGen; i++)
{
+ if(!fStandardComparison)
+ if(fGoodTracks[i].nhits < fMinHitsFromParticle) continue;
Double_t ptg=TMath::Sqrt(fGoodTracks[i].px*fGoodTracks[i].px + fGoodTracks[i].py*fGoodTracks[i].py);
+ if(ptg < fMinGoodPt || ptg > fMaxGoodPt) continue;
Double_t pzg=fGoodTracks[i].pz;
Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi();
//printf("filling particle with pt %f and dipangle %f\n",ptg,dipangle);
if(!track) continue;
Int_t nHits = track->GetNumberOfPoints();
if(nHits < fMinPointsOnTrack) break;
- if(track->GetPt()<fMinGoodPt) continue;
+ if(track->GetPt()<fMinGoodPt || track->GetPt() > fMaxGoodPt) continue;
if(fabs(track->GetPseudoRapidity())>0.9) continue;
fNFoundTracksPt->Fill(track->GetPt()); fNFoundTracksEta->Fill(track->GetPseudoRapidity());
}
-TNtupleD *AliL3Evaluate::CalculateResiduals(Char_t *datapath)
+TNtuple *AliL3Evaluate::GetNtuple()
+{
+ if(!fNtupleRes)
+ {
+ fNtupleRes = new TNtuple("ntuppel","Residuals","residual_trans:residual_long:zHit:pt:dipangle:beta:padrow:nHits");
+ fNtupleRes->SetDirectory(0);//Bug in older version of root.
+ }
+ return fNtupleRes;
+}
+
+void AliL3Evaluate::CalculateResiduals()
{
- TNtupleD *ntuppel=new TNtupleD("ntuppel","Residuals","residual_trans:residual_long:zHit:pt:dipangle:beta:padrow:nHits");
- ntuppel->SetDirectory(0);
- if(fTracks)
- delete fTracks;
-
- AliL3FileHandler *tfile = new AliL3FileHandler();
- char fname[256];
- sprintf(fname,"%s/tracks_0.raw",datapath);
- if(!tfile->SetBinaryInput(fname)){
- LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
- <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
- return 0;
- }
- fTracks = new AliL3TrackArray();
- tfile->Binary2TrackArray(fTracks);
- tfile->CloseBinaryInput();
- delete tfile;
+ TNtuple *ntuppel = GetNtuple();
for(Int_t i=0; i<fTracks->GetNTracks(); i++)
{
AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i);
if(!track) continue;
- if(track->GetNHits() < fMinPointsOnTrack) continue;
+ if(track->GetNHits() < fMinPointsOnTrack) break;
track->CalculateHelix();
UInt_t *hitnum = track->GetHitNumbers();
Int_t slice = (id>>25) & 0x7f;
Int_t patch = (id>>22) & 0x7;
UInt_t pos = id&0x3fffff;
-
- //if(slice!=1) continue;
+
+ //if(slice<18) continue;
AliL3SpacePointData *points = fClusters[slice][patch];
xyz[1] = points[pos].fY;
xyz[2] = points[pos].fZ;
padrow = points[pos].fPadRow;
- AliL3Transform::Global2Local(xyz,slice,kTRUE);
+ //AliL3Transform::Global2Local(xyz,slice,kTRUE);
+ AliL3Transform::Global2LocHLT(xyz,slice);
Float_t angle = 0;
AliL3Transform::Local2GlobalAngle(&angle,slice);
- track->CalculateReferencePoint(angle,AliL3Transform::Row2X(padrow));
+ if(!track->CalculateReferencePoint(angle,AliL3Transform::Row2X(padrow)))
+ {
+ LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Crossing point")
+ <<"Track does not crossing padrow "<<padrow<<" in slice "<<slice<<ENDLOG;
+ continue;
+ }
+
Float_t xyz_cross[3] = {track->GetPointX(),track->GetPointY(),track->GetPointZ()};
- AliL3Transform::Global2Local(xyz_cross,slice,kTRUE);
+ //AliL3Transform::Global2Local(xyz_cross,slice,kTRUE);
+ AliL3Transform::Global2LocHLT(xyz_cross,slice);
Double_t beta = track->GetCrossingAngle(padrow,slice);
Double_t yres = xyz_cross[1] - xyz[1];
Double_t zres = xyz_cross[2] - xyz[2];
-
Double_t dipangle = atan(track->GetTgl());
ntuppel->Fill(yres,zres,xyz_cross[2],track->GetPt(),dipangle,beta,padrow,track->GetNumberOfPoints());
}
}
- if(fTracks)
- delete fTracks;
-
- return ntuppel;
}
enum tagprimary {kPrimaryCharged = 0x4000};
return;
}
+ AliStack *astack=gAlice->Stack();
+
AliTPCClustersArray *arr=0;
for(Int_t event=0; event<nevent; event++)
{
for(Int_t m=0; m<num_of_offline; m++)
{
AliComplexCluster *cluster = (AliComplexCluster *)clusters->UncheckedAt(m);
+#ifdef use_newio
+ Int_t mcId = cluster->GetTrack(0);
+#else
Int_t mcId = cluster->fTracks[0];
-
+#endif
if(mcId <0) continue;
-
+
+#ifdef use_newio
+ if(cluster->GetY() < 1 || cluster->GetY() > AliL3Transform::GetNPads(padrow) - 2 ||
+ cluster->GetX() < 1 || cluster->GetX() > AliL3Transform::GetNTimeBins() - 2)
+ continue;
+#else
if(cluster->fY < 1 || cluster->fY > AliL3Transform::GetNPads(padrow) - 2 ||
cluster->fX < 1 || cluster->fX > AliL3Transform::GetNTimeBins() - 2)
continue;
-
+#endif
Float_t xyz_ex[3];
+#ifdef use_newio
+ AliL3Transform::Raw2Local(xyz_ex,cursec,currow,cluster->GetY(),cluster->GetX());
+#else
AliL3Transform::Raw2Local(xyz_ex,cursec,currow,cluster->fY,cluster->fX);
-
+#endif
//In function AliTPC::Hits2ExactClusters the time offset is not included,
//so we have to substract it again here.
if(slice<18)
if(param->GetPadRowRadii(cursec,currow)<230./250.*fabs(xyz_ex[2]))
continue;
- TParticle *part = gAlice->Particle(mcId);
+ TParticle *part = astack->Particle(mcId);
crosscount++;
if(part->Pt() < fMinGoodPt) continue;
sprintf(filename,"%s/alirunfile.root",path);
TFile *rfile = TFile::Open(filename);
gAlice = (AliRun*)rfile->Get("gAlice");
+
+ AliStack *astack=gAlice->Stack();
AliTPCParam *param = (AliTPCParam*)rfile->Get(AliL3Transform::GetParamName());
LoadData(event,sp);
rfile->cd();
gAlice->GetEvent(event);
- Int_t np = gAlice->GetNtrack();
+ Int_t np = astack->GetNtrack();
cout<<"Processing event "<<event<<" with "<<np<<" particles "<<endl;
dfile->cd();
sprintf(filename,"TreeD_75x40_100x60_150x60_%d",event);
} while (digits->Next());
for (Int_t j=0; j<np; j++)
{
- TParticle *part = gAlice->Particle(j);
+ TParticle *part = astack->Particle(j);
if(part->Pt() < fMinGoodPt) continue;
if(part->GetFirstMother() >= 0) continue;
if (count[j]>1) //at least two digits at this padrow
#endif
}
-Float_t AliL3Evaluate::GetCrossingAngle(TParticle *part,Int_t slice,Int_t padrow,Float_t *xyz)
+Float_t AliL3Evaluate::GetCrossingAngle(TParticle *part,Int_t slice,Int_t /*padrow*/,Float_t *xyz)
{
//Calculate the padrow crossing angle of the particle
Double_t decacut = 3.;
Double_t timecut = 0.;
Int_t nprch1=0;
- TParticle * part = gAlice->Particle(0);
+ AliStack *astack=gAlice->Stack();
+ TParticle * part = astack->Particle(0);
Double_t xori = part->Vx();
Double_t yori = part->Vy();
Double_t zori = part->Vz();
for(Int_t iprim = 0; iprim<nparticles; iprim++){ //loop on tracks
- part = gAlice->Particle(iprim);
- char *xxx=strstr(part->GetName(),"XXX");
+ part = astack->Particle(iprim);
+ const char * xxx=strstr(part->GetName(),"XXX");
if(xxx)continue;
TParticlePDG *ppdg = part->GetPDG();
}
if(ndau>0){
for(Int_t j=fidau;j<=lasdau;j++){
- TParticle *dau=gAlice->Particle(j);
+ TParticle *dau=astack->Particle(j);
Double_t distd=TMath::Sqrt((dau->Vx()-xori)*(dau->Vx()-xori)+(dau->Vy()-yori)*(dau->Vy()-yori)+(dau->Vz()-zori)*(dau->Vz()-zori));
if(distd<decacut)prmch=kFALSE; // eliminate if the decay is near the vertex
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff