[u/mrichter/AliRoot.git] / TPC / CalibMacros / CreateGainMap.C

/*

This macro creates a gain map for the TPC based on the results of the Krypton calibration.
The main steps are the following:

1. Define outlier-pads where the krypton calibration was not succesful
2. A parabolic fit for the whole chamber is performed
3. replace outliers with fitted values
4. normalize separately IROCs and OROCs

For more details see below.


example usage:
==============

TFile f("calibKr.root")
AliTPCCalPad * kryptonRaw = new AliTPCCalPad(*fitMean)
AliTPCCalPad * kryptonMean = new AliTPCCalPad(*spectrMean)
AliTPCCalPad * kryptonChi2 = new AliTPCCalPad(*fitNormChi2)
AliTPCCalPad * kryptonRMS = new AliTPCCalPad(*fitRMS)

.L CreateGainMap.C
AliTPCCalPad * final = CreateGainMap(kryptonRaw, kryptonRMS)

TFile *h = new TFile("GainMap.root", "RECREATE")
final.Write()

*/


AliTPCCalPad * CreateGainMap(AliTPCCalPad *krypFitMean, AliTPCCalPad *krypFitRMS, AliTPCCalPad *noiseMap = 0, AliTPCCalPad *krypSpectrMean = 0, 
			     AliTPCCalPad *krypChi2 = 0, AliTPCCalPad *pulser = 0, AliTPCCalPad *electrode = 0) {

  
  // Draw input map
  TCanvas *test3 = new TCanvas("ASIDE3", "Aoriginal");
  krypFitMean->MakeHisto2D()->Draw("colz");
  TCanvas *test4 = new TCanvas("CSIDE4", "Coriginal");
  krypFitMean->MakeHisto2D(1)->Draw("colz");
  TH1F * hDEBUG = new TH1F("bla", "fitRMS/fitMean",100,0,0.3);

  const Double_t kryptonMean = 0.05012;
  const Double_t kryptonSigma = 0.00386;

  TObjArray arrayRocFinal(72);

  //
  // Loop over all sectors
  //
  for(Int_t isector=0; isector < 72; isector++){

    AliTPCCalROC *rocKrypFitMean = krypFitMean->GetCalROC(isector);
    AliTPCCalROC *rocKrypFitRMS  = krypFitRMS->GetCalROC(isector);

    AliTPCCalROC *rocOutlier = new AliTPCCalROC(*rocKrypFitMean);

    // 1. define map of outliers: the 41keV peak is fitted and krypFitMean represents the peak position and krypFitRMS its sigma.
    //    In order to control the fit qualtiy krypFitRMS/krypFitMean is computed and plotted; it is roughly gaussian distributed
    //    with the mean resolution kryptonMean = 0.05012 and sigma kryptonSigma = 0.00386. If the deviation is larger than 4sigma
    //    the fit to the 41keV peak is considered as failed.

    Int_t nPointsFit = 0;
    for(Int_t iChannel=0; iChannel < rocOutlier->GetNchannels(); iChannel++){
      Double_t fitMean = rocKrypFitMean->GetValue(iChannel);
      Double_t fitRMS = rocKrypFitRMS->GetValue(iChannel);
      if (fitRMS < 0.001 || fitMean < 0.001) continue;
       hDEBUG->Fill(fitRMS/fitMean);
      if (TMath::Abs(fitRMS/fitMean - kryptonMean) < 4*kryptonSigma || fitRMS < 0.001 || fitMean < 0.001) {
	rocOutlier->SetValue(iChannel,0);
	nPointsFit++;
      }
    }
    //TCanvas *DEBUG = new TCanvas();
    //rocOutlier->MakeHisto2D()->Draw("colz");
    
    Double_t rocMedian = rocKrypFitMean->GetMedian(rocOutlier);
    

    // 2. make a parabolic fit for the whole chamber with excluded outliers
    
    TVectorD params;
    TMatrixD cov;
    Float_t chi2;
    rocKrypFitMean->GlobalFit(rocOutlier,0,params,cov,chi2);
    AliTPCCalROC *rocParabolicFit=AliTPCCalROC::CreateGlobalFitCalROC(params,isector);

    if (nPointsFit != 0) cout << "sector: "<< isector << " chi2: " << chi2/nPointsFit <<" median: "<< rocMedian <<" n: " << nPointsFit << endl;
    if (nPointsFit == 0) {
      AliTPCCalROC *rocFinal = new AliTPCCalROC(*rocParabolicFit); // What to do with sectors being switched off??
      arrayRocFinal.AddAt(rocFinal,isector);
      continue;
    }

    //
    // VERY IMPORTANT **** VERY IMPORTANT  **** VERY IMPORTANT  **** VERY IMPORTANT  **** VERY IMPORTANT  **** VERY IMPORTANT 
    //
    // if the fit is considered as failed, the mean value is taken for the whole chamber (TO AVOID THIS, RELEASE THIS CUT)
    // if cosmic data can be used for gain calibration, this cut should be strong (e.g. 5) !!!!!!!!!!

    if (chi2/nPointsFit > 5) {
      AliTPCCalROC *rocFinal = new AliTPCCalROC(*rocParabolicFit);
      if (rocMedian == 0 && isector == 51) rocMedian = 1075; // manual patch to be removed for sector 51!
      for(Int_t iChannel=0; iChannel < rocFinal->GetNchannels(); iChannel++) rocFinal->SetValue(iChannel,rocMedian);
      arrayRocFinal.AddAt(rocFinal,isector);
      continue;
    }

    // 3. replace outliers with fitted values

    AliTPCCalROC *rocFinal = new AliTPCCalROC(*rocParabolicFit);
    for(Int_t iChannel=0; iChannel < rocFinal->GetNchannels(); iChannel++){
      Double_t fitMean = rocKrypFitMean->GetValue(iChannel);
      Double_t fitRMS = rocKrypFitRMS->GetValue(iChannel);   
      if (fitRMS < 0.001 || fitMean < 0.001 || TMath::Abs(rocParabolicFit->GetValue(iChannel)/fitMean - 1) > 0.35) {
	rocFinal->SetValue(iChannel,rocParabolicFit->GetValue(iChannel));
	continue;
      }
      if (TMath::Abs(fitRMS/fitMean - kryptonMean) < 4*kryptonSigma) {
	rocFinal->SetValue(iChannel,rocKrypFitMean->GetValue(iChannel));
      }
    }
    

    // 4. Postprocessing: Set dead channels and very noisy channels (time dependent) to 0
    
    const Double_t noiseMin = 0.01;
    const Double_t noiseMax = 2;

    if (noiseMap) {
      AliTPCCalROC *rocNoise = noiseMap->GetCalROC(isector);
      for(Int_t iChannel=0; iChannel < rocFinal->GetNchannels(); iChannel++){
	Double_t noise = rocNoise->GetValue(iChannel);
	if (noise < noiseMin || noise > noiseMax) rocFinal->SetValue(iChannel, 0);
      }
    }
    // Fill an array of ROCs

    arrayRocFinal.AddAt(rocFinal,isector);

  }

  AliTPCCalPad *final = new AliTPCCalPad(&arrayRocFinal);
  
  // 4. normalize separately IROCs and OROCs to the mean of their medians

  Double_t meanMedianIROC;
  Double_t meanMedianOROC;
  Int_t n = 0;

  for(Int_t isector=0; isector < 36; isector++){ // IROCs
    AliTPCCalROC *rocFinal = final->GetCalROC(isector);
    if (rocFinal->GetMedian() != 0) {
      meanMedianIROC += rocFinal->GetMedian();
      n++;
    }
  }
  meanMedianIROC = meanMedianIROC/n;

  n = 0;
  for(Int_t isector=36; isector < 72; isector++){ // OROCs
    AliTPCCalROC *rocFinal = final->GetCalROC(isector);
    if (rocFinal->GetMedian() != 0) {
      meanMedianOROC += rocFinal->GetMedian();
      n++;
    }
  }
  meanMedianOROC = meanMedianOROC/n;

  for(Int_t isector=0; isector < 72; isector++){ // OROCs
    AliTPCCalROC *rocFinal = final->GetCalROC(isector);
    if (isector<36) rocFinal->Multiply(1./meanMedianIROC);
    if (isector>35) rocFinal->Multiply(1./meanMedianOROC);
  }

  // Draw results
  TCanvas *test = new TCanvas("ASIDE", "A");
  final->MakeHisto2D()->Draw("colz");
  TCanvas *test2 = new TCanvas("CSIDE", "C");
  final->MakeHisto2D(1)->Draw("colz");
  TCanvas *cDEBUG = new TCanvas();
  hDEBUG->Draw();

  //return results
  final->SetName("GainMap");
  return final;

}
Commit	Line	Data
6ba6d650	1	/*
	2
	3	This macro creates a gain map for the TPC based on the results of the Krypton calibration.
	4	The main steps are the following:
	5
	6	1. Define outlier-pads where the krypton calibration was not succesful
	7	2. A parabolic fit for the whole chamber is performed
	8	3. replace outliers with fitted values
	9	4. normalize separately IROCs and OROCs
	10
	11	For more details see below.
	12
	13
	14
	15	example usage:
	16	==============
	17
	18	TFile f("calibKr.root")
	19	AliTPCCalPad * kryptonRaw = new AliTPCCalPad(*fitMean)
	20	AliTPCCalPad * kryptonMean = new AliTPCCalPad(*spectrMean)
	21	AliTPCCalPad * kryptonChi2 = new AliTPCCalPad(*fitNormChi2)
	22	AliTPCCalPad * kryptonRMS = new AliTPCCalPad(*fitRMS)
	23
	24	.L CreateGainMap.C
	25	AliTPCCalPad * final = CreateGainMap(kryptonRaw, kryptonRMS)
	26
	27	TFile *h = new TFile("GainMap.root", "RECREATE")
	28	final.Write()
	29
	30	*/
	31
	32
	33	AliTPCCalPad * CreateGainMap(AliTPCCalPad krypFitMean, AliTPCCalPad krypFitRMS, AliTPCCalPad noiseMap = 0, AliTPCCalPad krypSpectrMean = 0,
	34	AliTPCCalPad krypChi2 = 0, AliTPCCalPad pulser = 0, AliTPCCalPad *electrode = 0) {
	35
	36
	37	// Draw input map
	38	TCanvas *test3 = new TCanvas("ASIDE3", "Aoriginal");
	39	krypFitMean->MakeHisto2D()->Draw("colz");
	40	TCanvas *test4 = new TCanvas("CSIDE4", "Coriginal");
	41	krypFitMean->MakeHisto2D(1)->Draw("colz");
	42	TH1F * hDEBUG = new TH1F("bla", "fitRMS/fitMean",100,0,0.3);
	43
	44	const Double_t kryptonMean = 0.05012;
	45	const Double_t kryptonSigma = 0.00386;
	46
	47	TObjArray arrayRocFinal(72);
	48
	49	//
	50	// Loop over all sectors
	51	//
	52	for(Int_t isector=0; isector < 72; isector++){
	53
	54	AliTPCCalROC *rocKrypFitMean = krypFitMean->GetCalROC(isector);
	55	AliTPCCalROC *rocKrypFitRMS = krypFitRMS->GetCalROC(isector);
	56
	57	AliTPCCalROC rocOutlier = new AliTPCCalROC(rocKrypFitMean);
	58
	59	// 1. define map of outliers: the 41keV peak is fitted and krypFitMean represents the peak position and krypFitRMS its sigma.
	60	// In order to control the fit qualtiy krypFitRMS/krypFitMean is computed and plotted; it is roughly gaussian distributed
	61	// with the mean resolution kryptonMean = 0.05012 and sigma kryptonSigma = 0.00386. If the deviation is larger than 4sigma
	62	// the fit to the 41keV peak is considered as failed.
	63
	64	Int_t nPointsFit = 0;
65	for(Int_t iChannel=0; iChannel < rocOutlier->GetNchannels(); iChannel++){
66	Double_t fitMean = rocKrypFitMean->GetValue(iChannel);
67	Double_t fitRMS = rocKrypFitRMS->GetValue(iChannel);
68	if (fitRMS < 0.001 \|\| fitMean < 0.001) continue;
69	hDEBUG->Fill(fitRMS/fitMean);
70	if (TMath::Abs(fitRMS/fitMean - kryptonMean) < 4*kryptonSigma \|\| fitRMS < 0.001 \|\| fitMean < 0.001) {
71	rocOutlier->SetValue(iChannel,0);
72	nPointsFit++;
73	}
74	}
75	//TCanvas *DEBUG = new TCanvas();
76	//rocOutlier->MakeHisto2D()->Draw("colz");
77
78	Double_t rocMedian = rocKrypFitMean->GetMedian(rocOutlier);
79
80
81	// 2. make a parabolic fit for the whole chamber with excluded outliers
82
83	TVectorD params;
84	TMatrixD cov;
85	Float_t chi2;
86	rocKrypFitMean->GlobalFit(rocOutlier,0,params,cov,chi2);
87	AliTPCCalROC *rocParabolicFit=AliTPCCalROC::CreateGlobalFitCalROC(params,isector);
88
89	if (nPointsFit != 0) cout << "sector: "<< isector << " chi2: " << chi2/nPointsFit <<" median: "<< rocMedian <<" n: " << nPointsFit << endl;
90	if (nPointsFit == 0) {
91	AliTPCCalROC rocFinal = new AliTPCCalROC(rocParabolicFit); // What to do with sectors being switched off??
92	arrayRocFinal.AddAt(rocFinal,isector);
93	continue;
94	}
95
96	//
97	// VERY IMPORTANT ** VERY IMPORTANT VERY IMPORTANT VERY IMPORTANT VERY IMPORTANT ** VERY IMPORTANT
98	//
99	// if the fit is considered as failed, the mean value is taken for the whole chamber (TO AVOID THIS, RELEASE THIS CUT)
100	// if cosmic data can be used for gain calibration, this cut should be strong (e.g. 5) !!!!!!!!!!
101
102	if (chi2/nPointsFit > 5) {
103	AliTPCCalROC rocFinal = new AliTPCCalROC(rocParabolicFit);
104	if (rocMedian == 0 && isector == 51) rocMedian = 1075; // manual patch to be removed for sector 51!
105	for(Int_t iChannel=0; iChannel < rocFinal->GetNchannels(); iChannel++) rocFinal->SetValue(iChannel,rocMedian);
106	arrayRocFinal.AddAt(rocFinal,isector);
107	continue;
108	}
109
110	// 3. replace outliers with fitted values
111
112	AliTPCCalROC rocFinal = new AliTPCCalROC(rocParabolicFit);
113	for(Int_t iChannel=0; iChannel < rocFinal->GetNchannels(); iChannel++){
114	Double_t fitMean = rocKrypFitMean->GetValue(iChannel);
115	Double_t fitRMS = rocKrypFitRMS->GetValue(iChannel);
116	if (fitRMS < 0.001 \|\| fitMean < 0.001 \|\| TMath::Abs(rocParabolicFit->GetValue(iChannel)/fitMean - 1) > 0.35) {
117	rocFinal->SetValue(iChannel,rocParabolicFit->GetValue(iChannel));
118	continue;
119	}
120	if (TMath::Abs(fitRMS/fitMean - kryptonMean) < 4*kryptonSigma) {
121	rocFinal->SetValue(iChannel,rocKrypFitMean->GetValue(iChannel));
122	}
123	}
124
125
126	// 4. Postprocessing: Set dead channels and very noisy channels (time dependent) to 0
127
128	const Double_t noiseMin = 0.01;
129	const Double_t noiseMax = 2;
130
131	if (noiseMap) {
132	AliTPCCalROC *rocNoise = noiseMap->GetCalROC(isector);
133	for(Int_t iChannel=0; iChannel < rocFinal->GetNchannels(); iChannel++){
134	Double_t noise = rocNoise->GetValue(iChannel);
135	if (noise < noiseMin \|\| noise > noiseMax) rocFinal->SetValue(iChannel, 0);
136	}
137	}
138	// Fill an array of ROCs
139
140	arrayRocFinal.AddAt(rocFinal,isector);
141
142	}
143
144	AliTPCCalPad *final = new AliTPCCalPad(&arrayRocFinal);
145
146	// 4. normalize separately IROCs and OROCs to the mean of their medians
147
148	Double_t meanMedianIROC;
149	Double_t meanMedianOROC;
150	Int_t n = 0;
151
152	for(Int_t isector=0; isector < 36; isector++){ // IROCs
153	AliTPCCalROC *rocFinal = final->GetCalROC(isector);
154	if (rocFinal->GetMedian() != 0) {
155	meanMedianIROC += rocFinal->GetMedian();
156	n++;
157	}
158	}
159	meanMedianIROC = meanMedianIROC/n;
160
161	n = 0;
162	for(Int_t isector=36; isector < 72; isector++){ // OROCs
163	AliTPCCalROC *rocFinal = final->GetCalROC(isector);
164	if (rocFinal->GetMedian() != 0) {
165	meanMedianOROC += rocFinal->GetMedian();
166	n++;
167	}
168	}
169	meanMedianOROC = meanMedianOROC/n;
170
171	for(Int_t isector=0; isector < 72; isector++){ // OROCs
172	AliTPCCalROC *rocFinal = final->GetCalROC(isector);
173	if (isector<36) rocFinal->Multiply(1./meanMedianIROC);
174	if (isector>35) rocFinal->Multiply(1./meanMedianOROC);
175	}
176
177	// Draw results
178	TCanvas *test = new TCanvas("ASIDE", "A");
179	final->MakeHisto2D()->Draw("colz");
180	TCanvas *test2 = new TCanvas("CSIDE", "C");
181	final->MakeHisto2D(1)->Draw("colz");
182	TCanvas *cDEBUG = new TCanvas();
183	hDEBUG->Draw();
184
185	//return results
186	final->SetName("GainMap");
187	return final;
188
189	}
190