ITS/UPGRADE/residuals.C

   1
   2
   3
   4 Int_t nLayers =0;
   5 Int_t GetLayer(AliTrackReference *ref, TObject *obj);
   6
   7 void residuals(){
   8
   9   gSystem->Load("libITSUpgradeBase.so");
  10   gSystem->Load("libITSUpgradeRec.so");
  11   gSystem->Load("libITSUpgradeSim.so");
  12
  13   gROOT->SetStyle("Plain");
  14   gStyle->SetOptStat(1111111);
  15   Int_t nbins=100;
  16   Int_t xmin=0;
  17   Int_t xmax=50000;//00*1e-09;
  18
  19   AliITSsegmentationUpgrade *seg = new AliITSsegmentationUpgrade();
  20
  21   if(!seg){
  22     printf("no segmentation info available... Exiting");
  23     return;
  24   }
  25
  26   nLayers = seg->GetNLayers();
  27
  28   Int_t nbins=100;
  29   Int_t xmin = -10;
  30   Int_t xmax= 10;
  31   Int_t zmin = -3000;
  32   Int_t zmax= 3000;
  33
  34   TH1D **hDiffx, **hDiffy, **hDiffz, **hDiffRphi, **hDiffR;
  35
  36   hDiffx    = new TH1D[nLayers];
  37   hDiffy    = new TH1D[nLayers];
  38   hDiffz    = new TH1D[nLayers];
  39   hDiffRphi = new TH1D[nLayers];
  40   hDiffR    = new TH1D[nLayers];
  41
  42   for(Int_t i=0; i<nLayers; i++){
  43     hDiffx[i] = new TH1D(Form("hDiffX%i",i),Form("#Delta X Glob (Clusters - track ref  [Layer %i])",i),100,-40,40);
  44     hDiffy[i] = new TH1D(Form("hDiffY%i",i),Form("#Delta Y Glob (Clusters - track ref  [Layer %i])",i),100,-40,40);
  45     if(i<4) hDiffz[i] = new TH1D(Form("hDiffZ%i",i),Form("#Delta Z Glob (Clusters - track ref  [Layer %i])",i),50,-300,300);
  46     else hDiffz[i] = new TH1D(Form("hDiffZ%i",i),Form("#Delta Z Glob (Clusters - track ref  [Layer %i])",i),100,-2000,2000);
  47     hDiffRphi[i] = new TH1D(Form("hDiffRphi%i",i),Form("#Delta R #phi [Layer %i]",i),200,-100,100);
  48     hDiffRphi[i]->SetXTitle("#mum");
  49     hDiffR[i]= new TH1D(Form("hR%i",i),Form("#Delta R (cluster-track ref) [Layer %i]",i),nbins,xmin/20.,xmax/20.);
  50     hDiffR[i]->SetXTitle("#mum");
  51   }
  52
  53   gAlice=NULL;
  54   AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
  55   runLoader->LoadgAlice();
  56
  57   gAlice = runLoader->GetAliRun();
  58
  59   runLoader->LoadHeader();
  60   runLoader->LoadKinematics();
  61   runLoader->LoadTrackRefs();
  62   runLoader->LoadRecPoints();
  63
  64   AliLoader *dl = runLoader->GetDetectorLoader("ITS");
  65
  66   //Trackf
  67   TTree *trackRefTree = 0x0;
  68   TClonesArray *trackRef = new TClonesArray("AliTrackReference",1000);
  69   //RECPOINT
  70   TTree *clusTree = 0x0;
  71   TClonesArray statITSCluster("AliITSRecPointU");
  72
  73
  74   AliITSsegmentationUpgrade *segmentation2=segmentation2 = new AliITSsegmentationUpgrade();
  75   // Loop over events
  76   for(Int_t i=0; i<runLoader->GetNumberOfEvents(); i++){
  77     runLoader->GetEvent(i);
  78     AliStack *stack = runLoader->Stack();
  79     trackRefTree=runLoader->TreeTR();
  80     TBranch *br = trackRefTree->GetBranch("TrackReferences");
  81     if(!br) {
  82       printf("no TR branch available , exiting \n");
  83       return;
  84     }
  85     br->SetAddress(&trackRef);
  86
  87     printf("event : %i \n",i);
  88     // init the clusters tree
  89     clusTree=dl->TreeR();
  90     TClonesArray *ITSCluster = &statITSCluster;
  91     TBranch* itsClusterBranch=clusTree->GetBranch("ITSRecPoints");
  92     if (!itsClusterBranch) {
  93       printf("can't get the branch with the ITS clusters ! \n");
  94       return;
  95     }
  96
  97     itsClusterBranch->SetAddress(&ITSCluster);
  98
  99     // init the trackRef tree
 100     trackRefTree=runLoader->TreeTR();
 101     trackRefTree->SetBranchAddress("TrackReferences",&trackRef);
 102
 103     //for(Int_t il=0; il<nLayers; il++){
 104     if (!clusTree->GetEvent(0))    continue;
 105     Int_t nCluster = ITSCluster->GetEntriesFast();
 106     for(Int_t in=0; in<nCluster; in++){
 107       AliITSRecPointU *recp = (AliITSRecPointU*)ITSCluster->UncheckedAt(in);
 108       Int_t il = recp->GetLayer();
 109       Double_t xr,yr,zr = 0.;
 110       Double_t xz[2];
 111       xz[0]= recp->GetDetLocalX(); //gets fXloc
 112       xz[1]= recp->GetDetLocalZ();   //gets fZloc
 113       segmentation2->DetToGlobal(il,recp->GetModule(),xz[0], xz[1],xr,yr,zr);
 114       for(Int_t iLabel =0; iLabel<3; iLabel++){
 115         if(recp->GetLabel(iLabel)<0) continue;
 116         trackRefTree->GetEntry(stack->TreeKEntry(recp->GetLabel(iLabel)));
 117         Int_t nref=trackRef->GetEntriesFast();
 118         for(Int_t iref =0; iref<nref; iref++){
 119           Double_t x,y,z=0.;
 120           Int_t labelTR=-999;
 121           AliTrackReference *trR = (AliTrackReference*)trackRef->At(iref);
 122           if(!trR) continue;
 123           if(trR->DetectorId()!=AliTrackReference::kITS) continue;
 124           if(!stack->IsPhysicalPrimary(recp->GetLabel(iLabel))) continue;
 125           Int_t lay = GetLayer(trR, seg);
 126           if(TMath::Abs(lay-il)>0.5) continue;
 127           x=trR->X();
 128           y=trR->Y();
 129           z=trR->Z();
 130           Double_t xx, zz;
 131           Int_t mod;
 132           segmentation2->GlobalToDet(il,x,y,x,xx,zz,mod);
 133           //printf("cluster module %i  - module %i \n",recp->GetModule(),mod);
 134           Double_t rTr = TMath::Sqrt(x*x+y*y);
 135
 136           Double_t difx, dify, difz=0.;
 137           Double_t dify=0.;
 138           Double_t difz=0.;
 139           Double_t deltaR=0.;
 140           difx = xr - x;
 141           dify = yr - y;
 142           difz = zr - z;
 143           Double_t phiGlob = TMath::ATan2(yr,xr);
 144           if(yr<0) phiGlob+=TMath::TwoPi();
 145           Double_t phiTr = TMath::ATan2(y,x);
 146           if(y<0) phiTr+=TMath::TwoPi();
 147           Double_t deltaPhi =  phiGlob-phiTr;
 148
 149           Double_t deltaRphi = (phiGlob*TMath::Sqrt(xr*xr+yr*yr)-phiTr*rTr)*1e+04;
 150           if(TMath::Abs(deltaRphi)/1e+04 > 3*(10000.*seg->GetCellSizeX(lay))) {
 151             // printf("   Layer %i   type %i  dRphi %f  segm %f  \n",lay,recp->GetType(),deltaRphi/1e+04,seg->GetCellSizeX(lay)*10000.);
 152             // printf("   Layer %i   type %i   Track (x,y,z) = (%f,%f,%f) Phi %f   Rphi %f         Cluster (x,y,z) = (%f,%f,%f)  Phi %f Rphi %f \n",lay,recp->GetType(),x,y,z,TMath::RadToDeg()*phiTr,phiTr*TMath::Sqrt(x*x+y*y),xr,yr,zr,TMath::RadToDeg()*phiGlob,phiGlob*TMath::Sqrt(xr*xr+yr*yr));
 153           }
 154           deltaR = (TMath::Sqrt(xr*xr+yr*yr)-TMath::Sqrt(x*x+y*y))*1e+04;
 155
 156           hDiffx[il]->Fill(difx*1e04);
 157           hDiffy[il]->Fill(dify*1e04);
 158           hDiffz[il]->Fill(difz*1e04);
 159           hDiffRphi[il]->Fill(deltaRphi);
 160           hDiffR[il]->Fill(deltaR);
 161
 162         }//loop entries fast
 163       }
 164     }//loop clusters
 165     //}//loop layer
 166   }//loop entries runloader
 167
 168   TCanvas *c3x = new TCanvas("cGx","Delta in X global",200,10,900,900);
 169   c3x->Divide(3,nLayers/3);
 170   TCanvas *c3y = new TCanvas("cGy","Delta in Y global",200,10,900,900);
 171   c3y->Divide(3,nLayers/3);
 172   TCanvas *c3z = new TCanvas("cGz","Delta in Z global",200,10,900,900);
 173   c3z->Divide(3,nLayers/3);
 174   TCanvas *c3rphi = new TCanvas("cRphi","Delta R-Phi",200,10,900,900);
 175   c3rphi->Divide(3,nLayers/3);
 176   TCanvas *c3r = new TCanvas("cR","Delta R-Phi",200,10,900,900);
 177   c3r->Divide(3,nLayers/3);
 178
 179
 180   for(Int_t iL=0; iL< nLayers; iL++){
 181     c3x->cd(iL+1);
 182     hDiffx[iL]->Draw();
 183     c3y->cd(iL+1);
 184     hDiffy[iL]->Draw();
 185     c3z->cd(iL+1);
 186     hDiffz[iL]->Draw();
 187     c3rphi->cd(iL+1);
 188     hDiffRphi[iL]->Draw();
 189     c3r->cd(iL+1);
 190     hDiffR[iL]->Draw();
 191   }
 192
 193 }// main macro read
 194
 195
 196 Int_t GetLayer(AliTrackReference *ref, TObject *obj){
 197
 198   AliITSsegmentationUpgrade *seg = (AliITSsegmentationUpgrade*)obj;
 199   Int_t ilayer =-1;
 200   for(Int_t ila=0; ila<6; ila++){
 201     if(ila==(nLayers-1)){
 202       if(ref->R()>seg->GetRadius(ila)+seg->GetThickness(ila)) ilayer=(nLayers-1);
 203     }
 204     else {
 205       if(ref->R()>seg->GetRadius(ila)+seg->GetThickness(ila) && ref->R()<seg->GetRadius(ila+1)) ilayer=ila;
 206     }
 207   }
 208   return ilayer;
 209 }
 210