[u/mrichter/AliRoot.git] / ITS / AliITSAnalizeSPDHits.C

void AliITSAnalizeSPDHits(TString hfn="galice.root",Int_t mod=-1,
                     Int_t evnt=-1){
  // Macro to analize hits in the SPD to chatch posible problems

  // Dynamically link some shared libs
  if (gClassTable->GetID("AliRun") < 0) {
    gROOT->LoadMacro("loadlibs.C");
    loadlibs();
  } else {
    if(gAlice){
      delete AliRunLoader::Instance();
      delete gAlice;
      gAlice=0;
    }
  }
  gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSstandard.C");
  // Set OCDB if needed
  AliCDBManager* man = AliCDBManager::Instance();
  if (!man->IsDefaultStorageSet()) {
    printf("Setting a local default storage and run number 0\n");
    man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
    man->SetRun(0);
  }else {
    printf("Using deafult storage \n");
  }
  // retrives geometry 
  TString geof(gSystem->DirName(hfn));
  geof += "/geometry.root";
  TGeoManager::Import(geof.Data());
  if (!gGeoManager) {
    cout<<"geometry not found\n";
    return -1;
  }


  AliRunLoader *rl = AccessFile(hfn); // Set up to read in Data
  Int_t retval = rl->LoadHeader();
  if (retval){
    cerr<<"AliITSPrintHits.C : LoadHeader returned error"<<endl;
    return;
  }

  AliITSLoader* ITSloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");

  if(!ITSloader){
    cerr<<"AliITSPrintHits.C :  ITS loader not found"<<endl;
    return;
  }

  ITSloader->LoadHits("read");
  AliITS *ITS = (AliITS*)gAlice->GetDetector("ITS");
  if(!ITS || ITS==0){
    cout << "Error: no ITS found. Aborting"<<endl;
    return;
  } // end if !ITS
  //cout << ITS << endl;
  AliITSInitGeometry *initgeom = new AliITSInitGeometry(ITS->GetMajorVersion(),ITS->GetMinorVersion());
  //cout << initgeom << endl;
  AliITSgeom *geom = initgeom->CreateAliITSgeom();
  if(!geom){
    cout << "Error: not AliITSgeom object found."<<endl;
    return;
  } // end if geom
  ITSloader->SetITSgeom(geom);

  Int_t evNumber1 = 0;
  Int_t evNumber2 = AliRunLoader::GetNumberOfEvents();
  if(evnt>=0){
    evNumber1 = evnt;
    evNumber2 = evnt+1;
  } // end if evnt>=0
  Int_t mod1 = 0;
  Int_t mod2 = geom->GetIndexMax();
  if(mod>=0){
    mod1 = mod;
    mod2 = mod+1;
  } // end if mod>=0
  AliITShit *hp = 0;

  Int_t i;
  Double_t zbin[161];
  for(i=0;i<161;i++) zbin[i] = 425.0E-4;
  zbin[32] = zbin[64] = zbin[96] = zbin[128] = 625.0E-4;
  zbin[0] = -3.536;
  for(i=1;i<161;i++) zbin[i] += zbin[i-1];
  TH1I *phiDistribution1 = new TH1I("phiDistribution1",
		   "Phi Distribution for the inner most layer of the SPD",
				    5027,0.0,360.0);
  TH2I *detCoverage1 = new TH2I("detCoverage1",
	      "The SPD layer 1 cell coverage summed over the first ladder per sector",
				160,(Double_t *)zbin,256,-0.64,+0.64);
  TH2I *detCoverage1a = new TH2I("detCoverage1a",
	      "The SPD layer 1 cell coverage sum over the second ladder per sector",
				 160,(Double_t *)zbin,256,-0.64,+0.64);
  TH2I *detThicknessZ1 = new TH2I("detThicknessZ1",
                    "Hit local y distribution as a function of z",
				  160,(Double_t*)zbin,200,-100.E-4,+100.E-4);
  TH2I *detThicknessX1 = new TH2I("detThicknessX1",
                    "Hit local y distribution as a function of x",
				  256,-0.64,+0.64,200,-100.E-4,+100.E-4);
  TH1I *phiDistribution2 = new TH1I("phiDistribution2",
		   "Phi Distribution for the outer most layer of the SPD",
				    8800,0.0,360.0);
  TH2I *detCoverage2 = new TH2I("detCoverage2",
	      "The SPD layer 2 cell coverage summed over all but the forth ladder per sector",
				160,(Double_t *)zbin,256,-0.64,+0.64);
  TH2I *detCoverage2a = new TH2I("detCoverage2a",
	      "The SPD layer 2 cell coverage summed over all forth ladder per sector",
				 160,(Double_t *)zbin,256,-0.64,+0.64);
  TH2I *detThicknessZ2 = new TH2I("detThicknessZ2",
                    "Hit local y distribution as a function of z",
				  160,(Double_t*)zbin,200,-100.E-4,+100.E-4);
  TH2I *detThicknessX2 = new TH2I("detThicknessX2",
                    "Hit local y distribution as a function of x",
				  256,-0.64,+0.64,200,-100.E-4,+100.E-4);
  TH2I *xySpace = new TH2I("xySpace",
			   "ALICE global coorinate location of hits in X,y",
			   100,-8.0,+8.0,100,-8.0,+8.0);
  TH2I *zrSpace = new TH2I("zrSpace",
			   "ALICE global coorinate location of hits in z,r",
			   100,-15.0,+15.0,100,0.0,+8.0);
  Double_t xg,yg,zg,xl,yl,zl,phi;
  Int_t nmodules,size=-1;
  Int_t event,m,i,i2,hit,trk,lay,lad,det;
  for(event = evNumber1; event < evNumber2; event++){
    //cout<<"Processing event "<<event<<endl;
    rl->GetEvent(event);
    ITS->InitModules(size,nmodules);
    ITS->FillModules(event,0,-1," "," ");
    for(m=mod1;m<mod2;m++)if(geom->GetModuleType(m)==(AliITSDetector)(0)){
      i2 = (ITS->GetModule(m))->GetNhits();
      //cout <<  "Event=" << event << " module=" << m <<
      //  " Number of Hits=" << i2 <<endl;
      for(i=0;i<i2;i++){
        trk = (ITS->GetModule(m))->GetHitTrackIndex(i);
        hit = (ITS->GetModule(m))->GetHitHitIndex(i);
        hp  = (ITS->GetModule(m))->GetHit(i);
	hp->GetPositionG(xg,yg,zg);
	hp->GetPositionL(xl,yl,zl);
	geom->GetModuleId(m,lay,lad,det);
	phi = TMath::ATan2(yg,xg)*TMath::RadToDeg();
	if(phi<0.0) phi+=360.0;
	switch(lay){
	case 1:
	  phiDistribution1->Fill(phi,1.0);
	  if(lad%2==1) detCoverage1->Fill(zl,xl,1.0);
	  if(lad%2==0) detCoverage1a->Fill(zl,xl,1.0);
	  detThicknessZ1->Fill(zl,yl,1.0);
	  detThicknessX1->Fill(xl,yl,1.0);
	  break;
	case 2:
	  phiDistribution2->Fill(phi,1.0);
	  if(lad%4!=0) detCoverage2->Fill(zl,xl,1.0);
	  if(lad%4==0) detCoverage2a->Fill(zl,xl,1.0);
	  detThicknessZ2->Fill(zl,yl,1.0);
	  detThicknessX2->Fill(xl,yl,1.0);
	break;
	default:
	} // end switch
      } // end for i
    } // end for m
    ITS->ClearModules();
  } // end for event
  //
  cout << "Creating pad"<<endl;
  TVirtualPad *pad=0;
  TCanvas *c0 = new TCanvas("c0","SPD Hits Tests",262,50,700,534);
  c0->Range(0,0,1,1);
  c0->Divide(3,2);
  //
  pad = c0->cd(1);
  phiDistribution2->Draw();
  pad = c0->cd(4);
  phiDistribution1->Draw();
  TPad *pad2 = c0->cd(2);
  delete pad2;
  c0->cd(0);
  TPad *pad21 = new TPad("pad21","CoverageZ",0.34,0.875,0.65,1.0);
  pad21->Draw();
  pad21->cd();
  pad21->Range(0,0,1,1);
  pad21->Modified();
  detCoverage2->Draw();
  c0->Update();
  c0->cd(0);
  TPad *pad22 = new TPad("pad22","CoverageZa",0.34,0.75,0.65,0.875);
  pad22->Draw();
  pad22->cd();
  pad22->Range(0,0,1,1);
  pad22->Modified();
  detCoverage2a->Draw();
  c0->Update();
  c0->cd(0);
  TPad *pad23 = new TPad("pad23","ThicknessZ",0.34,0.625,0.65,0.75);
  pad23->Draw();
  pad23->cd();
  pad23->Range(0,0,1,1);
  pad23->Modified();
  detThicknessZ2->Draw();
  c0->Update();
  c0->cd(0);
  TPad *pad24 = new TPad("pad24","ThicknessX",0.34,0.50,0.65,0.625);
  pad24->Draw();
  pad24->cd();
  pad24->Range(0,0,1,1);
  pad24->Modified();
  detThicknessX2->Draw();
  c0->Update();
  //
  c0->cd(0);
  TPad *pad5 = c0->cd(5);
  delete pad5;
  c0->cd(0);
  TPad *pad51 = new TPad("pad51","CoverageZ",0.34,0.375,0.65,0.50);
  pad51->Draw();
  pad51->cd();
  pad51->Range(0,0,1,1);
  pad51->Modified();
  detCoverage1->Draw();
  c0->Update();
  c0->cd(0);
  TPad *pad52 = new TPad("pad52","CoverageZa",0.34,0.25,0.65,0.375);
  pad52->Draw();
  pad52->cd();
  pad52->Range(0,0,1,1);
  pad52->Modified();
  detCoverage1a->Draw();
  c0->Update();
  c0->cd(0);
  TPad *pad53 = new TPad("pad53","ThicknessZ",0.34,0.125,0.65,0.25);
  pad53->Draw();
  pad53->cd();
  pad53->Range(0,0,1,1);
  pad53->Modified();
  detThicknessZ1->Draw();
  c0->Update();
  c0->cd(0);
  TPad *pad54 = new TPad("pad54","ThicknessX",0.34,0.00,0.65,0.125);
  pad54->Draw();
  pad54->cd();
  pad54->Range(0,0,1,1);
  pad54->Modified();
  detThicknessX1->Draw();
  c0->Update();
  //
  TList *tvTreeList = new TList;
  TTree *tvTree = (TTree *) gROOT->FindObject("TreeH");
  tvTreeList->Add(tvTree);
  pad = c0->cd(3);
  tvTree->Draw("ITS.fY:ITS.fX","ITS.fModule<240","", 100000, 0);
  pad = c0->cd(6);
  tvTree->Draw("ITS.fX*ITS.fX+ITS.fY*ITS.fY:ITS.fZ","ITS.fModule<240","", 
		 100000, 0);

}
Commit	Line	Data
b9202d71	1	void AliITSAnalizeSPDHits(TString hfn="galice.root",Int_t mod=-1,
	2	Int_t evnt=-1){
	3	// Macro to analize hits in the SPD to chatch posible problems
	4
	5	// Dynamically link some shared libs
	6	if (gClassTable->GetID("AliRun") < 0) {
	7	gROOT->LoadMacro("loadlibs.C");
	8	loadlibs();
	9	} else {
	10	if(gAlice){
33c3c91a	11	delete AliRunLoader::Instance();
b9202d71	12	delete gAlice;
	13	gAlice=0;
	14	}
	15	}
	16	gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSstandard.C");
	17	// Set OCDB if needed
	18	AliCDBManager* man = AliCDBManager::Instance();
	19	if (!man->IsDefaultStorageSet()) {
	20	printf("Setting a local default storage and run number 0\n");
162637e4	21	man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
b9202d71	22	man->SetRun(0);
	23	}else {
	24	printf("Using deafult storage \n");
	25	}
	26	// retrives geometry
	27	TString geof(gSystem->DirName(hfn));
	28	geof += "/geometry.root";
	29	TGeoManager::Import(geof.Data());
	30	if (!gGeoManager) {
	31	cout<<"geometry not found\n";
	32	return -1;
	33	}
	34
	35
	36	AliRunLoader *rl = AccessFile(hfn); // Set up to read in Data
	37	Int_t retval = rl->LoadHeader();
	38	if (retval){
	39	cerr<<"AliITSPrintHits.C : LoadHeader returned error"<<endl;
	40	return;
	41	}
	42
	43	AliITSLoader* ITSloader = (AliITSLoader*) rl->GetLoader("ITSLoader");
	44
	45	if(!ITSloader){
	46	cerr<<"AliITSPrintHits.C : ITS loader not found"<<endl;
	47	return;
	48	}
	49
	50	ITSloader->LoadHits("read");
	51	AliITS ITS = (AliITS)gAlice->GetDetector("ITS");
	52	if(!ITS \|\| ITS==0){
	53	cout << "Error: no ITS found. Aborting"<<endl;
	54	return;
	55	} // end if !ITS
	56	//cout << ITS << endl;
	57	AliITSInitGeometry *initgeom = new AliITSInitGeometry(ITS->GetMajorVersion(),ITS->GetMinorVersion());
	58	//cout << initgeom << endl;
	59	AliITSgeom *geom = initgeom->CreateAliITSgeom();
	60	if(!geom){
	61	cout << "Error: not AliITSgeom object found."<<endl;
	62	return;
	63	} // end if geom
	64	ITSloader->SetITSgeom(geom);
	65
	66	Int_t evNumber1 = 0;
880d6abe	67	Int_t evNumber2 = AliRunLoader::GetNumberOfEvents();
b9202d71	68	if(evnt>=0){
	69	evNumber1 = evnt;
	70	evNumber2 = evnt+1;
	71	} // end if evnt>=0
	72	Int_t mod1 = 0;
	73	Int_t mod2 = geom->GetIndexMax();
	74	if(mod>=0){
	75	mod1 = mod;
	76	mod2 = mod+1;
	77	} // end if mod>=0
	78	AliITShit *hp = 0;
	79
	80	Int_t i;
	81	Double_t zbin[161];
	82	for(i=0;i<161;i++) zbin[i] = 425.0E-4;
	83	zbin[32] = zbin[64] = zbin[96] = zbin[128] = 625.0E-4;
	84	zbin[0] = -3.536;
	85	for(i=1;i<161;i++) zbin[i] += zbin[i-1];
	86	TH1I *phiDistribution1 = new TH1I("phiDistribution1",
	87	"Phi Distribution for the inner most layer of the SPD",
	88	5027,0.0,360.0);
	89	TH2I *detCoverage1 = new TH2I("detCoverage1",
	90	"The SPD layer 1 cell coverage summed over the first ladder per sector",
	91	160,(Double_t *)zbin,256,-0.64,+0.64);
	92	TH2I *detCoverage1a = new TH2I("detCoverage1a",
	93	"The SPD layer 1 cell coverage sum over the second ladder per sector",
	94	160,(Double_t *)zbin,256,-0.64,+0.64);
	95	TH2I *detThicknessZ1 = new TH2I("detThicknessZ1",
	96	"Hit local y distribution as a function of z",
	97	160,(Double_t*)zbin,200,-100.E-4,+100.E-4);
	98	TH2I *detThicknessX1 = new TH2I("detThicknessX1",
	99	"Hit local y distribution as a function of x",
	100	256,-0.64,+0.64,200,-100.E-4,+100.E-4);
	101	TH1I *phiDistribution2 = new TH1I("phiDistribution2",
	102	"Phi Distribution for the outer most layer of the SPD",
	103	8800,0.0,360.0);
	104	TH2I *detCoverage2 = new TH2I("detCoverage2",
	105	"The SPD layer 2 cell coverage summed over all but the forth ladder per sector",
	106	160,(Double_t *)zbin,256,-0.64,+0.64);
	107	TH2I *detCoverage2a = new TH2I("detCoverage2a",
	108	"The SPD layer 2 cell coverage summed over all forth ladder per sector",
	109	160,(Double_t *)zbin,256,-0.64,+0.64);
	110	TH2I *detThicknessZ2 = new TH2I("detThicknessZ2",
	111	"Hit local y distribution as a function of z",
	112	160,(Double_t*)zbin,200,-100.E-4,+100.E-4);
	113	TH2I *detThicknessX2 = new TH2I("detThicknessX2",
	114	"Hit local y distribution as a function of x",
	115	256,-0.64,+0.64,200,-100.E-4,+100.E-4);
	116	TH2I *xySpace = new TH2I("xySpace",
	117	"ALICE global coorinate location of hits in X,y",
	118	100,-8.0,+8.0,100,-8.0,+8.0);
	119	TH2I *zrSpace = new TH2I("zrSpace",
	120	"ALICE global coorinate location of hits in z,r",
	121	100,-15.0,+15.0,100,0.0,+8.0);
	122	Double_t xg,yg,zg,xl,yl,zl,phi;
	123	Int_t nmodules,size=-1;
	124	Int_t event,m,i,i2,hit,trk,lay,lad,det;
	125	for(event = evNumber1; event < evNumber2; event++){
	126	//cout<<"Processing event "<<event<<endl;
	127	rl->GetEvent(event);
	128	ITS->InitModules(size,nmodules);
	129	ITS->FillModules(event,0,-1," "," ");
	130	for(m=mod1;m<mod2;m++)if(geom->GetModuleType(m)==(AliITSDetector)(0)){
	131	i2 = (ITS->GetModule(m))->GetNhits();
132	//cout << "Event=" << event << " module=" << m <<
133	// " Number of Hits=" << i2 <<endl;
134	for(i=0;i<i2;i++){
135	trk = (ITS->GetModule(m))->GetHitTrackIndex(i);
136	hit = (ITS->GetModule(m))->GetHitHitIndex(i);
137	hp = (ITS->GetModule(m))->GetHit(i);
138	hp->GetPositionG(xg,yg,zg);
139	hp->GetPositionL(xl,yl,zl);
140	geom->GetModuleId(m,lay,lad,det);
141	phi = TMath::ATan2(yg,xg)*TMath::RadToDeg();
142	if(phi<0.0) phi+=360.0;
143	switch(lay){
144	case 1:
145	phiDistribution1->Fill(phi,1.0);
146	if(lad%2==1) detCoverage1->Fill(zl,xl,1.0);
147	if(lad%2==0) detCoverage1a->Fill(zl,xl,1.0);
148	detThicknessZ1->Fill(zl,yl,1.0);
149	detThicknessX1->Fill(xl,yl,1.0);
150	break;
151	case 2:
152	phiDistribution2->Fill(phi,1.0);
153	if(lad%4!=0) detCoverage2->Fill(zl,xl,1.0);
154	if(lad%4==0) detCoverage2a->Fill(zl,xl,1.0);
155	detThicknessZ2->Fill(zl,yl,1.0);
156	detThicknessX2->Fill(xl,yl,1.0);
157	break;
158	default:
159	} // end switch
160	} // end for i
161	} // end for m
162	ITS->ClearModules();
163	} // end for event
164	//
165	cout << "Creating pad"<<endl;
166	TVirtualPad *pad=0;
167	TCanvas *c0 = new TCanvas("c0","SPD Hits Tests",262,50,700,534);
168	c0->Range(0,0,1,1);
169	c0->Divide(3,2);
170	//
171	pad = c0->cd(1);
172	phiDistribution2->Draw();
173	pad = c0->cd(4);
174	phiDistribution1->Draw();
175	TPad *pad2 = c0->cd(2);
176	delete pad2;
177	c0->cd(0);
178	TPad *pad21 = new TPad("pad21","CoverageZ",0.34,0.875,0.65,1.0);
179	pad21->Draw();
180	pad21->cd();
181	pad21->Range(0,0,1,1);
182	pad21->Modified();
183	detCoverage2->Draw();
184	c0->Update();
185	c0->cd(0);
186	TPad *pad22 = new TPad("pad22","CoverageZa",0.34,0.75,0.65,0.875);
187	pad22->Draw();
188	pad22->cd();
189	pad22->Range(0,0,1,1);
190	pad22->Modified();
191	detCoverage2a->Draw();
192	c0->Update();
193	c0->cd(0);
194	TPad *pad23 = new TPad("pad23","ThicknessZ",0.34,0.625,0.65,0.75);
195	pad23->Draw();
196	pad23->cd();
197	pad23->Range(0,0,1,1);
198	pad23->Modified();
199	detThicknessZ2->Draw();
200	c0->Update();
201	c0->cd(0);
202	TPad *pad24 = new TPad("pad24","ThicknessX",0.34,0.50,0.65,0.625);
203	pad24->Draw();
204	pad24->cd();
205	pad24->Range(0,0,1,1);
206	pad24->Modified();
207	detThicknessX2->Draw();
208	c0->Update();
209	//
210	c0->cd(0);
211	TPad *pad5 = c0->cd(5);
212	delete pad5;
213	c0->cd(0);
214	TPad *pad51 = new TPad("pad51","CoverageZ",0.34,0.375,0.65,0.50);
215	pad51->Draw();
216	pad51->cd();
217	pad51->Range(0,0,1,1);
218	pad51->Modified();
219	detCoverage1->Draw();
220	c0->Update();
221	c0->cd(0);
222	TPad *pad52 = new TPad("pad52","CoverageZa",0.34,0.25,0.65,0.375);
223	pad52->Draw();
224	pad52->cd();
225	pad52->Range(0,0,1,1);
226	pad52->Modified();
227	detCoverage1a->Draw();
228	c0->Update();
229	c0->cd(0);
230	TPad *pad53 = new TPad("pad53","ThicknessZ",0.34,0.125,0.65,0.25);
231	pad53->Draw();
232	pad53->cd();
233	pad53->Range(0,0,1,1);
234	pad53->Modified();
235	detThicknessZ1->Draw();
236	c0->Update();
237	c0->cd(0);
238	TPad *pad54 = new TPad("pad54","ThicknessX",0.34,0.00,0.65,0.125);
239	pad54->Draw();
240	pad54->cd();
241	pad54->Range(0,0,1,1);
242	pad54->Modified();
243	detThicknessX1->Draw();
244	c0->Update();
245	//
246	TList *tvTreeList = new TList;
247	TTree tvTree = (TTree ) gROOT->FindObject("TreeH");
248	tvTreeList->Add(tvTree);
249	pad = c0->cd(3);
250	tvTree->Draw("ITS.fY:ITS.fX","ITS.fModule<240","", 100000, 0);
251	pad = c0->cd(6);
252	tvTree->Draw("ITS.fXITS.fX+ITS.fYITS.fY:ITS.fZ","ITS.fModule<240","",
253	100000, 0);
254
880d6abe	255	}