[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;

}


void MakeCalibTree(char * inputKr="calibKr.root", char * inputCE ="fitCE.root", char * inputPulser=0){
  //
  //
  //
   AliTPCPreprocessorOnline * preprocesor = new AliTPCPreprocessorOnline;
   TFile f(inputKr);
   TFile fce(inputCE);
   AliTPCCalPad * kryptonMean     = (AliTPCCalPad*)f.Get("spectrMean");
   AliTPCCalPad * kryptonRMS      = (AliTPCCalPad*)f.Get("spectrRMS");
   AliTPCCalPad * kryptonMeanG    = (AliTPCCalPad*)f.Get("fitMean");
   AliTPCCalPad * kryptonRMSG     = (AliTPCCalPad*)f.Get("fitRMS");
   AliTPCCalPad * kryptonNormChi2 = (AliTPCCalPad*)f.Get("fitNormChi2");
   AliTPCCalPad * kryptonEntries  = (AliTPCCalPad*)f.Get("entries");
   AliTPCCalPad * ceqIn           = (AliTPCCalPad*)fce.Get("qIn");
   AliTPCCalPad * ceqF1             = (AliTPCCalPad*)fce.Get("qF1");
   AliTPCCalPad * ceqF2             = (AliTPCCalPad*)fce.Get("qF2");


   preprocesor->AddComponent(kryptonMean->Clone());
   preprocesor->AddComponent(kryptonRMS->Clone());
   preprocesor->AddComponent(kryptonMeanG->Clone());
   preprocesor->AddComponent(kryptonRMSG->Clone());
   preprocesor->AddComponent(kryptonNormChi2->Clone());
   preprocesor->AddComponent(kryptonEntries->Clone());
   //
   preprocesor->AddComponent(ceqIn->Clone());
   preprocesor->AddComponent(ceqF1->Clone());
   preprocesor->AddComponent(ceqF2->Clone());

   preprocesor->DumpToFile("gainTree.root");
   //
}

AliTPCCalibViewerGUI*viewer =0;
TTree * tree =0;

void LoadViewer(){
  //
  // Load calib Viewer
  //
  TObjArray * array = AliTPCCalibViewerGUI::ShowGUI("gainTree.root");
  AliTPCCalibViewerGUI* viewer = (AliTPCCalibViewerGUI*)array->At(0);
  makePad = viewer->GetViewer();
  tree = viewer->GetViewer()->GetTree();

  tree->SetAlias("krAccept0","abs(fitRMS.fElements/fitMean.fElements-0.06)<0.04");
  tree->SetAlias("krAccept1","abs(fitRMS.fElements)>30");
  tree->SetAlias("yedge","tan(10*pi/180.)*lx.fElements");
  tree->SetAlias("ceAccept1","abs(qIn.fElements/qIn_Median.fElements-1.5)<1.4&&qIn.fElements>3&&qIn_Median.fElements>3");

}

void Fit(){
  TF1 f1("f1","[0]*exp(-[1]*x)+[2]");
  f1.SetParameters(1,1,0.2);
  tree->Draw("1-qIn.fElements/qF1.fElements:yedge-abs(ly.fElements)>>his(50,1.5,5,100,-0.5,1.5)","ceAccept1&&sector>36"); 
  his->FitSlicesY();
  his_1->Fit(&f1);


}
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
1b0bf200	191
	192
	193	void MakeCalibTree(char * inputKr="calibKr.root", char * inputCE ="fitCE.root", char * inputPulser=0){
	194	//
	195	//
	196	//
	197	AliTPCPreprocessorOnline * preprocesor = new AliTPCPreprocessorOnline;
	198	TFile f(inputKr);
	199	TFile fce(inputCE);
	200	AliTPCCalPad * kryptonMean = (AliTPCCalPad*)f.Get("spectrMean");
	201	AliTPCCalPad * kryptonRMS = (AliTPCCalPad*)f.Get("spectrRMS");
	202	AliTPCCalPad * kryptonMeanG = (AliTPCCalPad*)f.Get("fitMean");
	203	AliTPCCalPad * kryptonRMSG = (AliTPCCalPad*)f.Get("fitRMS");
	204	AliTPCCalPad * kryptonNormChi2 = (AliTPCCalPad*)f.Get("fitNormChi2");
	205	AliTPCCalPad * kryptonEntries = (AliTPCCalPad*)f.Get("entries");
	206	AliTPCCalPad * ceqIn = (AliTPCCalPad*)fce.Get("qIn");
	207	AliTPCCalPad * ceqF1 = (AliTPCCalPad*)fce.Get("qF1");
	208	AliTPCCalPad * ceqF2 = (AliTPCCalPad*)fce.Get("qF2");
	209
	210
	211
	212	preprocesor->AddComponent(kryptonMean->Clone());
	213	preprocesor->AddComponent(kryptonRMS->Clone());
	214	preprocesor->AddComponent(kryptonMeanG->Clone());
	215	preprocesor->AddComponent(kryptonRMSG->Clone());
	216	preprocesor->AddComponent(kryptonNormChi2->Clone());
	217	preprocesor->AddComponent(kryptonEntries->Clone());
	218	//
	219	preprocesor->AddComponent(ceqIn->Clone());
	220	preprocesor->AddComponent(ceqF1->Clone());
	221	preprocesor->AddComponent(ceqF2->Clone());
	222
	223	preprocesor->DumpToFile("gainTree.root");
	224	//
	225	}
	226
	227	AliTPCCalibViewerGUI*viewer =0;
	228	TTree * tree =0;
	229
	230	void LoadViewer(){
	231	//
	232	// Load calib Viewer
	233	//
	234	TObjArray * array = AliTPCCalibViewerGUI::ShowGUI("gainTree.root");
	235	AliTPCCalibViewerGUI* viewer = (AliTPCCalibViewerGUI*)array->At(0);
	236	makePad = viewer->GetViewer();
	237	tree = viewer->GetViewer()->GetTree();
	238
	239	tree->SetAlias("krAccept0","abs(fitRMS.fElements/fitMean.fElements-0.06)<0.04");
	240	tree->SetAlias("krAccept1","abs(fitRMS.fElements)>30");
	241	tree->SetAlias("yedge","tan(10pi/180.)lx.fElements");
	242	tree->SetAlias("ceAccept1","abs(qIn.fElements/qIn_Median.fElements-1.5)<1.4&&qIn.fElements>3&&qIn_Median.fElements>3");
	243
	244	}
	245
	246	void Fit(){
	247	TF1 f1("f1","[0]exp(-[1]x)+[2]");
	248	f1.SetParameters(1,1,0.2);
	249	tree->Draw("1-qIn.fElements/qF1.fElements:yedge-abs(ly.fElements)>>his(50,1.5,5,100,-0.5,1.5)","ceAccept1&&sector>36");
	250	his->FitSlicesY();
	251	his_1->Fit(&f1);
	252
	253
	254	}