Modification for DA (Christian Lippmann)
[u/mrichter/AliRoot.git] / TPC / AliTPCCalibPedestal.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16
17 /* $Id$ */
18
19
20 //Root includes
21 #include <TH1F.h>
22 #include <TH2F.h>
23 #include <TString.h>
24 #include <TMath.h>
25 #include <TF1.h>
26 #include <TRandom.h>
27 #include <TDirectory.h>
28 #include <TFile.h>
29 //AliRoot includes
30 #include "AliRawReader.h"
31 #include "AliRawReaderRoot.h"
32 #include "AliRawReaderDate.h"
33 #include "AliTPCRawStream.h"
34 #include "AliTPCCalROC.h"
35 #include "AliTPCROC.h"
36 #include "AliMathBase.h"
37 #include "TTreeStream.h"
38
39 //date
40 #include "event.h"
41
42 //header file
43 #include "AliTPCCalibPedestal.h"
44
45
46 ///////////////////////////////////////////////////////////////////////////////////////
47 //          Implementation of the TPC pedestal and noise calibration
48 //
49 //   Origin: Jens Wiechula, Marian Ivanov   J.Wiechula@gsi.de, Marian.Ivanov@cern.ch
50 // 
51 // 
52 // *************************************************************************************
53 // *                                Class Description                                  *
54 // *************************************************************************************
55 //
56 // Working principle:
57 // ------------------
58 // Raw pedestal data is processed by calling one of the ProcessEvent(...) functions
59 // (see below). These in the end call the Update(...) function, where the data is filled
60 // into histograms.
61 //
62 // For each ROC one TH2F histo (ROC channel vs. ADC channel) is created when
63 // it is filled for the first time (GetHistoPedestal(ROC,kTRUE)). All histos are stored in the
64 // TObjArray fHistoPedestalArray.
65 //
66 // For a fast filling of the histogram the corresponding bin number of the channel and ADC channel
67 // is computed by hand and the histogram array is accessed directly via its pointer.
68 // ATTENTION: Doing so the the entry counter of the histogram is not increased
69 //            this means that e.g. the colz draw option gives an empty plot unless
70 //          calling 'histo->SetEntries(1)' before drawing.
71 //
72 // After accumulating the desired statistics the Analyse() function has to be called.
73 // Whithin this function the pedestal and noise values are calculated for each pad, using
74 // the fast gaus fit function  AliMathBase::FitGaus(...), and the calibration
75 // storage classes (AliTPCCalROC) are filled for each ROC.
76 // The calibration information is stored in the TObjArrays fCalRocArrayPedestal and fCalRocArrayRMS;
77 //
78 //
79 //
80 // User interface for filling data:
81 // --------------------------------
82 //
83 // To Fill information one of the following functions can be used:
84 //
85 // Bool_t ProcessEvent(eventHeaderStruct *event);
86 //   - process Date event
87 //   - use AliTPCRawReaderDate and call ProcessEvent(AliRawReader *rawReader)
88 //
89 // Bool_t ProcessEvent(AliRawReader *rawReader);
90 //  - process AliRawReader event
91 //   - use AliTPCRawStream to loop over data and call ProcessEvent(AliTPCRawStream *rawStream)
92 //
93 // Bool_t ProcessEvent(AliTPCRawStream *rawStream);
94 //   - process event from AliTPCRawStream
95 //   - call Update function for signal filling
96 //
97 // Int_t Update(const Int_t isector, const Int_t iRow, const Int_t
98 //              iPad,  const Int_t iTimeBin, const Float_t signal);
99 //   - directly  fill signal information (sector, row, pad, time bin, pad)
100 //     to the reference histograms
101 //
102 // It is also possible to merge two independently taken calibrations using the function
103 //
104 // void Merge(AliTPCCalibPedestal *ped)
105 //   - copy histograms in 'ped' if the do not exist in this instance
106 //   - Add histograms in 'ped' to the histograms in this instance if the allready exist
107 //   - After merging call Analyse again!
108 //
109 //
110 //
111 // -- example: filling data using root raw data:
112 // void fillPedestal(Char_t *filename)
113 // {
114 //    rawReader = new AliRawReaderRoot(fileName);
115 //    if ( !rawReader ) return;
116 //    AliTPCCalibPedestal *calib = new AliTPCCalibPedestal;
117 //    while (rawReader->NextEvent()){
118 //      calib->ProcessEvent(rawReader);
119 //    }
120 //    calib->Analyse();
121 //    calib->DumpToFile("PedestalData.root");
122 //    delete rawReader;
123 //    delete calib;
124 // }
125 //
126 //
127 // What kind of information is stored and how to retrieve them:
128 // ------------------------------------------------------------
129 //
130 // - Accessing the 'Reference Histograms' (pedestal distribution histograms):
131 //
132 // TH2F *GetHistoPedestal(Int_t sector);
133 //
134 // - Accessing the calibration storage objects:
135 //
136 // AliTPCCalROC *GetCalRocPedestal(Int_t sector);  - for the pedestal values
137 // AliTPCCalROC *GetCalRocNoise(Int_t sector);     - for the Noise values
138 //
139 // example for visualisation:
140 // if the file "PedestalData.root" was created using the above example one could do the following:
141 //
142 // TFile filePedestal("PedestalData.root")
143 // AliTPCCalibPedestal *ped = (AliTPCCalibPedestal*)filePedestal->Get("AliTPCCalibPedestal");
144 // ped->GetCalRocPedestal(0)->Draw("colz");
145 // ped->GetCalRocRMS(0)->Draw("colz");
146 //
147 // or use the AliTPCCalPad functionality:
148 // AliTPCCalPad padPedestal(ped->GetCalPadPedestal());
149 // AliTPCCalPad padNoise(ped->GetCalPadRMS());
150 // padPedestal->MakeHisto2D()->Draw("colz");  //Draw A-Side Pedestal Information
151 // padNoise->MakeHisto2D()->Draw("colz");  //Draw A-Side Noise Information
152 //
153
154 /*
155  example: fill pedestal with gausschen noise
156  AliTPCCalibPedestal ped;
157  ped.TestEvent();
158  ped.Analyse();
159  //Draw output;
160  TCanvas* c1 = new TCanvas;
161  c1->Divide(1,2);
162  c1->cd(1);
163  ped.GetHistoPedestal(0)->SetEntries(1); //needed in order for colz to work, reason: fast filling does not increase the entries counter
164  ped.GetHistoPedestal(0)->Draw("colz");
165  c1->cd(2);
166  ped.GetHistoPedestal(36)->SetEntries(1); //needed in order for colz to work, reason: fast filling does not increase the entries counter
167  ped.GetHistoPedestal(36)->Draw("colz");
168  TCanvas* c2 = new TCanvas;
169  c2->Divide(2,2);
170  c2->cd(1);
171  ped.GetCalRocPedestal(0)->Draw("colz");
172  c2->cd(2);
173  ped.GetCalRocRMS(0)->Draw("colz");
174  c2->cd(3);
175  ped.GetCalRocPedestal(36)->Draw("colz");
176  c2->cd(4);
177  ped.GetCalRocRMS(36)->Draw("colz");
178 */
179
180 //
181 // Time dependent pedestals:
182 //
183 // If wished there is the possibility to calculate for each channel and time bin
184 // the mean pedestal [pedestals(t)]. This is done by
185 //
186 // 1) setting SetTimeAnalysis(kTRUE),
187 // 2) processing the data by looping over the events using ProcessEvent(..)
188 // 3) calling the Analyse() and AnalyseTime(nevents) functions (providing nevents)
189 // 4) getting the pedestals(t) using   TArrayF **timePed = calibPedestal.GetTimePedestals();
190 // 5) looking at values using   timePed[row][pad].At(timebin)
191 //
192 // This functionality is intended to be used on an LDC bu the detector algorithm
193 // (TPCPEDESTALda) to generate a data set used for configuration of the pattern
194 // memory for baseline subtraction in the ALTROs. Later the information should also
195 // be stored as reference data.
196 //
197
198
199 ClassImp(AliTPCCalibPedestal)
200
201 AliTPCCalibPedestal::AliTPCCalibPedestal() : /*FOLD00*/
202   TObject(),
203   fFirstTimeBin(60),
204   fLastTimeBin(1000),
205   fAdcMin(1),
206   fAdcMax(100),
207   fOldRCUformat(kTRUE),
208   fTimeAnalysis(kFALSE),
209   fROC(AliTPCROC::Instance()),
210   fMapping(NULL),
211   fCalRocArrayPedestal(72),
212   fCalRocArrayRMS(72),
213   fHistoPedestalArray(72),
214   fTimeSignal(NULL)
215 {
216   //
217   // default constructor
218   //
219 }
220
221
222 //_____________________________________________________________________
223 AliTPCCalibPedestal::AliTPCCalibPedestal(const AliTPCCalibPedestal &ped) : /*FOLD00*/
224   TObject(ped),
225   fFirstTimeBin(ped.GetFirstTimeBin()),
226   fLastTimeBin(ped.GetLastTimeBin()),
227   fAdcMin(ped.GetAdcMin()),
228   fAdcMax(ped.GetAdcMax()),
229   fOldRCUformat(ped.fOldRCUformat),
230   fTimeAnalysis(ped.fTimeAnalysis),
231   fROC(AliTPCROC::Instance()),
232   fMapping(NULL),
233   fCalRocArrayPedestal(72),
234   fCalRocArrayRMS(72),
235   fHistoPedestalArray(72),
236   fTimeSignal(ped.fTimeSignal)
237 {
238   //
239   // copy constructor
240   //
241   for (Int_t iSec = 0; iSec < 72; ++iSec){
242     const AliTPCCalROC *calPed = (AliTPCCalROC*)ped.fCalRocArrayPedestal.UncheckedAt(iSec);
243     const AliTPCCalROC *calRMS = (AliTPCCalROC*)ped.fCalRocArrayRMS.UncheckedAt(iSec);
244     const TH2F         *hPed   = (TH2F*)ped.fHistoPedestalArray.UncheckedAt(iSec);
245     
246     if ( calPed != 0x0 ) fCalRocArrayPedestal.AddAt(new AliTPCCalROC(*calPed), iSec);
247     if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
248     
249     if ( hPed != 0x0 ){
250       TH2F *hNew = new TH2F(*hPed);
251       hNew->SetDirectory(0);
252       fHistoPedestalArray.AddAt(hNew,iSec);
253     }
254   }
255 }
256
257
258 //_____________________________________________________________________
259 AliTPCCalibPedestal& AliTPCCalibPedestal::operator = (const  AliTPCCalibPedestal &source)
260 {
261   //
262   // assignment operator
263   //
264   if (&source == this) return *this;
265   new (this) AliTPCCalibPedestal(source);
266
267   return *this;
268 }
269
270
271 //_____________________________________________________________________
272 AliTPCCalibPedestal::~AliTPCCalibPedestal() /*FOLD00*/
273 {
274   //
275   // destructor
276   //
277
278   fCalRocArrayPedestal.Delete();
279   fCalRocArrayRMS.Delete();
280   fHistoPedestalArray.Delete();
281
282   if ( fTimeSignal ) {
283     for (Int_t i = 0; i < 159; i++) {
284       delete [] fTimeSignal[i];
285       fTimeSignal[i] = 0;
286     }
287     delete [] fTimeSignal;
288     fTimeSignal = 0;
289   }
290
291   // do not delete fMapping, because we do not own it.
292
293 }
294
295
296 //_____________________________________________________________________
297 void AliTPCCalibPedestal::SetTimeAnalysis(Bool_t time)
298 {
299   //
300   // Use time dependent analysis: Pedestals are analysed as a function
301   // of the drift time. There is one mean value generated for each time
302   // bin and each channel. It can be used as reference data and for
303   // configuration of the ALTRO pattern memory for baseline subtraction.
304   //
305   // ATTENTION: Use only on LDC in TPCPEDESTALda! On a LDC we get data
306   // only from one sector. For the full TPC we would need a lot of
307   // memory (36*159*140*1024*4bytes = 3.3GB)!
308   //
309
310   fTimeAnalysis = time;
311
312   if ( !fTimeAnalysis ) return;
313
314   // prepare array for one sector (159*140*1024*4bytes = 92MB):
315   fTimeSignal = new TArrayF*[159];
316   for (Int_t i = 0; i < 159; i++) {  // padrows
317     fTimeSignal[i] = new TArrayF[140];
318     for (Int_t j = 0; j < 140; j++) {  // pads per row
319       fTimeSignal[i][j].Set(1024);
320       for (Int_t k = 0; k < 1024; k++) {  // time bins per pad
321         fTimeSignal[i][j].AddAt(0., k);
322       }
323     }      
324   }
325 }
326
327
328 //_____________________________________________________________________
329 Int_t AliTPCCalibPedestal::Update(const Int_t icsector, /*FOLD00*/
330                                   const Int_t icRow,
331                                   const Int_t icPad,
332                                   const Int_t icTimeBin,
333                                   const Float_t csignal)
334 {
335   //
336   // Signal filling method
337   //
338
339   // Time dependent pedestals
340   if ( fTimeAnalysis ) {
341     if ( icsector < 36 ) // IROC
342       fTimeSignal[icRow][icPad].AddAt(fTimeSignal[icRow][icPad].At(icTimeBin)+csignal, icTimeBin);
343     else 
344       fTimeSignal[icRow+63][icPad].AddAt(fTimeSignal[icRow+63][icPad].At(icTimeBin)+csignal, icTimeBin);
345   }
346   //return if we are out of the specified time bin or adc range
347   if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
348   if ( ((Int_t)csignal>fAdcMax) || ((Int_t)csignal<fAdcMin)  ) return 0;
349
350   Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
351
352   // fast filling method
353   // Attention: the entry counter of the histogram is not increased
354   //            this means that e.g. the colz draw option gives an empty plot
355   Int_t bin = (iChannel+1)*(fAdcMax-fAdcMin+2)+((Int_t)csignal-fAdcMin+1);
356
357   GetHistoPedestal(icsector,kTRUE)->GetArray()[bin]++;
358
359   return 0;
360 }
361
362
363 //_____________________________________________________________________
364 Bool_t AliTPCCalibPedestal::ProcessEvent(AliTPCRawStream *rawStream)
365 {
366   //
367   // Event Processing loop - AliTPCRawStream
368   //
369
370   rawStream->SetOldRCUFormat(fOldRCUformat);
371
372   Bool_t withInput = kFALSE;
373
374   while (rawStream->Next()) {
375
376     Int_t iSector  = rawStream->GetSector();      //  current ROC
377     Int_t iRow     = rawStream->GetRow();         //  current row
378     Int_t iPad     = rawStream->GetPad();         //  current pad
379     Int_t iTimeBin = rawStream->GetTime();        //  current time bin
380     Float_t signal = rawStream->GetSignal();      //  current ADC signal
381     
382     Update(iSector,iRow,iPad,iTimeBin,signal);
383     withInput = kTRUE;
384   }
385
386   return withInput;
387 }
388
389
390 //_____________________________________________________________________
391 Bool_t AliTPCCalibPedestal::ProcessEvent(AliRawReader *rawReader)
392 {
393   //
394   //  Event processing loop - AliRawReader
395   //
396
397   // if fMapping is NULL the rawstream will crate its own mapping
398   AliTPCRawStream rawStream(rawReader, (AliAltroMapping**)fMapping);
399   rawReader->Select("TPC");
400   return ProcessEvent(&rawStream);
401 }
402
403
404 //_____________________________________________________________________
405 Bool_t AliTPCCalibPedestal::ProcessEvent(eventHeaderStruct *event)
406 {
407   //
408   //  process date event
409   //
410
411   AliRawReader *rawReader = new AliRawReaderDate((void*)event);
412   Bool_t result=ProcessEvent(rawReader);
413   delete rawReader;
414   return result;
415 }
416
417
418 //_____________________________________________________________________
419 Bool_t AliTPCCalibPedestal::TestEvent() /*FOLD00*/
420 {
421   //
422   //  Test event loop
423   // fill one oroc and one iroc with random gaus
424   //
425
426     gRandom->SetSeed(0);
427
428     for (UInt_t iSec=0; iSec<72; ++iSec){
429         if (iSec%36>0) continue;
430         for (UInt_t iRow=0; iRow < fROC->GetNRows(iSec); ++iRow){
431             for (UInt_t iPad=0; iPad < fROC->GetNPads(iSec,iRow); ++iPad){
432                 for (UInt_t iTimeBin=0; iTimeBin<1024; ++iTimeBin){
433                     Float_t signal=(Int_t)(iRow+3+gRandom->Gaus(0,.7));
434                     if ( signal>0 )Update(iSec,iRow,iPad,iTimeBin,signal);
435                 }
436             }
437         }
438     }
439     return kTRUE;
440 }
441
442
443 //_____________________________________________________________________
444 TH2F* AliTPCCalibPedestal::GetHisto(Int_t sector, TObjArray *arr, /*FOLD00*/
445                                     Int_t nbinsY, Float_t ymin, Float_t ymax,
446                                     Char_t *type, Bool_t force)
447 {
448     //
449     // return pointer to Q histogram
450     // if force is true create a new histogram if it doesn't exist allready
451     //
452     if ( !force || arr->UncheckedAt(sector) )
453         return (TH2F*)arr->UncheckedAt(sector);
454
455     // if we are forced and histogram doesn't yes exist create it
456     Char_t name[255], title[255];
457
458     sprintf(name,"hCalib%s%.2d",type,sector);
459     sprintf(title,"%s calibration histogram sector %.2d;ADC channel;Channel (pad)",type,sector);
460
461     // new histogram with Q calib information. One value for each pad!
462     TH2F* hist = new TH2F(name,title,
463                           nbinsY, ymin, ymax,
464                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector)
465                          );
466     hist->SetDirectory(0);
467     arr->AddAt(hist,sector);
468     return hist;
469 }
470
471
472 //_____________________________________________________________________
473 TH2F* AliTPCCalibPedestal::GetHistoPedestal(Int_t sector, Bool_t force) /*FOLD00*/
474 {
475     //
476     // return pointer to T0 histogram
477     // if force is true create a new histogram if it doesn't exist allready
478     //
479     TObjArray *arr = &fHistoPedestalArray;
480     return GetHisto(sector, arr, fAdcMax-fAdcMin, fAdcMin, fAdcMax, "Pedestal", force);
481 }
482
483
484 //_____________________________________________________________________
485 AliTPCCalROC* AliTPCCalibPedestal::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) /*FOLD00*/
486 {
487     //
488     // return pointer to ROC Calibration
489     // if force is true create a new histogram if it doesn't exist allready
490     //
491     if ( !force || arr->UncheckedAt(sector) )
492         return (AliTPCCalROC*)arr->UncheckedAt(sector);
493
494     // if we are forced and the histogram doesn't yet exist create it
495
496     // new AliTPCCalROC for T0 information. One value for each pad!
497     AliTPCCalROC *croc = new AliTPCCalROC(sector);
498     arr->AddAt(croc,sector);
499     return croc;
500 }
501
502
503 //_____________________________________________________________________
504 AliTPCCalROC* AliTPCCalibPedestal::GetCalRocPedestal(Int_t sector, Bool_t force) /*FOLD00*/
505 {
506     //
507     // return pointer to Carge ROC Calibration
508     // if force is true create a new histogram if it doesn't exist allready
509     //
510     TObjArray *arr = &fCalRocArrayPedestal;
511     return GetCalRoc(sector, arr, force);
512 }
513
514
515 //_____________________________________________________________________
516 AliTPCCalROC* AliTPCCalibPedestal::GetCalRocRMS(Int_t sector, Bool_t force) /*FOLD00*/
517 {
518     //
519     // return pointer to signal width ROC Calibration
520     // if force is true create a new histogram if it doesn't exist allready
521     //
522     TObjArray *arr = &fCalRocArrayRMS;
523     return GetCalRoc(sector, arr, force);
524 }
525
526
527 //_____________________________________________________________________
528 void AliTPCCalibPedestal::Merge(AliTPCCalibPedestal *ped)
529 {
530   //
531   //  Merge reference histograms of sig to the current AliTPCCalibSignal
532   //
533
534   // merge histograms
535   for (Int_t iSec=0; iSec<72; ++iSec){
536     TH2F *hRefPedMerge   = ped->GetHistoPedestal(iSec);
537
538     if ( hRefPedMerge ){
539       TDirectory *dir = hRefPedMerge->GetDirectory(); hRefPedMerge->SetDirectory(0);
540       TH2F *hRefPed   = GetHistoPedestal(iSec);
541       if ( hRefPed ) hRefPed->Add(hRefPedMerge);
542       else {
543         TH2F *hist = new TH2F(*hRefPedMerge);
544         hist->SetDirectory(0);
545         fHistoPedestalArray.AddAt(hist, iSec);
546       }
547       hRefPedMerge->SetDirectory(dir);
548     }
549   }
550
551   // merge array
552   // ...
553
554 }
555
556
557 //_____________________________________________________________________
558 void AliTPCCalibPedestal::Analyse() /*FOLD00*/
559 {
560   //
561   //  Calculate calibration constants
562   //
563
564   Int_t nbinsAdc = fAdcMax-fAdcMin;
565
566   TVectorD param(3);
567   TMatrixD dummy(3,3);
568
569   Float_t *array_hP=0;  
570
571   for (Int_t iSec=0; iSec<72; ++iSec){
572     TH2F *hP = GetHistoPedestal(iSec);
573     if ( !hP ) continue;
574
575     AliTPCCalROC *rocPedestal = GetCalRocPedestal(iSec,kTRUE);
576     AliTPCCalROC *rocRMS      = GetCalRocRMS(iSec,kTRUE);
577
578     array_hP = hP->GetArray();
579     UInt_t nChannels = fROC->GetNChannels(iSec);
580
581     for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
582       Int_t offset = (nbinsAdc+2)*(iChannel+1)+1;
583       Double_t ret = AliMathBase::FitGaus(array_hP+offset,nbinsAdc,fAdcMin,fAdcMax,&param,&dummy);
584       // if the fitting failed set noise and pedestal to 0
585       if ( ret == -4 ) {
586         param[1]=0;
587         param[2]=0;
588       }
589       rocPedestal->SetValue(iChannel,param[1]);
590       rocRMS->SetValue(iChannel,param[2]);
591     }
592   }
593 }
594
595
596 //_____________________________________________________________________
597 void AliTPCCalibPedestal::AnalyseTime(Int_t nevents)
598 {
599   //
600   // Calculate for each channel and time bin the mean pedestal. This
601   // is used on LDC by TPCPEDESTALda to generate data used for configuration
602   // of the pattern memory for baseline subtraction in the ALTROs.
603   //
604
605   if ( nevents <= 0 ) return;
606   if ( fTimeAnalysis ) {
607     for (Int_t i = 0; i < 159; i++) {  // padrows
608       for (Int_t j = 0; j < 140; j++) {  // pads per row
609         for (Int_t k = 0; k < 1024; k++) {  // time bins per pad
610           fTimeSignal[i][j].AddAt(fTimeSignal[i][j].At(k)/(Float_t)nevents, k);
611         }
612       }
613     }
614   }
615 }
616
617
618 //_____________________________________________________________________
619 void AliTPCCalibPedestal::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append) /*FOLD00*/
620 {
621   //
622   //  Write class to file
623   //
624
625   TString sDir(dir);
626   TString option;
627
628   if ( append )
629     option = "update";
630   else
631     option = "recreate";
632
633   TDirectory *backup = gDirectory;
634   TFile f(filename,option.Data());
635   f.cd();
636   if ( !sDir.IsNull() ){
637     f.mkdir(sDir.Data());
638     f.cd(sDir);
639   }
640   this->Write();
641   f.Close();
642
643   if ( backup ) backup->cd();
644 }