-#ifndef __CINT__
- #include <iostream.h>
-
- #include "AliTRDtracker.h"
- #include "AliTRDcluster.h"
- #include "AliTRDhit.h"
- #include "AliTRDv1.h"
- #include "AliTRDgeometry.h"
- #include "AliTRDparameter.h"
-
- #include "alles.h"
- #include "TFile.h"
- #include "TStopwatch.h"
-
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <Riostream.h>
+
+#include "AliTRDtracker.h"
+#include "AliTRDclusterMI.h"
+#include "AliTRDhit.h"
+#include "AliTRDv1.h"
+#include "AliTRDgeometry.h"
+#include "AliTRDgeometryDetail.h"
+#include "AliTRDparameter.h"
+#include "AliTRDclusterCorrection.h"
+#include "alles.h"
+#include "TFile.h"
+#include "TStopwatch.h"
+#include "Rtypes.h"
+#include "TTree.h"
+
+#include "AliRunLoader.h"
+#include "AliStack.h"
+#include "TF1.h"
+#include "AliTrackReference.h"
#endif
+#include "AliTRDclusterErrors.h"
-void AliTRDclusterErrors() {
-
- Int_t No_of_tracks_to_analyze = 40;
-
- const Int_t tbpp = 15;
- const Int_t nPlanes = 6;
- const Int_t ntb = tbpp * nPlanes;
-
- TH1F *hy = new TH1F("delta R*phi","Cluster displacement in R*phi",200,-1.,1.);
- TH1F *hyp = new TH1F("delta R*phi pos","delta R*phi, positive",200,-1.,1.);
- TH1F *hym = new TH1F("delta R*phi neg","delta R*phi, negative",200,-1.,1.);
-
- TH1F *hyn = new TH1F("Norm., d(R*phi)","Norm. cluster displacement in R*phi",400,-8.,8.);
- TH1F *hz = new TH1F("delta Z","Cluster displacement in Z",300,-10.,50.);
- TH2F *hy2 = new TH2F("Amp vs delta R*phi","Amplitude versus delta R*phi",200,-5.,5.,200,0.,600.);
- TH2F *herr = new TH2F("sigmaY vs delta R*phi","sigmaY vs delta R*phi",200,-1,1,200,0.,0.1);
- TH2F *hy3 = new TH2F("Position within pad vs delta R*phi","Position within pad vs delta R*phi",200,-1.,1.,200,-0.5,1.5);
- TH2F *hy4 = new TH2F("local tb vs delta R*phi","local tb vs delta R*phi",200,-1.,1.,20,-2.5,17.5);
-
- hy->SetXTitle("Displacement, cm");
- hyn->SetXTitle("Displacement, SigmaY");
- hy2->SetXTitle("Displacement, cm");
- hy2->SetYTitle("Amplitude");
- hz->SetXTitle("Displacement, cm");
- hy3->SetXTitle("Displacement, cm");
- hy3->SetYTitle("Position, cm");
- hy4->SetXTitle("Displacement, cm");
- hy4->SetYTitle("local time bin");
-
- /*
- // Dynamically link some shared libs
- if (gClassTable->GetID("AliRun") < 0) {
- gROOT->LoadMacro("loadlibs.C");
- loadlibs();
- cout << "Loaded shared libraries" << endl;
- }
- */
-
- // Load clusters
- Char_t *clusterfile = "AliTRDclusters.root";
- Int_t nEvent = 0;
- TObjArray carray(2000);
- TObjArray *ClustersArray = &carray;
- TFile *geofile =TFile::Open("AliTRDclusters.root");
+ClassImp(AliTRDCI)
+ClassImp(AliTRDExactPoint)
+ClassImp(AliTRDClusterErrAnal)
+ClassImp(AliTRDClusterErrDraw)
- AliTRDparameter *par = (AliTRDparameter*) geofile->Get("TRDparameter");
- AliTRDgeometry *fGeo = (AliTRDgeometry*) geofile->Get("TRDgeometry");
- AliTRDtracker *Tracker = new AliTRDtracker(geofile);
- Tracker->SetEventNumber(nEvent);
- Tracker->ReadClusters(ClustersArray,clusterfile);
- Int_t nClusters = carray.GetEntriesFast();
- printf("Total number of clusters %d \n", nClusters);
-
- Char_t *alifile = "galice.root";
- // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits
- TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile);
- if (!gafl) {
- cout << "Open the ALIROOT-file " << alifile << endl;
- gafl = new TFile(alifile);
+AliTRDClusterErrAnal::AliTRDClusterErrAnal(Char_t *chloader )
+{
+ //
+ //SET Input loaders
+ if (gAlice){
+ delete gAlice->GetRunLoader();
+ delete gAlice;
+ gAlice = 0x0;
+ }
+ fRunLoader = AliRunLoader::Open(chloader);
+ if (fRunLoader == 0x0){
+ cerr<<"Can not open session"<<endl;
+ return;
}
- else {
- cout << alifile << " is already open" << endl;
+ fTRDLoader = fRunLoader->GetLoader("TRDLoader");
+ if (fTRDLoader == 0x0){
+ cerr<<"Can not get TRD Loader"<<endl;
+ return ;
+ }
+ if (fRunLoader->LoadgAlice()){
+ cerr<<"Error occured while l"<<endl;
+ return;
}
+ fRunLoader->CdGAFile();
+ fTracker = new AliTRDtracker(gFile);
+ fParam = (AliTRDparameter*) gFile->Get("TRDparameter");
+ fGeometry = new AliTRDgeometryDetail();
+ fTRD = (AliTRDv1*) gAlice->GetDetector("TRD");
+ //
+ AliTRDCI * clinfo = new AliTRDCI();
+ fFileA = new TFile("trdclusteranal.root","recreate");
+ fFileA->cd();
+ fTreeA = new TTree("trdcl","trdcl");
+ fTreeA->Branch("trdcl","AliTRDCI",&clinfo);
+ delete clinfo;
+}
- // Get AliRun object from file or create it if not on file
- gAlice = (AliRun*) gafl->Get("gAlice");
- if (gAlice)
- cout << "AliRun object found on file" << endl;
- else
- gAlice = new AliRun("gAlice","Alice test program");
-
- AliTRDv1 *fTRD = (AliTRDv1*) gAlice->GetDetector("TRD");
-
- // Import the Trees for the event nEvent in the file
- Int_t nparticles = gAlice->GetEvent(nEvent);
-
- TObjArray *particles=gAlice->Particles();
-
- TTree *hitTree = gAlice->TreeH();
- Int_t nBytes = 0;
-
- // Get the number of entries in the hit tree
- // (Number of primary particles creating a hit somewhere)
- Int_t nTrack = (Int_t) hitTree->GetEntries();
- cout << " Found " << nTrack << " primary particles with hits" << endl;
- No_of_tracks_to_analyze = TMath::Min(No_of_tracks_to_analyze, nTrack);
-
-
- // Loop through particles and fill histoes
-
- Float_t hitY[ntb];
- Float_t hitZ[ntb];
- Float_t hitO[ntb];
- Float_t clusterY[ntb];
- Float_t clusterZ[ntb];
- Float_t clusterQ[ntb];
- Float_t clusterSigmaY[ntb];
- Float_t pos[3];
- Float_t rot[3];
- Float_t global[3];
- Int_t det = 0;
- Int_t track_index[3];
-
- printf("\n");
-
- for (Int_t ii=0; ii<No_of_tracks_to_analyze; ii++) {
- printf("track %d out of %d \n", ii+1 , No_of_tracks_to_analyze);
-
- for(Int_t plane = 0; plane < nPlanes; plane++) {
- for(Int_t ltb = 14; ltb > -1; ltb--) {
-
- if(ii >= nTrack) continue;
+void AliTRDClusterErrAnal::SetIO(Int_t event)
+{
+ //
+ //set input output for given event
+ fRunLoader->SetEventNumber(event);
+ fRunLoader->LoadHeader();
+ fRunLoader->LoadKinematics();
+ fRunLoader->LoadTrackRefs();
+ fTRDLoader->LoadHits();
+ fTRDLoader->LoadRecPoints("read");
+ //
+ fStack = fRunLoader->Stack();
+ fHitTree = fTRDLoader->TreeH();
+ fClusterTree = fTRDLoader->TreeR();
+ fReferenceTree = fRunLoader->TreeTR();
+ //
+}
- TParticle *p = gAlice->Particle(ii);
- if (p->GetFirstMother()>=0) continue;
- TParticlePDG *pdg = p->GetPDG();
- Float_t charge=pdg->Charge();
- if(TMath::Abs(charge) < 0.5) continue;
+void AliTRDClusterErrAnal::LoadClusters()
+{
+ //
+ //
+ // Load clusters
+ TObjArray *ClusterArray = new TObjArray(400);
+ TObjArray carray(2000);
+ TBranch *branch=fClusterTree->GetBranch("TRDcluster");
+ if (!branch) {
+ Error("ReadClusters","Can't get the branch !");
+ return;
+ }
+ branch->SetAddress(&ClusterArray);
+ Int_t nentries = (Int_t)fClusterTree->GetEntries();
+ for (Int_t i=0;i<nentries;i++){
+ fClusterTree->GetEvent(i);
+ Int_t nCluster = ClusterArray->GetEntriesFast();
+ for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
+ AliTRDcluster* c = (AliTRDcluster*)ClusterArray->UncheckedAt(iCluster);
+ carray.AddLast(c);
+ ClusterArray->RemoveAt(iCluster);
+ }
+ }
+ Int_t nClusters = carray.GetEntriesFast();
+ printf("Total number of clusters %d \n", nClusters);
+ //
+ //
+ //SORT clusters
+ //
+ Int_t all=0;
+ for (Int_t i=0;i<nClusters;i++){
+ AliTRDcluster *cl = (AliTRDcluster *) carray.UncheckedAt(i);
+ Int_t plane = fGeometry->GetPlane(cl->GetDetector());
+ if (plane>=12) continue;
+ Int_t time = cl->GetLocalTimeBin();
+ if (time>=100) continue;
+ Int_t sector = fGeometry->GetSector(cl->GetDetector());
+ if (sector>=18){
+ printf("problem1\n");
+ continue;
+ }
+ fClusters[plane][time][sector].AddLast(cl);
+ all++;
+ }
+}
- Int_t gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
- Double_t x = Tracker->GetX(0,plane,ltb);
+void AliTRDClusterErrAnal::SortReferences()
+{
+ //
+ //
+ //
+ printf("Sorting references\n");
+ fReferences = new TObjArray;
+ Int_t ntracks = fStack->GetNtrack();
+ fReferences->Expand(ntracks);
+ Int_t nentries = (Int_t)fReferenceTree->GetEntries();
+ TClonesArray * arr = new TClonesArray("AliTrackReference");
+ TBranch * br = fReferenceTree->GetBranch("TRD");
+ br->SetAddress(&arr);
+ //
+ TClonesArray *labarr=0;
+ Int_t nreferences=0;
+ Int_t nreftracks=0;
+ for (Int_t iprim=0;iprim<nentries;iprim++){
+ if (br->GetEntry(iprim)){
+ for (Int_t iref=0;iref<arr->GetEntriesFast();iref++){
+ AliTrackReference *ref =(AliTrackReference*)arr->At(iref);
+ if (!ref) continue;
+ Int_t lab = ref->GetTrack();
+ if ( (lab<0) || (lab>ntracks)) continue;
+ //
+ if (fReferences->At(lab)==0) {
+ labarr = new TClonesArray("AliTrackReference");
+ fReferences->AddAt(labarr,lab);
+ nreftracks++;
+ }
+ TClonesArray &larr = *labarr;
+ new(larr[larr.GetEntriesFast()]) AliTrackReference(*ref);
+ nreferences++;
+ }
+ }
+ }
+ printf("Total number of references = \t%d\n", nreferences);
+ printf("Total number of tracks with references = \t%d\n", nreftracks);
+ printf("End - Sorting references\n");
+
+}
- // loop through clusters
+AliTrackReference * AliTRDClusterErrAnal::FindNearestReference(Int_t lab, Float_t pos[3], Float_t dmax)
+{
+ //
+ //
+ //
+ if (fReferences->At(lab)==0) return 0;
+ AliTrackReference *nearest=0;
+ TClonesArray * arr = (TClonesArray *)fReferences->At(lab);
+ for (Int_t iref =0;iref<arr->GetEntriesFast();iref++){
+ AliTrackReference * ref = ( AliTrackReference *)arr->UncheckedAt(iref);
+ if (!ref) continue;
+ Float_t delta = (pos[0]-ref->X())*(pos[0]-ref->X());
+ delta += (pos[1]-ref->Y())*(pos[1]-ref->Y());
+ delta += (pos[2]-ref->Z())*(pos[2]-ref->Z());
+ delta = TMath::Sqrt(delta);
+ if (delta<dmax){
+ dmax=delta;
+ nearest = ref;
+ }
+ }
+ return nearest;
+}
- Bool_t cluster_found = kFALSE;
- Int_t nw = 0;
+void AliTRDClusterErrAnal::MakeExactPoints(Int_t trackmax)
+{
+ //
+ //make exact points:-)
+ //
+ //
+
+ fExactPoints.Delete();
+ fExactPoints.Expand(fStack->GetNtrack());
+ //
+ Double_t fSum=0;
+ Double_t fSumQ =0;
+ Double_t fSumX=0;
+ Double_t fSumX2=0;
+ Double_t fSumXY=0;
+ Double_t fSumXZ=0;
+ Double_t fSumY=0;
+ Double_t fSumZ=0;
+ //
+ Int_t entries = Int_t(fHitTree->GetEntries());
+ printf("Number of primary entries\t%d\n",entries);
+ entries = TMath::Min(trackmax,entries);
+ Int_t nallpoints = 0;
+
+ Int_t nalltracks =0;
+ Int_t pointspertrack =0;
+
+ for (Int_t entry=0;entry<entries; entry++){
+ gAlice->ResetHits();
+ fHitTree->GetEvent(entry);
+ Int_t lastlabel = -1;
+ Int_t lastdetector = -1;
+ Int_t lasttimebin = -1;
+ Float_t lastpos[3];
+ //
+ for(AliTRDhit *hit = (AliTRDhit *) fTRD->FirstHit(-1); hit;
+ hit = (AliTRDhit *) fTRD->NextHit()) {
+ //
+ Int_t label = hit->Track();
+ TParticle * particle = fStack->Particle(label);
+ if (!particle) continue;
+ if (particle->Pt()<0.05) continue;
+ Int_t detector = hit->GetDetector();
+ Int_t plane = fGeometry->GetPlane(detector);
+ //
+ //
+ if (hit->GetCharge()==0) continue;
+ Float_t pos[3] = {hit->X(),hit->Y(),hit->Z()};
+ Int_t indexes[3];
+ fGeometry->Global2Detector(pos,indexes,fParam);
+ //
+ Float_t rot[3];
+ fGeometry->Rotate(detector,pos,rot);
+ //rot[0] *=-1;
+ // rot[1] *=-1;
+ //
+ //
+ Float_t time0 = fParam->GetTime0(plane);
+ Int_t timebin = Int_t(TMath::Nint(((time0 - rot[0])/fParam->GetTimeBinSize())+ fParam->GetTimeBefore())+0.1);
+ if (timebin<0) continue;
+ //
+ //
+ if (label!=lastlabel || detector != lastdetector || lasttimebin !=timebin){
+ //
+ if (label!=lastlabel){
+ fExactPoints.AddAt(new TClonesArray("AliTRDExactPoint",0),label);
+ //printf("new particle\t%d\n",label);
+ nalltracks++;
+ // printf("particle\t%d- hits\t%d\n",lastlabel, pointspertrack);
+ pointspertrack=0;
+ }
+
+ if ( (fSum>1) && lasttimebin>=0 && lasttimebin<fParam->GetTimeMax() ){
+ //if we have enough info for given layer time bin - store it
+ AliTrackReference * ref = FindNearestReference(lastlabel,lastpos,4.);
+ Float_t rotmom[3];
+ Float_t rotpos[3];
+ Float_t refangle=0;
+ if (ref){
+ Float_t mom[3] = {ref->Px(),ref->Py(),ref->Pz()};
+ Float_t refpos[3] = {ref->X(),ref->Y(),ref->Z()};
+ fGeometry->Rotate(detector,mom,rotmom);
+ fGeometry->Rotate(detector,refpos,rotpos);
+ refangle = rotmom[1]/rotmom[0];
+
+ }
+
+ Double_t ay,by,az,bz;
+ Double_t det = fSum*fSumX2-fSumX*fSumX;
+ if (TMath::Abs(det)> 0.000000000000001) {
+ by = (fSum*fSumXY-fSumX*fSumY)/det;
+ ay = (fSumX2*fSumY-fSumX*fSumXY)/det;
+
+ }else{
+ ay =fSumXY/fSumX;
+ by =0;
+ }
+ if (TMath::Abs(det)> 0.000000000000001) {
+ bz = (fSum*fSumXZ-fSumX*fSumZ)/det;
+ az = (fSumX2*fSumZ-fSumX*fSumXZ)/det;
+ }else{
+ az =fSumXZ/fSumX;
+ bz =0;
+ }
+ //
+ Float_t lastplane = fGeometry->GetPlane(lastdetector);
+ Float_t time0 = fParam->GetTime0(lastplane);
+ Float_t xcenter0 = time0 - (lasttimebin - fParam->GetTimeBefore()+0.5)*fParam->GetTimeBinSize();
+ Float_t xcenter = fTracker->GetX(0,lastplane,lasttimebin);
+ if (TMath::Abs(xcenter-xcenter0)>0.001){
+ printf("problem");
+ }
+
+ Float_t ty = ay + by * xcenter;
+ Float_t tz = az + bz * xcenter;
+ //
+
+ TClonesArray * arr = (TClonesArray *) fExactPoints.At(label);
+ TClonesArray & larr= *arr;
+ Int_t arrent = arr->GetEntriesFast();
+ AliTRDExactPoint * point = new (larr[arrent]) AliTRDExactPoint;
+ nallpoints++;
+
+ if (ref){
+ point->SetReference(ref);
+ point->fTRefAngleY = rotmom[1]/rotmom[0];
+ }
+ point->fTX = xcenter;
+ point->fTY = ty;
+ point->fTZ = tz;
+ point->fTAY = by;
+ point->fTAZ = bz;
+ //
+ point->fGx = lastpos[0];
+ point->fGy = lastpos[1];
+ point->fGz = lastpos[2];
+
+ //
+ point->fDetector = lastdetector;
+ point->fLocalTimeBin = lasttimebin;
+ point->fPlane = fGeometry->GetPlane(lastdetector);
+ point->fSector = fGeometry->GetSector(lastdetector);
+ point->fPlaneMI = indexes[0];
+ //
+ point->fTPrim = fSum;
+ point->fTQ = fSumQ;
+ //
+ }
+ lastdetector = detector;
+ lastlabel = label;
+ lasttimebin = timebin;
+ fSum=fSumQ=fSumX=fSumX2=fSumXY=fSumXZ=fSumY=fSumZ=0.;
+ }
+ //
+ lastpos[0] = hit->X();
+ lastpos[1] = hit->Y();
+ lastpos[2] = hit->Z();
+ fSum++;
+ fSumQ +=hit->GetCharge();
+ fSumX +=rot[0];
+ fSumX2 +=rot[0]*rot[0];
+ fSumXY +=rot[0]*rot[1];
+ fSumXZ +=rot[0]*rot[2];
+ fSumY +=rot[1];
+ fSumZ +=rot[2];
+ pointspertrack++;
+ }
+ }
+ //
+ printf("Found %d exact points\n",nallpoints);
+}
- for (Int_t i = 0; i < nClusters; i++) {
- AliTRDcluster *cl = (AliTRDcluster *) carray.UncheckedAt(i);
- nw = cl->GetDetector();
- nw = fGeo->GetPlane(nw);
- if(nw != plane) continue;
- for(Int_t j=0; j<3; j++) track_index[j] = cl->GetLabel(j);
- if((track_index[0] != ii) &&
- (track_index[1] != ii) &&
- (track_index[2] != ii)) continue;
- nw = cl->GetLocalTimeBin();
- if(nw != ltb) continue;
+Int_t AliTRDClusterErrAnal::Analyze(Int_t trackmax) {
- clusterY[gtb] = cl->GetY();
- clusterZ[gtb] = cl->GetZ();
- clusterQ[gtb] = cl->GetQ();
+ //
+ // comparison works with both cluster types MI and old also
+ //dummy cluster to be fill if not cluster info
+ AliTRDclusterMI clmi;
+ TClass * classmi = clmi.IsA();
+ //
+ //SetOutput
+ AliTRDCI * clinfo = new AliTRDCI();
+ TBranch * clbr = fTreeA->GetBranch("trdcl");
+ clbr->SetAddress(&clinfo);
- clusterSigmaY[gtb] = TMath::Sqrt(cl->GetSigmaY2());
- cluster_found = kTRUE;
- break;
+ SetIO(0);
+ SortReferences();
+ MakeExactPoints(trackmax);
+ LoadClusters();
+ //
+ trackmax = fStack->GetNtrack();
+ //
+ // Get the number of entries in the hit tree
+ // (Number of primary particles creating a hit somewhere)
+ Int_t nTrack = (Int_t)fExactPoints.GetEntries();
+ printf("Found %d charged in TRD in first %d particles", nTrack, trackmax);
+ //
+
+ for (Int_t itrack = 0; itrack<trackmax; itrack++){
+ TClonesArray *arrpoints = (TClonesArray*)fExactPoints.At(itrack);
+
+ if (!arrpoints) continue;
+ //printf("new particle\t%d\n",itrack);
+ TParticle * particle = fStack->Particle(itrack);
+ if (!particle) continue;
+ //printf("founded in kine tree \t%d\n",itrack);
+ Int_t npoints = arrpoints->GetEntriesFast();
+ if (npoints<10) continue;
+ //printf("have enough points \t%d\t%d\n",itrack,npoints);
+
+ for (Int_t ipoint=0;ipoint<npoints;ipoint++){
+ AliTRDExactPoint * point = (AliTRDExactPoint*)arrpoints->UncheckedAt(ipoint);
+ if (!point) continue;
+ //
+ Int_t sec = fGeometry->GetSector(point->fDetector);
+ if (sec>18){
+ printf("problem2\n");
+ }
+ TObjArray & cllocal = fClusters[point->fPlane][point->fLocalTimeBin][sec];
+ Int_t nclusters = cllocal.GetEntriesFast();
+ Float_t maxdist = 10;
+ AliTRDcluster * nearestcluster =0;
+ //find nearest cluster to hit with given label
+ for (Int_t icluster =0; icluster<nclusters; icluster++){
+ AliTRDcluster * cluster = (AliTRDcluster*)cllocal.UncheckedAt(icluster);
+ if (!cluster) continue;
+ if ( (cluster->GetLabel(0)!= itrack) && (cluster->GetLabel(1)!= itrack)&&(cluster->GetLabel(2)!= itrack))
+ continue;
+ Float_t dist = TMath::Abs(cluster->GetY()-point->fTY);
+ if (TMath::Abs(cluster->GetZ()-point->fTZ)>5.5) continue;
+ if (dist<maxdist){
+ maxdist = dist;
+ nearestcluster = cluster;
+ }
+ }
+ //
+ clinfo->fEp = *point;
+ clinfo->fP = *particle;
+ if (!nearestcluster) {
+ clinfo->fStatus=1;
+ clinfo->fCl = clmi;
+ }
+ else{
+ clinfo->fStatus=0;
+ if (nearestcluster->IsA()==classmi){
+ clinfo->fCl =*((AliTRDclusterMI*)nearestcluster);
}
+ else{
+ clinfo->fCl = *nearestcluster;
+ }
+ //
+ Float_t dz = clinfo->fCl.GetZ()-point->fTZ;
+ Double_t h01 = sin(TMath::Pi() / 180.0 * fParam->GetTiltingAngle());
+ clinfo->fTDistZ = dz;
+ clinfo->fDYtilt = h01*dz*((point->fPlane%2)*2.-1.);
+ //
+ clinfo->fNTracks =1;
+ if (nearestcluster->GetLabel(1)>=0) clinfo->fNTracks++;
+ if (nearestcluster->GetLabel(2)>=0) clinfo->fNTracks++;
+ clinfo->Update();
+ }
+ //
+ fTreeA->Fill();
+ }
+ }
+
+
+ fFileA->cd();
+ fTreeA->Write();
+ fFileA->Close();
+ return 0;
+}
- if(!cluster_found) continue;
+AliTRDclusterCorrection* AliTRDClusterErrDraw::MakeCorrection(TTree * tree, Float_t offset)
+{
+ //
+ //
+ // make corrections
+ AliTRDclusterCorrection * cor = new AliTRDclusterCorrection;
+ cor->SetOffsetAngle(offset);
+ for (Int_t iplane=0;iplane<6;iplane++)
+ for (Int_t itime=0;itime<15;itime++)
+ for (Int_t iangle=0; iangle<20;iangle++){
+ Float_t angle = cor->GetAngle(iangle);
+ TH1F delta("delta","delta",30,-0.3,0.3);
+ char selection[100]="fStatus==0&&fNTracks<2";
+ char selectionall[1000];
+ sprintf(selectionall,"%s&&abs(fEp.fTAY-%f)<0.2&&fEp.fPlane==%d&&fCl.fTimeBin==%d&&fCl.fQ>20",
+ selection,angle,iplane,itime);
+ printf("\n%s",selectionall);
+ tree->Draw("fEp.fTY-fCl.fY+fDYtilt>>delta",selectionall);
+ gPad->Update();
+ printf("\nplane\t%d\ttime%d\tangle%f",iplane,itime,angle);
+ printf("\tentries%f\tmean\t%f\t%f",delta.GetEntries(),delta.GetMean(),delta.GetRMS());
+ cor->SetCorrection(iplane,itime,angle,delta.GetMean(),delta.GetRMS());
+ }
+ TFile * f = new TFile("TRDcorrection.root","new");
+ if (!f) f = new TFile("TRDcorrection.root","recreate");
+ f->cd();
+ cor->Write("TRDcorrection");
+ f->Close();
+ return cor;
+}
- gAlice->ResetHits();
+TH1F * AliTRDClusterErrDraw::ResDyVsAmp(TTree* tree, const char* selection, Float_t t0, Float_t ampmin, Float_t ampmax)
+{
+ //
+ //
+ TH2F hisdy("resy","resy",10,ampmin,ampmax,30,-0.3,0.3);
+ char expression[1000];
+ sprintf(expression,"fEp.fTY-fCl.fY+fDYtilt+%.4f*fEp.fTAY:fCl.fQ>>resy",t0);
+ char selectionall[1000];
+ sprintf(selectionall,"fStatus==0&&%s",selection);
+ printf("%s\n",expression);
+ printf("%s\n",selectionall);
+ tree->Draw(expression,selectionall);
+ return CreateResHisto(&hisdy);
+}
- // nBytes += hitTree->GetEvent(nPrimaries - ii - 1);
- nBytes += hitTree->GetEvent(nTrack - ii - 1);
- // Loop through the TRD hits
- Bool_t found_hit = kFALSE;
+TH1F * AliTRDClusterErrDraw::ResDyVsRelPos(TTree* tree, const char* selection, Float_t t0, Float_t min, Float_t max)
+{
+ //
+ //
+ min *=128;
+ max *=128;
+ TH2F hisdy("resy","resy",10,min,max,30,-0.3,0.3);
+ char expression[1000];
+ sprintf(expression,"fEp.fTY-fCl.fY+fDYtilt+%.4f*fEp.fTAY:fCl.fRelPos>>resy",t0);
+ char selectionall[1000];
+ sprintf(selectionall,"fStatus==0&&%s",selection);
+ printf("%s\n",expression);
+ printf("%s\n",selectionall);
+ tree->Draw(expression,selectionall);
+ return CreateResHisto(&hisdy);
+}
- for(AliTRDhit *hit = (AliTRDhit *) fTRD->FirstHit(-1);
- hit;
- hit = (AliTRDhit *) fTRD->NextHit()) {
- nw = hit->Track();
- if(nw != ii) continue;
- det = hit->GetDetector();
- nw = fGeo->GetPlane(det);
- if(nw != plane) continue;
+TH1F * AliTRDClusterErrDraw::ResDyVsAngleY(TTree* tree, const char* selection, Float_t t0, Float_t min, Float_t max)
+{
+ //
+ //
+ TH2F hisdy("resy","resy",10,min,max,30,-0.3,0.3);
+ char expression[1000];
+ sprintf(expression,"fEp.fTY-fCl.fY+fDYtilt+%f*fEp.fTAY:fEp.fTAY>>resy",t0);
+ char selectionall[1000];
+ sprintf(selectionall,"fStatus==0&&%s",selection);
- pos[0]=hit->X();
- pos[1]=hit->Y();
- pos[2]=hit->Z();
- fGeo->Rotate(det,pos,rot);
+ tree->Draw(expression,selectionall);
+ return CreateResHisto(&hisdy);
+}
- if(TMath::Abs(rot[0]-x) > 0.01) continue;
- hitY[gtb] = rot[1];
- hitZ[gtb] = rot[2];
+void AliTRDClusterErrDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
+{
+ histo->GetXaxis()->SetTitle(xAxisTitle);
+ histo->GetYaxis()->SetTitle(yAxisTitle);
+}
- Float_t col0 = par->GetCol0(plane);
- Float_t colPadSize = par->GetColPadSize(plane);
- Float_t colH = (Int_t ((rot[1] - col0)/colPadSize)) * colPadSize;
- hitO[gtb] = (rot[1] - col0) - colH;
- found_hit = kTRUE;
- break;
- }
- if(!found_hit) continue;
+TH1F* AliTRDClusterErrDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec)
+{
+ Int_t nBins = hGen->GetNbinsX();
+ TH1F* hEff = (TH1F*) hGen->Clone("hEff");
+ hEff->SetTitle("");
+ hEff->SetStats(kFALSE);
+ hEff->SetMinimum(0.);
+ hEff->SetMaximum(110.);
+
+ for (Int_t iBin = 0; iBin <= nBins; iBin++) {
+ Double_t nGen = hGen->GetBinContent(iBin);
+ Double_t nRec = hRec->GetBinContent(iBin);
+ if (nGen > 0) {
+ Double_t eff = nRec/nGen;
+ hEff->SetBinContent(iBin, 100. * eff);
+ Double_t error = sqrt((eff*(1.-eff)+0.01) / nGen);
+ // if (error == 0) error = sqrt(0.1/nGen);
+ //
+ if (error == 0) error = 0.0001;
+ hEff->SetBinError(iBin, 100. * error);
+ } else {
+ hEff->SetBinContent(iBin, 100. * 0.5);
+ hEff->SetBinError(iBin, 100. * 0.5);
+ }
+ }
- /*
- printf("gtb: %d, x: %f, rot[0]: %f, Yhit: %f, Ycl: %f\n",
- gtb, x, rot[0], rot[1], clusterY[gtb]);
- printf("\n Zhit - Zcl = %f - %f = %f\n",
- rot[2], clusterZ[gtb], rot[2] - clusterZ[gtb]);
+ return hEff;
+}
-
- printf("found hit within dx = %f - %f \n",rot[0],x);
- printf("pos: %f, %f, %f \n",pos[0],pos[1],pos[2]);
- printf("rot: %f, %f, %f \n",rot[0],rot[1],rot[2]);
- printf("cluster: %d, %f, %f \n",gtb,clusterY[gtb],clusterZ[gtb]);
- */
+TH1F* AliTRDClusterErrDraw::CreateResHisto(TH2F* hRes2, Bool_t draw, Bool_t drawBinFits,
+ Bool_t overflowBinFits)
+{
+ TVirtualPad* currentPad = gPad;
+ TAxis* axis = hRes2->GetXaxis();
+ Int_t nBins = axis->GetNbins();
+ TH1F* hRes, *hMean;
+ if (axis->GetXbins()->GetSize()){
+ hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
+ hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
+ }
+ else{
+ hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
+ hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
- hy->Fill(hitY[gtb]-clusterY[gtb]);
- if(charge > 0) hyp->Fill(hitY[gtb]-clusterY[gtb]);
- else if(charge < 0) hym->Fill(hitY[gtb]-clusterY[gtb]);
+ }
+ hRes->SetStats(false);
+ hRes->SetOption("E");
+ hRes->SetMinimum(0.);
+ //
+ hMean->SetStats(false);
+ hMean->SetOption("E");
+
+ // create the fit function
+ //TKFitGaus* fitFunc = new TKFitGaus("resFunc");
+ TF1 * fitFunc = new TF1("G","[3]+[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
+
+ fitFunc->SetLineWidth(2);
+ fitFunc->SetFillStyle(0);
+ // create canvas for fits
+ TCanvas* canBinFits = NULL;
+ Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
+ Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
+ Int_t ny = (nPads-1) / nx + 1;
+ if (drawBinFits) {
+ canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
+ if (canBinFits) delete canBinFits;
+ canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
+ canBinFits->Divide(nx, ny);
+ }
- if((clusterQ[gtb]>10)&&(clusterSigmaY[gtb]>0))
- hyn->Fill((hitY[gtb]-clusterY[gtb])/clusterSigmaY[gtb]);
- hz->Fill(hitZ[gtb]-clusterZ[gtb]);
- hy2->Fill(hitY[gtb]-clusterY[gtb],clusterQ[gtb]);
- hy3->Fill(hitY[gtb]-clusterY[gtb],hitO[gtb]);
- hy4->Fill(hitY[gtb]-clusterY[gtb],(Float_t)(tbpp - 1 - gtb%tbpp));
- herr->Fill(hitY[gtb]-clusterY[gtb],clusterSigmaY[gtb]);
+ // loop over x bins and fit projection
+ Int_t dBin = ((overflowBinFits) ? 1 : 0);
+ for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
+ if (drawBinFits) canBinFits->cd(bin + dBin);
+ TH1D* hBin = hRes2->ProjectionY("hBin", bin, bin);
+ //
+ if (hBin->GetEntries() > 10) {
+ fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS(),0.02*hBin->GetMaximum());
+ hBin->Fit(fitFunc,"s");
+ Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
+
+ if (sigma > 0.){
+ hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
+ hMean->SetBinContent(bin, fitFunc->GetParameter(1));
+ }
+ else{
+ hRes->SetBinContent(bin, 0.);
+ hMean->SetBinContent(bin,0);
}
+ hRes->SetBinError(bin, fitFunc->GetParError(2));
+ hMean->SetBinError(bin, fitFunc->GetParError(1));
+
+ //
+ //
+
+ } else {
+ hRes->SetBinContent(bin, 0.);
+ hRes->SetBinError(bin, 0.);
+ hMean->SetBinContent(bin, 0.);
+ hMean->SetBinError(bin, 0.);
}
+
+
+ if (drawBinFits) {
+ char name[256];
+ if (bin == 0) {
+ sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
+ } else if (bin == nBins+1) {
+ sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
+ } else {
+ sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
+ axis->GetTitle(), axis->GetBinUpEdge(bin));
+ }
+ canBinFits->cd(bin + dBin);
+ hBin->SetTitle(name);
+ hBin->SetStats(kTRUE);
+ hBin->DrawCopy("E");
+ canBinFits->Update();
+ canBinFits->Modified();
+ canBinFits->Update();
+ }
+
+ delete hBin;
+ }
+
+ delete fitFunc;
+ currentPad->cd();
+ if (draw){
+ currentPad->Divide(1,2);
+ currentPad->cd(1);
+ hRes->Draw();
+ currentPad->cd(2);
+ hMean->Draw();
}
- gStyle->SetOptStat(1);
- gStyle->SetOptFit(1);
-
- TCanvas* c = new TCanvas("c", "c", 110, 110, 810, 840);
- c->SetFillColor(10);
- c->Divide(2,2);
- c->cd(1); hy->SetLineWidth(2); hy->SetFillColor(29); hy->Draw();
- c->cd(2); hz->SetLineWidth(2); hz->SetFillColor(29); hz->Draw();
- c->cd(3); hyn->Draw();
- c->cd(4); hy4->Draw();
-
- TCanvas* c1 = new TCanvas("c1", "c1", 210, 210, 910, 940);
- c1->SetFillColor(10);
- c1->Divide(2,2);
- c1->cd(1); hyp->SetLineWidth(2); hyp->SetFillColor(29); hyp->Fit("gaus");
- c1->cd(2); hym->SetLineWidth(2); hym->SetFillColor(29); hym->Fit("gaus");
- c1->cd(3); hy3->Draw();
- c1->cd(4); hy2->Draw();
-
+ return hRes;
}
-
-
-
-
-
-
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff