TPC/CalibMacros/CreateGainMap.C

   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
 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(10*pi/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 }