[u/mrichter/AliRoot.git] / TRD / Macros / AliTRDanalyzeTestBeam.C


void AliTRDanalyzeTestBeam(Int_t run, Int_t begin, Int_t end) {
  
  gROOT->SetStyle("Plain");
  //gStyle->SetPadTopMargin(0.02);
  //gStyle->SetPadRightMargin(0.02);
  //gStyle->SetPadLeftMargin(0.1);

  gStyle->SetPalette(1);
  gStyle->SetOptStat(0);
  gStyle->SetOptDate();
  
  TGaxis::SetMaxDigits(3);

  TStopwatch st;
  st.Start();

  const Int_t N = 640;

  // declare histograms
  const int nbins = 100;
  const int nstart = 0;

  TH1D *mSi1L = new TH1D("mSi1L", ";number of Si1 fired pads", 50, nstart-0.5, nstart+nbins-0.5);
  TH1D *mSi2L = new TH1D("mS21L", ";number of Si2 fired pads", 50, nstart-0.5, nstart+nbins-0.5);
 
  TProfile *mSi1ChP = new TProfile("mSi1ChP",";Si1 pad number;signal", 1280, -0.5, 1279.5, 0, 200, "s");
  TProfile *mSi2ChP = new TProfile("mS21ChP",";Si2 pad number;signal", 1280, -0.5, 1279.5, 0, 200, "s");
  
  TH1D *mSi1N = new TH1D("mSi1N", "Noise Dist Si1;ADC", 100, 0, 50);
  TH1D *mSi2N = new TH1D("mSi2N", "Noise Dist Si2;ADC", 100, 0, 50);
  
  TH1D *mSiCh[4];
  mSiCh[0] = new TH1D("mSi1ChX", ";Si1X pad amplitude (ADC)", 250, -0.5, 249.5);
  mSiCh[1] = new TH1D("mSi1ChY", ";Si1Y pad amplitude (ADC)", 250, -0.5, 249.5);
  mSiCh[2] = new TH1D("mSi2ChX", ";Si2X pad amplitude (ADC)", 250, -0.5, 249.5);
  mSiCh[3] = new TH1D("mSi2ChY", ";Si2Y pad amplitude (ADC)", 250, -0.5, 249.5);
  
  TH1D *mSiFullCh[4];
  mSiFullCh[0] = new TH1D("mSi1fChX", "Si1X;max amplitude (ADC)", 300, -0.5, 299.5);
  mSiFullCh[1] = new TH1D("mSi1fChY", "Si1Y;max amplitude (ADC)", 300, -0.5, 299.5);
  mSiFullCh[2] = new TH1D("mSi2fChX", "Si2X;max amplitude (ADC)", 300, -0.5, 299.5);
  mSiFullCh[3] = new TH1D("mSi2fChY", "Si2Y;max amplitude (ADC)", 300, -0.5, 299.5);  

  TH2D *mPos[2];
  mPos[0] = new TH2D("posSi1", ";x Si1 (mm);y Si1 (mm)", 128, 0, 32, 128, 0, 32);
  mPos[1] = new TH2D("posSi2", ";x Si2 (mm);y Si2 (mm)", 128, 0, 32, 128, 0, 32);

  TH2D *mCor[2];
  mCor[0] = new TH2D("corX", ";Si1 X (mm);Si2 X (mm)", 128, 0, 32, 128, 0, 32);
  mCor[1] = new TH2D("corY", ";Si1 Y (mm);Si2 Y (mm)", 128, 0, 32, 128, 0, 32);

  TH2D *mChCor[2];
  mChCor[0] = new TH2D("ChCorSi1", ";Si1 amp X;Si1 amp Y", 100, 0, 200, 100, 0, 200);
  mChCor[1] = new TH2D("ChCorSi2", ";Si2 amp X;Si2 amp Y", 100, 0, 200, 100, 0, 200);
  
  gStyle->SetOptStat(11);
  TH1D *mPb   = new TH1D("mPb", ";Amp Pb (ADC)", 150, -0.5, 4499.5);
  TH1D *mCher = new TH1D("mCher", ";Amp Cherenkov (ADC)", 150, -0.5, 4499.5);
  TH2D *mPbCher = new TH2D("mPbCher", ";amp Cherenkov;amp Pb", 150, -0.5, 4499.5, 150, -0.5, 4599.5);
  // gStyle->SetOptStat(0);

  // needed by the AliTRDRawStreamTB
  AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
  AliTRDcalibDB *calib = AliTRDcalibDB::Instance();
  calib->SetRun(0);

  // TRD data monitoring
  TH1D *mDet = new TH1D("mDet", ";chamber", 20, -0.5, 19.5);
  TH1D *mROB = new TH1D("mROB", ";rob", 20, -0.5, 19.5);
  TH1D *mTRDsig = new TH1D("mTRDsig", ";trdSignal", 100, -0.5, 299.5);
  
  //AliLog::SetClassDebugLevel("AliTRDRawStreamTB", 10);

  int counter = 0;
  // for(Int_t run = 365; run < 386; run++) {
  //for(Int_t run = 369; run < 382; run++) {
  //for(int run = 387; run < 389; run++) {

  // if (run == 389) continue;
    cout << run << endl;

    for(Int_t fn=begin; fn<end; fn++) {
      
      // connect to data 
      //const char *base="/Users/radomski/data/1GeV/";
      const char *base = "./";
      const char *filename = Form("%s/run%d_gdc_daq09.%03d.raw",base,run,fn);
      if (gSystem->AccessPathName(filename)) continue;
      cout << filename << endl;

      AliTRDtestBeam *data = new AliTRDtestBeam(filename);
      
      // process data
      while (data->NextEvent()) {
      
	if (!(counter%1000)) cout << "Event = " << counter << endl;
	counter++;
	
	/*
	AliTRDRawStreamTB *tb = data->GetTRDrawStream();
	while(tb->Next()) {
	  mROB->Fill(tb->GetROB());
	  mDet->Fill(tb->GetDet());
	  int *sig = tb->GetSignals();
	  mTRDsig->Fill(sig[0]);
	  mTRDsig->Fill(sig[1]);
	  mTRDsig->Fill(sig[2]);
	}
	delete tb;
	*/

	mCher->Fill(data->GetCher());
	mPb->Fill(data->GetPb());
	mPbCher->Fill(data->GetCher(), data->GetPb());

	mSi1L->Fill(data->GetNSi1());
	mSi2L->Fill(data->GetNSi2());
	
	for(int i=0; i<data->GetNSi1(); i++) {
	  Int_t q = data->GetSi1Charge(i);
	  Int_t a = data->GetSi1Address(i);
	  if (a == 0) continue; // noisy channels
	  mSi1ChP->Fill(a, q);
	  if (a < N) mSiCh[0]->Fill(q);
	  else mSiCh[1]->Fill(q);	       
	}
	
	for(int i=0; i<data->GetNSi2(); i++) {
	  Int_t q = data->GetSi2Charge(i);
	  Int_t a = data->GetSi2Address(i);
	  if (a == 0 || a == 1279) continue; // noisy channels
	  mSi2ChP->Fill(a, q);
	  if (a < N) mSiCh[2]->Fill(q);
	  else mSiCh[3]->Fill(q);	       
	}
	
	mSiFullCh[0]->Fill(data->GetQx(0));
	mSiFullCh[1]->Fill(data->GetQy(0));
	mSiFullCh[2]->Fill(data->GetQx(1));
	mSiFullCh[3]->Fill(data->GetQy(1));
	
	for(int k=0; k<2; k++)
	  mChCor[k]->Fill(data->GetQx(k), data->GetQy(k));
	
	/*
	  if (data->GetQx(0) < 20) continue;
	  if (data->GetQx(1) < 20) continue;
	  if (data->GetQy(0) < 20) continue;
	  if (data->GetQy(1) < 20) continue;
	*/
	
	for(int k=0; k<2; k++)
	  mPos[k]->Fill(data->GetX(k), data->GetY(k));
      
	mCor[0]->Fill(data->GetX(0), data->GetX(1));
	mCor[1]->Fill(data->GetY(0), data->GetY(1));
      }
      
      delete data;
    }
    //}

  // process histograms
  for(int i=1; i<1281; i++) mSi1N->Fill(mSi1ChP->GetBinError(i));
  for(int i=1; i<1281; i++) mSi2N->Fill(mSi2ChP->GetBinError(i));

  // display
  cout << "Number of Events = " << counter << endl;
  
  /**/
  TCanvas *c = new TCanvas("siliconSignal", "silicon signal");
  c->Divide(2,2, 0.01, 0.01);
  c->cd(1);
  mSi1L->Draw();
  c->cd(2);
  mSi2L->Draw();
  c->cd(3);
  mSi1ChP->Draw();
  c->cd(4);
  mSi2ChP->Draw();
  /* */

  /**/
  c = new TCanvas();
  c->Divide(2,2,0.01,0.01);
  c->cd(1);
  mSi1N->Draw();
  c->cd(2);
  mSi2N->Draw();
  /**/

  /**/
  // pads
  c = new TCanvas("siPads", "Silicon Pads");
  c->Divide(2,2, 0.01, 0.01);
  for(int i=0; i<4; i++) {
    c->cd(i+1);
    gPad->SetLogy();
    mSiCh[i]->Draw();
  }
  
  // clusters
  c = new TCanvas("siCluster", "silicon clusters");
  c->Divide(2,2, 0.01, 0.01);
  for(int i=0; i<4; i++) {
    c->cd(i+1);
    gPad->SetLogy();
    mSiFullCh[i]->Draw();
  }

  // position and correlation
  c = new TCanvas("siPosition", "reconstructed position");
  c->Divide(2,2, 0.01, 0.01);
  for(int i=0; i<2; i++) {
    c->cd(1+i);
    mPos[i]->Draw("col");
    c->cd(3+i);
    mCor[i]->Draw("col");  
  }
  
  c = new TCanvas("siCharge", "si charge correlation");
  c->Divide(2,2, 0.01, 0.01);
  c->cd(1);
  gPad->SetLogz();
  mChCor[0]->Draw("col");
  c->cd(2);
  gPad->SetLogz();
  mChCor[1]->Draw("col");
  /**/

  new TCanvas();
  gPad->SetLogy();
  mCher->Draw();
  
  // electron sample
  int bin = mCher->FindBin(500.);
  double ef = (mCher->Integral(bin, 151)/ mCher->GetSum());
  TLatex *l = new TLatex(2e3, 0.02*mCher->GetSum(), Form("Electron fraction = %.2f ", ef));
  l->Draw();
  
  new TCanvas();
  gPad->SetLogy();
  mPb->Draw();

  new TCanvas();
  gPad->SetLogz();
  mPbCher->Draw("colz");
  
  /*
  c = new TCanvas();
  c->Divide(2,2,0.01, 0.01);
  c->cd(1);
  gPad->SetLogy();
  mTRDsig->Draw();

  c->cd(2);
  mROB->Draw();

  c->cd(3);
  mDet->Draw();
  */
  
  st.Stop();
  st.Print();
}

