141557110142b9ec2e9d82cdb074f5addefa7309
[u/mrichter/AliRoot.git] / TRD / AliTRDCalibPadStatus.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 /* $Id$ */
17
18 ////////////////////////////////////////////////////////////////////////////
19 //                                                                        //
20 // Example: fill pedestal with Gaussian noise                             //
21 //                                                                        //
22 //  AliTRDCalibPadStatus ped;                                             //
23 //  ped.TestEvent(numberofevent);                                         //
24 //                                                                        //
25 //  // Method without histo                                               //
26 //  ped.Analyse();                                                        //
27 //                                                                        //
28 //  // Create the histo of the AliTRDCalROC                               //
29 //  TH2F * histo2dm = ped.GetCalRocMean(0,kFALSE)->MakeHisto2D();         //
30 //  histo2dm->Scale(10.0);                                                //
31 //  TH1F * histo1dm = ped.GetCalRocMean(0,kFALSE)->MakeHisto1D();         //
32 //  histo1dm->Scale(10.0);                                                //
33 //  TH2F * histo2ds = ped.GetCalRocSquares(0,kFALSE)->MakeHisto2D();      //
34 //  histo2ds->Scale(10.0);                                                //
35 //  TH1F * histo1ds = ped.GetCalRocSquares(0,kFALSE)->MakeHisto1D();      //
36 //  histo1ds->Scale(10.0)                                                 //
37 //                                                                        //
38 //  // Draw output                                                        //
39 //  TCanvas* c1 = new TCanvas;                                            //
40 //  c1->Divide(2,2);                                                      //
41 //  c1->cd(1);                                                            //
42 //  histo2dm->Draw("colz");                                               //
43 //  c1->cd(2);                                                            //
44 //  histo1dm->Draw();                                                     //
45 //  c1->cd(3);                                                            //
46 //  histo2ds->Draw("colz");                                               //
47 //  c1->cd(4);                                                            //
48 //  histo1ds->Draw();                                                     //
49 //                                                                        //
50 //  // Method with histo                                                  //
51 //  ped.AnalyseHisto();                                                   //
52 //                                                                        //
53 //  // Take the histo                                                     //
54 //  TH1F *histo = ped.GetHisto(31);                                       //
55 //  histo->SetEntries(1);                                                 //
56 //  histo->Draw();                                                        //
57 //
58 // Authors:
59 //   R. Bailhache (R.Bailhache@gsi.de, rbailhache@ikf.uni-frankfurt.de)
60 //   J. Book (jbook@ikf.uni-frankfurt.de)
61 //                                                                                                    //
62 ////////////////////////////////////////////////////////////////////////////
63
64
65 //Root includes
66 #include <TObjArray.h>
67 #include <TH2F.h>
68 #include <TString.h>
69 #include <TMath.h>
70 #include <TRandom.h>
71
72 //#include <TRandom.h>
73 #include <TDirectory.h>
74 #include <TFile.h>
75
76 //AliRoot includes
77 #include <AliMathBase.h>
78 #include "AliRawReader.h"
79 #include "AliRawReaderRoot.h"
80 #include "AliRawReaderDate.h"
81
82 //header file
83 #include "AliLog.h"
84 #include "AliTRDCalibPadStatus.h"
85 #include "AliTRDrawStreamBase.h"
86 #include "AliTRDgeometry.h"
87 #include "AliTRDCommonParam.h"
88 #include "./Cal/AliTRDCalROC.h"
89 #include "./Cal/AliTRDCalPadStatus.h"
90 #include "./Cal/AliTRDCalDet.h"
91 #include "./Cal/AliTRDCalPad.h"
92 #include "./Cal/AliTRDCalSingleChamberStatus.h"
93
94 #include "AliTRDrawFastStream.h"
95 #include "AliTRDdigitsManager.h"
96 #include "AliTRDdigitsParam.h"
97 #include "AliTRDSignalIndex.h"
98 #include "AliTRDarraySignal.h"
99 #include "AliTRDarrayADC.h"
100 #include "AliTRDfeeParam.h"
101
102 #ifdef ALI_DATE
103 #include "event.h"
104 #endif
105
106 ClassImp(AliTRDCalibPadStatus) /*FOLD00*/
107
108 //_____________________________________________________________________
109 AliTRDCalibPadStatus::AliTRDCalibPadStatus() : /*FOLD00*/
110   TObject(),
111   fGeo(0),
112   fAdcMin(0),
113   fAdcMax(21),
114   fDetector(-1),
115   fNumberOfTimeBins(0),
116   fCalRocArrayMean(540),
117   fCalRocArrayRMS(540),
118   fCalRocArrayMeand(540),
119   fCalRocArrayRMSd(540),
120   fHistoArray(540)
121 {
122     //
123     // default constructor
124     //
125
126   fGeo = new AliTRDgeometry();
127
128 }
129
130 //_____________________________________________________________________
131 AliTRDCalibPadStatus::AliTRDCalibPadStatus(const AliTRDCalibPadStatus &ped) : /*FOLD00*/
132   TObject(ped),
133   fGeo(0),
134   fAdcMin(ped.GetAdcMin()),
135   fAdcMax(ped.GetAdcMax()),
136   fDetector(ped.fDetector),
137   fNumberOfTimeBins(ped.fNumberOfTimeBins),
138   fCalRocArrayMean(540),
139   fCalRocArrayRMS(540),
140   fCalRocArrayMeand(540),
141   fCalRocArrayRMSd(540),
142   fHistoArray(540)
143 {
144     //
145     // copy constructor
146     //
147     for (Int_t idet = 0; idet < 540; idet++){
148         const AliTRDCalROC *calRocMean  = (AliTRDCalROC*)ped.fCalRocArrayMean.UncheckedAt(idet);
149         const AliTRDCalROC *calRocRMS   = (AliTRDCalROC*)ped.fCalRocArrayRMS.UncheckedAt(idet);
150         const AliTRDCalROC *calRocMeand = (AliTRDCalROC*)ped.fCalRocArrayMeand.UncheckedAt(idet);
151         const AliTRDCalROC *calRocRMSd  = (AliTRDCalROC*)ped.fCalRocArrayRMSd.UncheckedAt(idet);
152         const TH2F         *hped        = (TH2F*)ped.fHistoArray.UncheckedAt(idet);
153     
154         if ( calRocMean != 0x0 ) fCalRocArrayMean.AddAt(new AliTRDCalROC(*calRocMean), idet);
155         if ( calRocRMS != 0x0 )  fCalRocArrayRMS.AddAt(new AliTRDCalROC(*calRocRMS), idet);
156
157         if ( calRocMeand != 0x0 ) fCalRocArrayMeand.AddAt(new AliTRDCalROC(*calRocMeand), idet);
158         if ( calRocRMSd != 0x0 )  fCalRocArrayRMSd.AddAt(new AliTRDCalROC(*calRocRMSd), idet);
159
160         if ( hped != 0x0 ){
161           TH2F *hNew = new TH2F(*hped);
162           hNew->SetDirectory(0);
163           fHistoArray.AddAt(hNew,idet);
164         }
165         
166     }
167     if (fGeo) {
168       delete fGeo;
169     }
170     fGeo = new AliTRDgeometry();
171 }
172
173 //_____________________________________________________________________
174 AliTRDCalibPadStatus& AliTRDCalibPadStatus::operator = (const  AliTRDCalibPadStatus &source)
175 {
176   //
177   // assignment operator
178   //
179   if (&source == this) return *this;
180   new (this) AliTRDCalibPadStatus(source);
181
182   return *this;
183 }
184 //_____________________________________________________________________
185 AliTRDCalibPadStatus::~AliTRDCalibPadStatus() /*FOLD00*/
186 {
187   //
188   // destructor
189   //
190   fCalRocArrayMean.Delete();
191   fCalRocArrayRMS.Delete();
192   fCalRocArrayMeand.Delete();
193   fCalRocArrayRMSd.Delete();
194   fHistoArray.Delete();
195   if (fGeo) {
196     delete fGeo;
197   }
198 }
199 //_____________________________________________________________________
200 void AliTRDCalibPadStatus::Destroy()
201 {
202   //
203   // Destroy
204   //
205   fCalRocArrayMean.Delete();
206   fCalRocArrayRMS.Delete();
207   fCalRocArrayMeand.Delete();
208   fCalRocArrayRMSd.Delete();
209   fHistoArray.Delete();
210 }
211 //_____________________________________________________________________
212 Int_t AliTRDCalibPadStatus::UpdateHisto(const Int_t icdet, /*FOLD00*/
213                                         const Int_t icRow,
214                                         const Int_t icCol,
215                                         const Int_t csignal,
216                                         const Int_t crowMax,
217                                         const Int_t ccold,
218                                         const Int_t icMcm)
219 {
220   //
221   // Signal filling methode 
222   //
223   Int_t nbchannel = icRow+icCol*crowMax;
224
225   // now the case of double read channel
226   if(ccold > 0){
227     nbchannel = (((ccold-1)*8+ icMcm)*crowMax+icRow)+144*crowMax;
228     //printf("nbchannel %d, ccold %d, icMcm %d, crowMax %d, icRow %d\n",nbchannel,ccold,icMcm,crowMax,icRow);
229   }
230   
231   // fast filling methode.
232   // Attention: the entry counter of the histogram is not increased
233   //            this means that e.g. the colz draw option gives an empty plot
234   
235   Int_t bin = 0;
236   
237   if ( !(((Int_t)csignal>=fAdcMax ) || ((Int_t)csignal<fAdcMin)) )
238     bin = (nbchannel+1)*(fAdcMax-fAdcMin+2)+((Int_t)csignal-fAdcMin+1);
239   
240   //GetHisto(icdet,kTRUE)->Fill(csignal,nbchannel);
241   
242   GetHisto(icdet,kTRUE)->GetArray()[bin]++;
243   
244   return 0;
245 }
246 //_____________________________________________________________________
247 Int_t AliTRDCalibPadStatus::UpdateHisto2(const Int_t icdet, /*FOLD00*/
248                                          const Int_t icRow,
249                                          const Int_t icCol,
250                                          const Int_t csignal,
251                                          const Int_t crowMax,
252                                          const Int_t ccold,
253                                          const Int_t icMcm,
254                                          const Int_t icRob
255                                          )
256 {
257   //
258   // Signal filling methode 
259   //
260   Int_t nbchannel = icRow+icCol*crowMax;
261   Int_t mCm = icMcm%4;
262   Int_t rOb = icRob%2;
263
264   // now the case of double read channel
265   if(ccold > 0){
266     nbchannel = (((ccold-1)*8+ (mCm+rOb*4))*crowMax+icRow)+144*crowMax;
267     //printf("nbchannel %d, ccold %d, icMcm %d, crowMax %d, icRow %d\n",nbchannel,ccold,icMcm,crowMax,icRow);
268   }
269   
270   // fast filling methode.
271   // Attention: the entry counter of the histogram is not increased
272   //            this means that e.g. the colz draw option gives an empty plot
273   
274   Int_t bin = 0;
275   
276   if ( !(((Int_t)csignal>=fAdcMax ) || ((Int_t)csignal<fAdcMin)) )
277     bin = (nbchannel+1)*(fAdcMax-fAdcMin+2)+((Int_t)csignal-fAdcMin+1);
278
279   //GetHisto(icdet,kTRUE)->Fill(csignal,nbchannel);
280
281   GetHisto(icdet,kTRUE)->GetArray()[bin]++;
282   
283   return 0;
284 }
285 //_____________________________________________________________________
286 Int_t AliTRDCalibPadStatus::ProcessEvent(AliTRDrawStreamBase *rawStream, Bool_t nocheck)
287 {
288   //
289   // Event Processing loop - AliTRDRawStreamCosmic
290   // 0 time bin problem or zero suppression
291   // 1 no input
292   // 2 input
293   // 
294
295   //
296   // Raw version number: 
297   // [3,31] non zero suppressed
298   // 2,4 and [32,63] zero suppressed 
299   //
300
301   Int_t withInput = 1;
302
303   rawStream->SetSharedPadReadout(kTRUE);
304
305   if(!nocheck) {
306
307     // Check the raw version and if all have the same number of timebins. 
308
309     while (rawStream->Next()) {
310
311       Int_t rawversion = rawStream->GetRawVersion();                     //  current raw version
312       //printf("Raw version is %d\n",rawversion);
313
314       // Could eventually change, have to check with time    
315       if((rawversion < 3) || (rawversion > 31)) {
316         AliInfo(Form("this is not no-zero-suppressed data, the version is %d",rawversion));
317         return 0;
318       }
319       Int_t idetector  = rawStream->GetDet();                            //  current detector
320       Int_t iRow       = rawStream->GetRow();                            //  current row
321       Int_t iRowMax    = rawStream->GetMaxRow();                         //  current rowmax
322       Int_t iCol       = rawStream->GetCol();                            //  current col
323       Int_t iADC       = 21-rawStream->GetADC();                         //  current ADC
324       
325       // It goes in the opposite direction
326       Int_t col        = 0;
327       if(iADC == 1) col = 1;
328       else {
329         col = TMath::Max(0,(Int_t)(iADC-19));
330         if(col > 0) col++;
331       }
332       Int_t mcm        = (Int_t)(iCol/18);                               //  current group of 18 col pads
333       if(col > 1) mcm -= 1;      
334       if(col ==1) mcm += 1;
335
336       // printf to check
337       //Bool_t shared = rawStream->IsCurrentPadShared();                  
338       //printf("ADC %d, iCol %d, col %d, mcm %d, shared %d\n",iADC,iCol,col,mcm,(Int_t)shared);
339
340       // Take the signal
341       Int_t *signal    = rawStream->GetSignals();                        //  current ADC signal
342       Int_t nbtimebin  = rawStream->GetNumberOfTimeBins();               //  number of time bins read from data
343
344       if((fDetector != -1) && (nbtimebin != fNumberOfTimeBins)) {
345         AliInfo(Form("the number of time bins is %d, is different from the previous one %d",nbtimebin,fNumberOfTimeBins));
346         return 0;
347       }
348       fNumberOfTimeBins = nbtimebin;
349       fDetector         = idetector;      
350
351       for(Int_t k = 0; k < fNumberOfTimeBins; k++){
352         if(signal[k]>0 && iCol != -1) UpdateHisto(idetector,iRow,iCol,signal[k],iRowMax,col,mcm);
353       }
354       
355       withInput = 2;
356     }
357   }
358   else {
359
360     while (rawStream->Next()) {
361     
362       Int_t idetector  = rawStream->GetDet();                            //  current detector
363       Int_t iRow       = rawStream->GetRow();                            //  current row
364       Int_t iRowMax    = rawStream->GetMaxRow();                         //  current rowmax
365       Int_t iCol       = rawStream->GetCol();                            //  current col
366       Int_t iADC       = 21-rawStream->GetADC();                            //  current ADC
367
368       // It goes in the opposite direction      
369       Int_t col        = 0;
370       if(iADC == 1) col = 1;
371       else {
372         col = TMath::Max(0,(Int_t)(iADC-19));
373         if(col > 0) col++;
374       }
375       Int_t mcm        = (Int_t)(iCol/18);                               //  current group of 18 col pads
376       if(col > 1) mcm -= 1;      
377       if(col ==1) mcm += 1;
378
379       // Take the signal
380       Int_t *signal    = rawStream->GetSignals();                        //  current ADC signal
381       Int_t nbtimebin = rawStream->GetNumberOfTimeBins();               //  number of time bins read from data
382       
383       
384       //printf("det %d, row %d, signal[0] %d, signal[1] %d, signal [2] %d\n", idetector, iRow, signal[0], signal[1], signal[2]);
385  
386       for(Int_t k = 0; k < nbtimebin; k++){
387         if(signal[k]>0 && iCol != -1) {
388           UpdateHisto(idetector,iRow,iCol,signal[k],iRowMax,col,mcm);
389           //printf("Update with det %d, row %d, col %d, signal %d, rowmax %d, col %d, mcm %d\n",idetector,iRow,iCol,signal[n],iRowMax,col,mcm);
390         }
391       }
392       
393       withInput = 2;
394     }
395   }
396   
397   return withInput;
398 }
399 //_____________________________________________________________________
400 Int_t AliTRDCalibPadStatus::ProcessEvent(AliRawReader *rawReader, Bool_t nocheck)
401 {
402   //
403   //  Event processing loop - AliRawReader
404   //
405
406   Int_t result;
407   
408   rawReader->Select("TRD");
409   
410   AliTRDrawStreamBase *pstream = AliTRDrawStreamBase::GetRawStream(rawReader);
411  
412   result = ProcessEvent(pstream, nocheck);
413
414   delete pstream;
415
416   return result;
417 }
418
419 //_________________________________________________________________________
420 Int_t AliTRDCalibPadStatus::ProcessEvent(
421 #ifdef ALI_DATE
422                                           const eventHeaderStruct *event,
423                                           Bool_t nocheck
424 #else
425                                           const eventHeaderStruct* /*event*/,
426                                           Bool_t /*nocheck*/
427             
428 #endif 
429                                           )
430 {
431   //
432   //  process date event
433   //
434 #ifdef ALI_DATE
435     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
436     Bool_t result=ProcessEvent(rawReader, nocheck);
437     delete rawReader;
438     return result;
439 #else
440     Fatal("AliTRDCalibPadStatus", "this class was compiled without DATE");
441     return 0;
442 #endif
443
444 }
445
446 //_____________________________________________________________________
447 Int_t AliTRDCalibPadStatus::ProcessEvent2(AliRawReader *rawReader)
448 {
449   //
450   // Event Processing loop - AliTRDRawStreamCosmic
451   // 0 time bin problem or zero suppression
452   // 1 no input
453   // 2 input
454   // Raw version number: 
455   // [3,31] non zero suppressed
456   // 2,4 and [32,63] zero suppressed 
457   //
458   
459   Int_t withInput = 1;
460
461   AliTRDrawFastStream *rawStream = new AliTRDrawFastStream(rawReader);
462   rawStream->SetNoErrorWarning();
463   rawStream->SetSharedPadReadout(kTRUE);
464
465   AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager(kTRUE);
466   digitsManager->CreateArrays();
467   
468   AliTRDfeeParam *feeParam = AliTRDfeeParam::Instance();
469
470   Int_t det    = 0;
471   while ((det = rawStream->NextChamber(digitsManager, NULL, NULL)) >= 0) { //idetector
472     if (digitsManager->GetIndexes(det)->HasEntry()) {//QA
473       //        printf("there is ADC data on this chamber!\n");
474       
475        AliTRDarrayADC *digits = (AliTRDarrayADC *) digitsManager->GetDigits(det); //mod
476       if (digits->HasData()) { //array
477         
478         AliTRDSignalIndex   *indexes = digitsManager->GetIndexes(det);
479         if (indexes->IsAllocated() == kFALSE) {
480           AliError("Indexes do not exist!");
481           break;
482         }
483         Int_t iRow  = 0;
484         Int_t iCol  = 0;
485         indexes->ResetCounters();
486
487         while (indexes->NextRCIndex(iRow, iCol)) { //column,row
488           
489
490           AliTRDdigitsParam *digitParam = (AliTRDdigitsParam *)digitsManager->GetDigitsParam();
491           
492           Int_t mcm          = 0;     // MCM from AliTRDfeeParam
493           Int_t rob          = 0;     // ROB from AliTRDfeeParam
494           Int_t extCol       = 0;     // extended column from  AliTRDfeeParam  
495           mcm = feeParam->GetMCMfromPad(iRow,iCol);
496           rob = feeParam->GetROBfromPad(iRow,iCol);
497
498           Int_t idetector  = det;                            //  current detector
499           Int_t iRowMax    = rawStream->GetMaxRow();         //  current rowmax
500
501           Int_t adc        = 20 - (iCol%18) -1;                 //  current adc
502           Int_t col        = 0;                              //  col!=0 ->Shared Pad
503           extCol = feeParam->GetExtendedPadColFromADC(rob,mcm,adc);
504           //printf("  iCol %d  iRow %d  iRowMax %d  rob %d  mcm %d  adc %d  extCol %d\n",iCol,iRow,iRowMax,rob,mcm,adc,extCol);   
505           
506           // Signal for regular pads
507           Int_t nbtimebin  = digitParam->GetNTimeBins(idetector);  //  number of time bins read from data         
508           for(Int_t k = 0; k < nbtimebin; k++){
509             Short_t signal = 0;
510             signal = digits->GetData(iRow,iCol,k);
511             
512             if(signal>0) {
513               UpdateHisto2(idetector,iRow,iCol,signal,iRowMax,col,mcm,rob);
514             }
515           }
516
517           
518           
519           if((adc==3-1 || adc==20-1 || adc==19-1) && (iCol > 1 && iCol <142) /* && fSharedPadsOn*/ ) { //SHARED PADS
520             
521               switch(adc) {
522               case 2:  
523                 adc = 20;                                       //shared Pad adc 
524                 mcm = feeParam->GetMCMfromSharedPad(iRow,iCol); //shared Pad mcm 
525                 col =  1;
526                 break;
527               case 19:  
528                 adc = 1;                                        //shared Pad adc  
529                 mcm = feeParam->GetMCMfromSharedPad(iRow,iCol); //shared Pad mcm  
530                 col =  2;
531                 break;
532               case 18: 
533                 adc =  0;                                       //shared Pad adc  
534                 mcm = feeParam->GetMCMfromSharedPad(iRow,iCol); //shared Pad mcm 
535                 col =  3;
536                 break;
537               }
538               rob = feeParam->GetROBfromSharedPad(iRow,iCol);     //shared Pad rob 
539               
540             
541             extCol = feeParam->GetExtendedPadColFromADC(rob,mcm,adc);     //extended pad col via the shared pad rob,mcm and adc
542             
543             //printf("SHARED PAD ---  iCol %d  iRow %d  rob %d  mcm %d  adc %d  extCol %d  col %d\n",iCol,iRow,rob,mcm,adc,extCol,col);
544             for(Int_t k = 0; k < nbtimebin; k++){
545               Short_t signal = 0;
546               signal = digits->GetDataByAdcCol(iRow,extCol,k);
547               
548               if(signal>0) {
549                 UpdateHisto2(idetector,iRow,iCol,signal,iRowMax,col,mcm,rob);
550               }
551             }
552           } //shared pads end
553             
554           
555           withInput = 2;
556         }//column,row
557
558       }//array
559     }//QA
560     digitsManager->ClearArrays(det);
561   }//idetector
562   delete digitsManager;
563   delete rawStream;
564   return withInput;
565 }
566
567 //_____________________________________________________________________
568 Bool_t AliTRDCalibPadStatus::TestEventHisto(Int_t nevent, Int_t sm, Int_t ch) /*FOLD00*/
569 {
570   //
571   //  Test event loop
572   // fill one oroc and one iroc with random gaus
573   //
574
575   gRandom->SetSeed(0);
576
577     for (Int_t ism=sm; ism<sm+1; ism++){
578         for (Int_t ich=ch; ich < ch+1; ich++){
579             for (Int_t ipl=0; ipl < 6; ipl++){
580               for(Int_t irow = 0; irow < fGeo->GetRowMax(ipl,ich,ism); irow++){
581                 for(Int_t icol = 0; icol < fGeo->GetColMax(ipl); icol++){
582                   for (Int_t iTimeBin=0; iTimeBin<(30*nevent); iTimeBin++){
583                     Int_t signal=TMath::Nint(gRandom->Gaus(10,1.5));
584                     if ( signal>0 )UpdateHisto((ipl+ich*6+ism*6*5),irow,icol,signal,fGeo->GetRowMax(ipl,ich,ism),0,0);
585                   }
586                 }
587               }
588             }
589         }
590     }
591     return kTRUE;
592 }
593
594 //_____________________________________________________________________
595 TH2F* AliTRDCalibPadStatus::GetHisto(Int_t det, TObjArray *arr, /*FOLD00*/
596                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
597                                   const Char_t *type, Bool_t force)
598 {
599     //
600     // return pointer to histogram
601     // if force is true create a new histogram if it doesn't exist allready
602     //
603     if ( !force || arr->UncheckedAt(det) )
604         return (TH2F*)arr->UncheckedAt(det);
605
606     // if we are forced and histogram doesn't yes exist create it
607     Char_t name[255], title[255];
608
609     sprintf(name,"hCalib%s%.3d",type,det);
610     sprintf(title,"%s calibration histogram detector %.2d;ADC channel;Channel (pad)",type,det);
611
612    
613     Int_t nbchannels = fGeo->GetRowMax(GetLayer(det),GetStack(det),GetSector(det))*fGeo->GetColMax(GetLayer(det));
614     
615     // we will add 3*8*rowMax channels at the end for the double counted
616     nbchannels += 3*8*(fGeo->GetRowMax(GetLayer(det),GetStack(det),GetSector(det)));
617
618
619     // new histogram with calib information. One value for each pad!
620     TH2F* hist = new TH2F(name,title,
621                           nbinsY, ymin, ymax,
622                           nbchannels,0,nbchannels
623                           );
624     hist->SetDirectory(0);
625     arr->AddAt(hist,det);
626     return hist;
627 }
628
629 //_____________________________________________________________________
630 TH2F* AliTRDCalibPadStatus::GetHisto(Int_t det, Bool_t force) /*FOLD00*/
631 {
632     //
633     // return pointer to histogram
634     // if force is true create a new histogram if it doesn't exist allready
635     //
636     TObjArray *arr = &fHistoArray;
637     return GetHisto(det, arr, fAdcMax-fAdcMin, fAdcMin-0.5, fAdcMax-0.5, "Pedestal", force);
638 }
639
640 //_____________________________________________________________________
641 AliTRDCalROC* AliTRDCalibPadStatus::GetCalRoc(Int_t det, TObjArray* arr, Bool_t force) /*FOLD00*/
642 {
643     //
644     // return pointer to ROC Calibration
645     // if force is true create a new AliTRDCalROC if it doesn't exist allready
646     //
647     if ( !force || arr->UncheckedAt(det) )
648         return (AliTRDCalROC*)arr->UncheckedAt(det);
649
650     // if we are forced and histogram doesn't yes exist create it
651
652     // new AliTRDCalROC. One value for each pad!
653     AliTRDCalROC *croc = new AliTRDCalROC(GetLayer(det),GetStack(det));
654     arr->AddAt(croc,det);
655     return croc;
656 }
657 //_____________________________________________________________________
658 AliTRDCalROC* AliTRDCalibPadStatus::GetCalRocMean(Int_t det, Bool_t force) /*FOLD00*/
659 {
660     //
661     // return pointer to Carge ROC Calibration
662     // if force is true create a new histogram if it doesn't exist allready
663     //
664     TObjArray *arr = &fCalRocArrayMean;
665     return GetCalRoc(det, arr, force);
666 }
667
668 //_____________________________________________________________________
669 AliTRDCalROC* AliTRDCalibPadStatus::GetCalRocRMS(Int_t det, Bool_t force) /*FOLD00*/
670 {
671     //
672     // return pointer to Carge ROC Calibration
673     // if force is true create a new histogram if it doesn't exist allready
674     //
675     TObjArray *arr = &fCalRocArrayRMS;
676     return GetCalRoc(det, arr, force);
677 }
678 //_____________________________________________________________________
679 AliTRDCalROC* AliTRDCalibPadStatus::GetCalRocMeand(Int_t det, Bool_t force) /*FOLD00*/
680 {
681     //
682     // return pointer to Carge ROC Calibration
683     // if force is true create a new histogram if it doesn't exist allready
684     //
685     TObjArray *arr = &fCalRocArrayMeand;
686     return GetCalRoc(det, arr, force);
687 }
688
689 //_____________________________________________________________________
690 AliTRDCalROC* AliTRDCalibPadStatus::GetCalRocRMSd(Int_t det, Bool_t force) /*FOLD00*/
691 {
692     //
693     // return pointer to Carge ROC Calibration
694     // if force is true create a new histogram if it doesn't exist allready
695     //
696     TObjArray *arr = &fCalRocArrayRMSd;
697     return GetCalRoc(det, arr, force);
698 }
699
700 //_____________________________________________________________________
701 void AliTRDCalibPadStatus::AnalyseHisto() /*FOLD00*/
702 {
703     //
704     //  Calculate calibration constants
705     //
706
707     Int_t nbinsAdc = fAdcMax-fAdcMin;
708
709     TVectorD param(4);
710     TMatrixD dummy(3,3);
711
712     Float_t *arrayHP=0;
713
714
715     for (Int_t idet=0; idet<540; idet++){
716         TH2F *hP = GetHisto(idet);
717         if ( !hP ) {
718           continue;
719         }
720
721         //printf("Entries for %d\n",idet);
722
723         AliTRDCalROC *rocMean     = GetCalRocMean(idet,kTRUE);
724         AliTRDCalROC *rocRMS      = GetCalRocRMS(idet,kTRUE);
725
726         arrayHP = hP->GetArray();
727         Int_t nChannels = rocMean->GetNchannels();
728
729         for (Int_t iChannel=0; iChannel<nChannels; iChannel++){
730             Int_t offset = (nbinsAdc+2)*(iChannel+1)+1;
731             Double_t ret = AliMathBase::FitGaus(arrayHP+offset,nbinsAdc,fAdcMin-0.5,fAdcMax-0.5,&param,&dummy);
732             // if the fitting failed set noise and pedestal to 0
733             if ((ret==-4) || (ret==-1) || (ret==-2)) {
734                 param[1]=0.0;
735                 param[2]=0.0;
736             }
737             if((param[1]/10.0) > 65534.0) param[1] = 0.0;
738             if((param[2]/10.0) > 65534.0) param[2] = 0.0;
739             rocMean->SetValue(iChannel,param[1]/10.0);
740             rocRMS->SetValue(iChannel,param[2]/10.0);
741         }
742
743         // here we analyse doubled read channels
744         
745         AliTRDCalROC *rocMeand     = GetCalRocMeand(idet,kTRUE);
746         AliTRDCalROC *rocRMSd      = GetCalRocRMSd(idet,kTRUE);
747
748         Int_t nrows = rocMeand->GetNrows();
749         Int_t shift = 144*nrows;
750         Int_t total = shift+3*8*nrows; 
751
752         for (Int_t iChannel=shift; iChannel<total; iChannel++){
753             Int_t offset = (nbinsAdc+2)*(iChannel+1)+1;
754             Double_t ret = AliMathBase::FitGaus(arrayHP+offset,nbinsAdc,fAdcMin-0.5,fAdcMax-0.5,&param,&dummy);
755             // if the fitting failed set noise and pedestal to 0
756             if ((ret==-4) || (ret==-1) || (ret==-2)) {
757                 param[1]=0.0;
758                 param[2]=0.0;
759             }
760             if((param[1]/10.0) > 65534.0) param[1] = 0.0;
761             if((param[2]/10.0) > 65534.0) param[2] = 0.0;
762             
763             // here we have to recalculate backward
764             Int_t nb   = iChannel-shift;
765             Int_t row  = nb%nrows;
766             Int_t j    = (Int_t)(nb/nrows);
767             Int_t imcm = j%8;
768             Int_t icol = (Int_t)(j/8);
769             
770             Int_t finalcol = 18*imcm;
771             if(icol > 0) icol += 17;
772             else icol = -1;
773             finalcol += icol;
774
775             Int_t channel = row+finalcol*nrows;
776
777             //printf("iChannel %d, nrows %d, finalcol %d, row %d, channel %d\n",iChannel,nrows,finalcol,row,channel);
778             if((finalcol < 0) || (finalcol >= 144)) continue;
779
780             rocMeand->SetValue(channel,param[1]/10.0);
781             rocRMSd->SetValue(channel,param[2]/10.0);
782         }
783
784     }
785       
786 }
787
788 //_______________________________________________________________________________________
789 AliTRDCalPadStatus* AliTRDCalibPadStatus::CreateCalPadStatus()
790 {
791   //
792   // Create Pad Status out of Mean and RMS values
793   // The chamber without data are masked, this is the corrected in the preprocessor
794   //
795
796   AliTRDCalPadStatus* obj = new AliTRDCalPadStatus("padstatus", "padstatus");
797   
798   for (Int_t idet=0; idet<540; ++idet)
799     {
800       AliTRDCalSingleChamberStatus *calROC = obj->GetCalROC(idet);
801
802
803       if ( !GetCalRocMean(idet)) {
804         for(Int_t k = 0; k < calROC->GetNchannels(); k++){
805           calROC->SetStatus(k,AliTRDCalPadStatus::kNotConnected);
806         }
807         continue;
808       }
809       
810
811       //Take the stuff
812       AliTRDCalROC *calRocMean    = new AliTRDCalROC(*( (AliTRDCalROC *) GetCalRocMean(idet)));
813       AliTRDCalROC *calRocRMS     = new AliTRDCalROC(*( (AliTRDCalROC *) GetCalRocRMS(idet)));
814
815       //Take the stuff second chance
816       AliTRDCalROC *calRocMeand    = new AliTRDCalROC(*( (AliTRDCalROC *) GetCalRocMeand(idet)));
817       AliTRDCalROC *calRocRMSd     = new AliTRDCalROC(*( (AliTRDCalROC *) GetCalRocRMSd(idet)));
818
819       calRocRMS->Unfold();
820       calRocRMSd->Unfold();
821
822      
823       //Range
824       Int_t row      = calROC->GetNrows();
825       Int_t col      = calROC->GetNcols();
826       
827       Double_t rmsmean       = calRocMean->GetRMS()*10.0;
828       Double_t meanmean      = calRocMean->GetMean()*10.0;
829       Double_t meansquares   = calRocRMS->GetMean();
830
831       
832       for(Int_t irow = 0; irow < row; irow++){
833         
834         // for bridged pads
835         Float_t meanprevious = 0.0;
836         Float_t rmsprevious  = 0.0; 
837         Float_t mean         = 0.0;
838         Float_t rms          = 0.0;
839         
840         for(Int_t icol = 0; icol < col; icol++){
841           
842           mean     = calRocMean->GetValue(icol,irow)*10.0;
843           rms      = calRocRMS->GetValue(icol,irow);
844
845           if(icol > 0) {
846             meanprevious     = calRocMean->GetValue((icol -1),irow)*10.0;
847             rmsprevious      = calRocRMS->GetValue((icol - 1),irow);
848           }
849           
850           Bool_t pb = kFALSE;
851           // masked if two noisy
852           if((rms <= 0.0001) || (TMath::Abs(mean-meanmean)>(5*rmsmean)) || (TMath::Abs(rms)>(5.0*TMath::Abs(meansquares)))) {
853             
854             pb = kTRUE;
855             // look at second chance
856             Float_t meand     = calRocMeand->GetValue(icol,irow)*10.0;
857             Float_t rmsd      = calRocRMSd->GetValue(icol,irow);
858             
859             if((rmsd <= 0.0001) || (TMath::Abs(meand-meanmean)>(5*rmsmean)) || (TMath::Abs(rmsd)>(5.0*TMath::Abs(meansquares)))) {
860               if((rmsd <= 0.0001) && (rms <= 0.0001)) {
861                 calROC->SetStatus(icol,irow,AliTRDCalPadStatus::kNotConnected);
862               }
863               else {
864                 calROC->SetStatus(icol, irow, AliTRDCalPadStatus::kMasked);
865               }
866             }
867             else {
868               calROC->SetStatus(icol, irow, AliTRDCalPadStatus::kReadSecond);
869             }
870           }
871
872
873           // bridge if previous pad found something
874           if(!pb) {
875             if((meanprevious == mean) && (rmsprevious == rms) && (mean > 0.0001)) {
876               //printf("mean previous %f, mean %f, rms %f, rmsprevious %f, col %d\n",meanprevious,mean,rms,rmsprevious,icol);
877               calROC->SetStatus(icol -1 ,irow, AliTRDCalPadStatus::kPadBridgedRight);
878               calROC->SetStatus(icol ,irow, AliTRDCalPadStatus::kPadBridgedLeft);
879             }       
880           }
881
882         }
883       }
884
885       delete calRocMean;
886       delete calRocRMS;
887       delete calRocMeand;
888       delete calRocRMSd;
889
890
891     }
892   
893   return obj;
894   
895 }
896 //_______________________________________________________________________________________
897 AliTRDCalPad* AliTRDCalibPadStatus::CreateCalPad()
898 {
899   //
900   // Create Pad Noise out of RMS values
901   //
902
903   AliTRDCalPad* obj = new AliTRDCalPad("PadNoise", "PadNoise");
904   
905   
906   for (Int_t det=0; det<AliTRDgeometry::kNdet; ++det)  {
907     
908     AliTRDCalROC *calROC22 = obj->GetCalROC(det);
909
910     AliTRDCalROC *calRocRMS     = ((AliTRDCalROC *)GetCalRocRMS(det,kTRUE));
911    
912     for(Int_t k = 0; k < calROC22->GetNchannels(); k++){
913       calROC22->SetValue(k,calRocRMS->GetValue(k));
914     }
915
916   }
917   
918   return obj;
919   
920 }
921
922 //_______________________________________________________________________________________
923 AliTRDCalDet* AliTRDCalibPadStatus::CreateCalDet() const
924 {
925   //
926   // Create Det Noise correction factor
927   //
928
929   AliTRDCalDet* obj = new AliTRDCalDet("DetNoise", "DetNoise (correction factor)");
930
931   for(Int_t l = 0; l < 540; l++){
932     obj->SetValue(l,10.0);
933   }
934   
935   return obj;
936   
937 }
938
939 //_____________________________________________________________________
940 void AliTRDCalibPadStatus::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append) /*FOLD00*/
941 {
942     //
943     //  Write class to file
944     //
945
946     TString sDir(dir);
947     TString option;
948
949     if ( append )
950         option = "update";
951     else
952         option = "recreate";
953
954     TDirectory *backup = gDirectory;
955     TFile f(filename,option.Data());
956     f.cd();
957     if ( !sDir.IsNull() ){
958         f.mkdir(sDir.Data());
959         f.cd(sDir);
960     }
961     this->Write();
962     f.Close();
963
964     if ( backup ) backup->cd();
965 }
966
967 //_____________________________________________________________________
968 void AliTRDCalibPadStatus::SetCalRocMean(AliTRDCalROC *mean, Int_t det) /*FOLD00*/
969 {
970     //
971     //  Put the AliTRDCalROC in the array fCalRocArrayMean
972     //
973
974
975   AliTRDCalROC *rocMean = GetCalRocMean(det,kTRUE);
976   
977   Int_t nChannels = rocMean->GetNchannels();
978   
979   for (Int_t iChannel=0; iChannel<nChannels; iChannel++){
980     
981     rocMean->SetValue(iChannel,mean->GetValue(iChannel));
982     
983   }
984   
985 }
986
987 //_____________________________________________________________________
988 void AliTRDCalibPadStatus::SetCalRocRMS(AliTRDCalROC *rms, Int_t det) /*FOLD00*/
989 {
990     //
991     //  Put the AliTRDCalROC in the array fCalRocArrayRMS
992     //
993
994
995   AliTRDCalROC *rocRms = GetCalRocRMS(det,kTRUE);
996   
997   Int_t nChannels = rocRms->GetNchannels();
998   
999   for (Int_t iChannel=0; iChannel<nChannels; iChannel++){
1000     
1001     rocRms->SetValue(iChannel,rms->GetValue(iChannel));
1002     
1003   }
1004   
1005 }
1006 //_____________________________________________________________________
1007 void AliTRDCalibPadStatus::SetCalRocMeand(AliTRDCalROC *mean, Int_t det) /*FOLD00*/
1008 {
1009     //
1010     //  Put the AliTRDCalROC in the array fCalRocArrayMean
1011     //
1012
1013
1014   AliTRDCalROC *rocMean = GetCalRocMeand(det,kTRUE);
1015   
1016   Int_t nChannels = rocMean->GetNchannels();
1017   
1018   for (Int_t iChannel=0; iChannel<nChannels; iChannel++){
1019     
1020     rocMean->SetValue(iChannel,mean->GetValue(iChannel));
1021     
1022   }
1023   
1024 }
1025
1026 //_____________________________________________________________________
1027 void AliTRDCalibPadStatus::SetCalRocRMSd(AliTRDCalROC *rms, Int_t det) /*FOLD00*/
1028 {
1029     //
1030     //  Put the AliTRDCalROC in the array fCalRocArrayRMS
1031     //
1032
1033
1034   AliTRDCalROC *rocRms = GetCalRocRMSd(det,kTRUE);
1035   
1036   Int_t nChannels = rocRms->GetNchannels();
1037   
1038   for (Int_t iChannel=0; iChannel<nChannels; iChannel++){
1039     
1040     rocRms->SetValue(iChannel,rms->GetValue(iChannel));
1041     
1042   }
1043   
1044 }
1045 //_____________________________________________________________________________
1046 Int_t AliTRDCalibPadStatus::GetLayer(Int_t d) const
1047 {
1048   //
1049   // Reconstruct the layer number from the detector number
1050   //
1051
1052   return ((Int_t) (d % 6));
1053
1054 }
1055
1056 //_____________________________________________________________________________
1057 Int_t AliTRDCalibPadStatus::GetStack(Int_t d) const
1058 {
1059   //
1060   // Reconstruct the chamber number from the detector number
1061   //
1062
1063   return ((Int_t) (d % 30) / 6);
1064
1065 }
1066
1067 //_____________________________________________________________________________
1068 Int_t AliTRDCalibPadStatus::GetSector(Int_t d) const
1069 {
1070   //
1071   // Reconstruct the sector number from the detector number
1072   //
1073
1074   return ((Int_t) (d / 30));
1075
1076 }
1077