[u/mrichter/AliRoot.git] / ITS / UPGRADE / residuals.C



Int_t nLayers =0;
Int_t GetLayer(AliTrackReference *ref, TObject *obj);

void residuals(){

  gSystem->Load("libITSUpgradeBase.so");
  gSystem->Load("libITSUpgradeRec.so");
  gSystem->Load("libITSUpgradeSim.so");

  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(1111111);
  Int_t nbins=100;
  Int_t xmin=0;
  Int_t xmax=50000;//00*1e-09;

  AliITSsegmentationUpgrade *seg = new AliITSsegmentationUpgrade();
  
  if(!seg){
    printf("no segmentation info available... Exiting");
    return;
  }

  nLayers = seg->GetNLayers();

  Int_t nbins=100;
  Int_t xmin = -10;
  Int_t xmax= 10;
  Int_t zmin = -3000;
  Int_t zmax= 3000;

  TH1D **hDiffx, **hDiffy, **hDiffz, **hDiffRphi, **hDiffR;
  
  hDiffx    = new TH1D[nLayers];
  hDiffy    = new TH1D[nLayers];
  hDiffz    = new TH1D[nLayers];
  hDiffRphi = new TH1D[nLayers];
  hDiffR    = new TH1D[nLayers];
    
  for(Int_t i=0; i<nLayers; i++){
    hDiffx[i] = new TH1D(Form("hDiffX%i",i),Form("#Delta X Glob (Clusters - track ref  [Layer %i])",i),100,-40,40);
    hDiffy[i] = new TH1D(Form("hDiffY%i",i),Form("#Delta Y Glob (Clusters - track ref  [Layer %i])",i),100,-40,40);
    if(i<4) hDiffz[i] = new TH1D(Form("hDiffZ%i",i),Form("#Delta Z Glob (Clusters - track ref  [Layer %i])",i),50,-300,300);
    else hDiffz[i] = new TH1D(Form("hDiffZ%i",i),Form("#Delta Z Glob (Clusters - track ref  [Layer %i])",i),100,-2000,2000);
    hDiffRphi[i] = new TH1D(Form("hDiffRphi%i",i),Form("#Delta R #phi [Layer %i]",i),200,-100,100);
    hDiffRphi[i]->SetXTitle("#mum");
    hDiffR[i]= new TH1D(Form("hR%i",i),Form("#Delta R (cluster-track ref) [Layer %i]",i),nbins,xmin/20.,xmax/20.);
    hDiffR[i]->SetXTitle("#mum");
  }

  gAlice=NULL;
  AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
  runLoader->LoadgAlice();

  gAlice = runLoader->GetAliRun();

  runLoader->LoadHeader();
  runLoader->LoadKinematics();
  runLoader->LoadTrackRefs();
  runLoader->LoadRecPoints();

  AliLoader *dl = runLoader->GetDetectorLoader("ITS");

  //Trackf
  TTree *trackRefTree = 0x0; 
  TClonesArray *trackRef = new TClonesArray("AliTrackReference",1000);
  //RECPOINT
  TTree *clusTree = 0x0;
  TClonesArray statITSCluster("AliITSRecPointU");


  AliITSsegmentationUpgrade *segmentation2=segmentation2 = new AliITSsegmentationUpgrade();
  // Loop over events
  for(Int_t i=0; i<runLoader->GetNumberOfEvents(); i++){
    runLoader->GetEvent(i); 
    AliStack *stack = runLoader->Stack(); 
    trackRefTree=runLoader->TreeTR();
    TBranch *br = trackRefTree->GetBranch("TrackReferences");
    if(!br) {
      printf("no TR branch available , exiting \n");
      return;
    }
    br->SetAddress(&trackRef);

    printf("event : %i \n",i);
    // init the clusters tree 
    clusTree=dl->TreeR();
    TClonesArray *ITSCluster = &statITSCluster;
    TBranch* itsClusterBranch=clusTree->GetBranch("ITSRecPoints");
    if (!itsClusterBranch) {
      printf("can't get the branch with the ITS clusters ! \n");
      return;
    }

    itsClusterBranch->SetAddress(&ITSCluster);

    // init the trackRef tree 
    trackRefTree=runLoader->TreeTR();
    trackRefTree->SetBranchAddress("TrackReferences",&trackRef);

    //for(Int_t il=0; il<nLayers; il++){
    if (!clusTree->GetEvent(0))    continue;
    Int_t nCluster = ITSCluster->GetEntriesFast();
    for(Int_t in=0; in<nCluster; in++){
      AliITSRecPointU *recp = (AliITSRecPointU*)ITSCluster->UncheckedAt(in);
      Int_t il = recp->GetLayer();  
      Double_t xr,yr,zr = 0.;
      Double_t xz[2];
      xz[0]= recp->GetDetLocalX(); //gets fXloc
      xz[1]= recp->GetDetLocalZ();   //gets fZloc
      segmentation2->DetToGlobal(il,recp->GetModule(),xz[0], xz[1],xr,yr,zr);
      for(Int_t iLabel =0; iLabel<12; iLabel++){
	if(recp->GetTrackID(iLabel)<0) continue; 
	trackRefTree->GetEntry(stack->TreeKEntry(recp->GetTrackID(iLabel)));  	
	Int_t nref=trackRef->GetEntriesFast();
	for(Int_t iref =0; iref<nref; iref++){
	  Double_t x,y,z=0.;
	  Int_t labelTR=-999;
	  AliTrackReference *trR = (AliTrackReference*)trackRef->At(iref);
	  if(!trR) continue;
	  if(trR->DetectorId()!=AliTrackReference::kITS) continue;
	  if(!stack->IsPhysicalPrimary(recp->GetTrackID(iLabel))) continue;   
	  Int_t lay = GetLayer(trR, seg);    
	  if(TMath::Abs(lay-il)>0.5) continue;          
	  x=trR->X();
	  y=trR->Y();
	  z=trR->Z();
	  Double_t xx, zz;
	  Int_t mod;
	  segmentation2->GlobalToDet(il,x,y,x,xx,zz,mod);
	  //printf("cluster module %i  - module %i \n",recp->GetModule(),mod);
	  Double_t rTr = TMath::Sqrt(x*x+y*y); 

	  Double_t difx, dify, difz=0.;
	  Double_t dify=0.;
	  Double_t difz=0.;
	  Double_t deltaR=0.;
	  difx = xr - x;
	  dify = yr - y;
	  difz = zr - z;
	  Double_t phiGlob = TMath::ATan2(yr,xr);
	  if(yr<0) phiGlob+=TMath::TwoPi();
	  Double_t phiTr = TMath::ATan2(y,x);
	  if(y<0) phiTr+=TMath::TwoPi();
	  Double_t deltaPhi =  phiGlob-phiTr;
	      
	  Double_t deltaRphi = (phiGlob*TMath::Sqrt(xr*xr+yr*yr)-phiTr*rTr)*1e+04;
	  if(TMath::Abs(deltaRphi)/1e+04 > 3*(10000.*seg->GetCellSizeX(lay))) {
	    // printf("   Layer %i   type %i  dRphi %f  segm %f  \n",lay,recp->GetType(),deltaRphi/1e+04,seg->GetCellSizeX(lay)*10000.);
	    // 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));
	  }              
	  deltaR = (TMath::Sqrt(xr*xr+yr*yr)-TMath::Sqrt(x*x+y*y))*1e+04;

	  hDiffx[il]->Fill(difx*1e04);
	  hDiffy[il]->Fill(dify*1e04);
	  hDiffz[il]->Fill(difz*1e04);
	  hDiffRphi[il]->Fill(deltaRphi);
	  hDiffR[il]->Fill(deltaR);
           
	}//loop entries fast 
      } 
    }//loop clusters
    //}//loop layer
  }//loop entries runloader
  
  TCanvas *c3x = new TCanvas("cGx","Delta in X global",200,10,900,900);
  c3x->Divide(3,nLayers/3);
  TCanvas *c3y = new TCanvas("cGy","Delta in Y global",200,10,900,900);
  c3y->Divide(3,nLayers/3);
  TCanvas *c3z = new TCanvas("cGz","Delta in Z global",200,10,900,900);
  c3z->Divide(3,nLayers/3);
  TCanvas *c3rphi = new TCanvas("cRphi","Delta R-Phi",200,10,900,900);
  c3rphi->Divide(3,nLayers/3);
  TCanvas *c3r = new TCanvas("cR","Delta R-Phi",200,10,900,900);
  c3r->Divide(3,nLayers/3);
  

  for(Int_t iL=0; iL< nLayers; iL++){
    c3x->cd(iL+1);
    hDiffx[iL]->Draw();
    c3y->cd(iL+1);
    hDiffy[iL]->Draw();
    c3z->cd(iL+1);
    hDiffz[iL]->Draw();
    c3rphi->cd(iL+1);
    hDiffRphi[iL]->Draw();
    c3r->cd(iL+1);
    hDiffR[iL]->Draw();
  }
 
}// main macro read