//Int_t SelectEvent() {
//}
Commit	Line	Data
1a6c0f30	1
	2	void AliTRDanalyzeTestBeam(Int_t run, Int_t begin, Int_t end) {
	3
	4	gROOT->SetStyle("Plain");
	5	//gStyle->SetPadTopMargin(0.02);
	6	//gStyle->SetPadRightMargin(0.02);
	7	//gStyle->SetPadLeftMargin(0.1);
	8
	9	gStyle->SetPalette(1);
	10	gStyle->SetOptStat(0);
	11	gStyle->SetOptDate();
	12
	13	TGaxis::SetMaxDigits(3);
	14
	15	TStopwatch st;
	16	st.Start();
	17
	18	const Int_t N = 640;
	19
	20	// declare histograms
	21	const int nbins = 100;
	22	const int nstart = 0;
	23
	24	TH1D *mSi1L = new TH1D("mSi1L", ";number of Si1 fired pads", 50, nstart-0.5, nstart+nbins-0.5);
	25	TH1D *mSi2L = new TH1D("mS21L", ";number of Si2 fired pads", 50, nstart-0.5, nstart+nbins-0.5);
	26
	27	TProfile *mSi1ChP = new TProfile("mSi1ChP",";Si1 pad number;signal", 1280, -0.5, 1279.5, 0, 200, "s");
	28	TProfile *mSi2ChP = new TProfile("mS21ChP",";Si2 pad number;signal", 1280, -0.5, 1279.5, 0, 200, "s");
	29
	30	TH1D *mSi1N = new TH1D("mSi1N", "Noise Dist Si1;ADC", 100, 0, 50);
	31	TH1D *mSi2N = new TH1D("mSi2N", "Noise Dist Si2;ADC", 100, 0, 50);
	32
	33	TH1D *mSiCh[4];
	34	mSiCh[0] = new TH1D("mSi1ChX", ";Si1X pad amplitude (ADC)", 250, -0.5, 249.5);
	35	mSiCh[1] = new TH1D("mSi1ChY", ";Si1Y pad amplitude (ADC)", 250, -0.5, 249.5);
	36	mSiCh[2] = new TH1D("mSi2ChX", ";Si2X pad amplitude (ADC)", 250, -0.5, 249.5);
	37	mSiCh[3] = new TH1D("mSi2ChY", ";Si2Y pad amplitude (ADC)", 250, -0.5, 249.5);
	38
	39	TH1D *mSiFullCh[4];
	40	mSiFullCh[0] = new TH1D("mSi1fChX", "Si1X;max amplitude (ADC)", 300, -0.5, 299.5);
	41	mSiFullCh[1] = new TH1D("mSi1fChY", "Si1Y;max amplitude (ADC)", 300, -0.5, 299.5);
	42	mSiFullCh[2] = new TH1D("mSi2fChX", "Si2X;max amplitude (ADC)", 300, -0.5, 299.5);
	43	mSiFullCh[3] = new TH1D("mSi2fChY", "Si2Y;max amplitude (ADC)", 300, -0.5, 299.5);
	44
	45	TH2D *mPos[2];
	46	mPos[0] = new TH2D("posSi1", ";x Si1 (mm);y Si1 (mm)", 128, 0, 32, 128, 0, 32);
	47	mPos[1] = new TH2D("posSi2", ";x Si2 (mm);y Si2 (mm)", 128, 0, 32, 128, 0, 32);
	48
	49	TH2D *mCor[2];
	50	mCor[0] = new TH2D("corX", ";Si1 X (mm);Si2 X (mm)", 128, 0, 32, 128, 0, 32);
	51	mCor[1] = new TH2D("corY", ";Si1 Y (mm);Si2 Y (mm)", 128, 0, 32, 128, 0, 32);
	52
	53	TH2D *mChCor[2];
	54	mChCor[0] = new TH2D("ChCorSi1", ";Si1 amp X;Si1 amp Y", 100, 0, 200, 100, 0, 200);
	55	mChCor[1] = new TH2D("ChCorSi2", ";Si2 amp X;Si2 amp Y", 100, 0, 200, 100, 0, 200);
	56
	57	gStyle->SetOptStat(11);
	58	TH1D *mPb = new TH1D("mPb", ";Amp Pb (ADC)", 150, -0.5, 4499.5);
	59	TH1D *mCher = new TH1D("mCher", ";Amp Cherenkov (ADC)", 150, -0.5, 4499.5);
	60	TH2D *mPbCher = new TH2D("mPbCher", ";amp Cherenkov;amp Pb", 150, -0.5, 4499.5, 150, -0.5, 4599.5);
	61	// gStyle->SetOptStat(0);
	62
	63	// needed by the AliTRDRawStreamTB
162637e4	64	AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
1a6c0f30	65	AliTRDcalibDB *calib = AliTRDcalibDB::Instance();
	66	calib->SetRun(0);
	67
	68	// TRD data monitoring
	69	TH1D *mDet = new TH1D("mDet", ";chamber", 20, -0.5, 19.5);
	70	TH1D *mROB = new TH1D("mROB", ";rob", 20, -0.5, 19.5);
	71	TH1D *mTRDsig = new TH1D("mTRDsig", ";trdSignal", 100, -0.5, 299.5);
	72
	73	//AliLog::SetClassDebugLevel("AliTRDRawStreamTB", 10);
	74
	75	int counter = 0;
	76	// for(Int_t run = 365; run < 386; run++) {
	77	//for(Int_t run = 369; run < 382; run++) {
	78	//for(int run = 387; run < 389; run++) {
	79
	80	// if (run == 389) continue;
	81	cout << run << endl;
	82
	83	for(Int_t fn=begin; fn<end; fn++) {
	84
	85	// connect to data
	86	//const char *base="/Users/radomski/data/1GeV/";
	87	const char *base = "./";
	88	const char *filename = Form("%s/run%d_gdc_daq09.%03d.raw",base,run,fn);
	89	if (gSystem->AccessPathName(filename)) continue;
	90	cout << filename << endl;
	91
	92	AliTRDtestBeam *data = new AliTRDtestBeam(filename);
	93
	94	// process data
	95	while (data->NextEvent()) {
	96
	97	if (!(counter%1000)) cout << "Event = " << counter << endl;
	98	counter++;
	99
	100	/*
	101	AliTRDRawStreamTB *tb = data->GetTRDrawStream();
	102	while(tb->Next()) {
	103	mROB->Fill(tb->GetROB());
	104	mDet->Fill(tb->GetDet());
	105	int *sig = tb->GetSignals();
	106	mTRDsig->Fill(sig[0]);
	107	mTRDsig->Fill(sig[1]);
	108	mTRDsig->Fill(sig[2]);
	109	}
	110	delete tb;
	111	*/
	112
	113	mCher->Fill(data->GetCher());
	114	mPb->Fill(data->GetPb());
	115	mPbCher->Fill(data->GetCher(), data->GetPb());
	116
	117	mSi1L->Fill(data->GetNSi1());
	118	mSi2L->Fill(data->GetNSi2());
	119
	120	for(int i=0; i<data->GetNSi1(); i++) {
	121	Int_t q = data->GetSi1Charge(i);
	122	Int_t a = data->GetSi1Address(i);
	123	if (a == 0) continue; // noisy channels
	124	mSi1ChP->Fill(a, q);
	125	if (a < N) mSiCh[0]->Fill(q);
	126	else mSiCh[1]->Fill(q);
	127	}
	128
129	for(int i=0; i<data->GetNSi2(); i++) {
130	Int_t q = data->GetSi2Charge(i);
131	Int_t a = data->GetSi2Address(i);
132	if (a == 0 \|\| a == 1279) continue; // noisy channels
133	mSi2ChP->Fill(a, q);
134	if (a < N) mSiCh[2]->Fill(q);
135	else mSiCh[3]->Fill(q);
136	}
137
138	mSiFullCh[0]->Fill(data->GetQx(0));
139	mSiFullCh[1]->Fill(data->GetQy(0));
140	mSiFullCh[2]->Fill(data->GetQx(1));
141	mSiFullCh[3]->Fill(data->GetQy(1));
142
143	for(int k=0; k<2; k++)
144	mChCor[k]->Fill(data->GetQx(k), data->GetQy(k));
145
146	/*
147	if (data->GetQx(0) < 20) continue;
148	if (data->GetQx(1) < 20) continue;
149	if (data->GetQy(0) < 20) continue;
150	if (data->GetQy(1) < 20) continue;
151	*/
152
153	for(int k=0; k<2; k++)
154	mPos[k]->Fill(data->GetX(k), data->GetY(k));
155
156	mCor[0]->Fill(data->GetX(0), data->GetX(1));
157	mCor[1]->Fill(data->GetY(0), data->GetY(1));
158	}
159
160	delete data;
161	}
162	//}
163
164	// process histograms
165	for(int i=1; i<1281; i++) mSi1N->Fill(mSi1ChP->GetBinError(i));
166	for(int i=1; i<1281; i++) mSi2N->Fill(mSi2ChP->GetBinError(i));
167
168	// display
169	cout << "Number of Events = " << counter << endl;
170
171	/**/
172	TCanvas *c = new TCanvas("siliconSignal", "silicon signal");
173	c->Divide(2,2, 0.01, 0.01);
174	c->cd(1);
175	mSi1L->Draw();
176	c->cd(2);
177	mSi2L->Draw();
178	c->cd(3);
179	mSi1ChP->Draw();
180	c->cd(4);
181	mSi2ChP->Draw();
182	/* */
183
184	/**/
185	c = new TCanvas();
186	c->Divide(2,2,0.01,0.01);
187	c->cd(1);
188	mSi1N->Draw();
189	c->cd(2);
190	mSi2N->Draw();
191	/**/
192
193	/**/
194	// pads
195	c = new TCanvas("siPads", "Silicon Pads");
196	c->Divide(2,2, 0.01, 0.01);
197	for(int i=0; i<4; i++) {
198	c->cd(i+1);
199	gPad->SetLogy();
200	mSiCh[i]->Draw();
201	}
202
203	// clusters
204	c = new TCanvas("siCluster", "silicon clusters");
205	c->Divide(2,2, 0.01, 0.01);
206	for(int i=0; i<4; i++) {
207	c->cd(i+1);
208	gPad->SetLogy();
209	mSiFullCh[i]->Draw();
210	}
211
212	// position and correlation
213	c = new TCanvas("siPosition", "reconstructed position");
214	c->Divide(2,2, 0.01, 0.01);
215	for(int i=0; i<2; i++) {
216	c->cd(1+i);
217	mPos[i]->Draw("col");
218	c->cd(3+i);
219	mCor[i]->Draw("col");
220	}
221
222	c = new TCanvas("siCharge", "si charge correlation");
223	c->Divide(2,2, 0.01, 0.01);
224	c->cd(1);
225	gPad->SetLogz();
226	mChCor[0]->Draw("col");
227	c->cd(2);
228	gPad->SetLogz();
229	mChCor[1]->Draw("col");
230	/**/
231
232	new TCanvas();
233	gPad->SetLogy();
234	mCher->Draw();
235
236	// electron sample
237	int bin = mCher->FindBin(500.);
238	double ef = (mCher->Integral(bin, 151)/ mCher->GetSum());
239	TLatex l = new TLatex(2e3, 0.02mCher->GetSum(), Form("Electron fraction = %.2f ", ef));
240	l->Draw();
241
242	new TCanvas();
243	gPad->SetLogy();
244	mPb->Draw();
245
246	new TCanvas();
247	gPad->SetLogz();
248	mPbCher->Draw("colz");
249
250	/*
251	c = new TCanvas();
252	c->Divide(2,2,0.01, 0.01);
253	c->cd(1);
254	gPad->SetLogy();
255	mTRDsig->Draw();
256
257	c->cd(2);
258	mROB->Draw();
259
260	c->cd(3);
261	mDet->Draw();
262	*/
263
264	st.Stop();
265	st.Print();
266	}
267
268	//Int_t SelectEvent() {
269	//}