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