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