Int_t GetLayer(AliTrackReference *ref, TObject *obj){
 
  AliITSsegmentationUpgrade *seg = (AliITSsegmentationUpgrade*)obj;
  Int_t ilayer =-1;
  for(Int_t ila=0; ila<6; ila++){
    if(ila==(nLayers-1)){
      if(ref->R()>seg->GetRadius(ila)+seg->GetThickness(ila)) ilayer=(nLayers-1);
    }
    else {
      if(ref->R()>seg->GetRadius(ila)+seg->GetThickness(ila) && ref->R()<seg->GetRadius(ila+1)) ilayer=ila;
    }
  }
  return ilayer;
}
Commit	Line	Data
556b741c	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);
9b615954	114	for(Int_t iLabel =0; iLabel<12; iLabel++){
	115	if(recp->GetTrackID(iLabel)<0) continue;
	116	trackRefTree->GetEntry(stack->TreeKEntry(recp->GetTrackID(iLabel)));
556b741c	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;
9b615954	124	if(!stack->IsPhysicalPrimary(recp->GetTrackID(iLabel))) continue;
556b741c	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(xx+yy);
	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 = (phiGlobTMath::Sqrt(xrxr+yryr)-phiTrrTr)*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,phiTrTMath::Sqrt(xx+yy),xr,yr,zr,TMath::RadToDeg()phiGlob,phiGlobTMath::Sqrt(xrxr+yryr));
	153	}
	154	deltaR = (TMath::Sqrt(xrxr+yryr)-TMath::Sqrt(xx+yy))*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