PHOS/DA/CMakeLists.txt
[u/mrichter/AliRoot.git] / PHOS / PHOSbase / AliPHOSCpvGainCalibDA.cxx
1 #include "AliPHOSCpvGainCalibDA.h" 
2 #include "AliPHOSCpvParam.h" 
3 #include "AliPHOSCpvRawDigiProducer.h"
4 #include "AliPHOSDigit.h"
5 #include <fstream>
6 #include <iostream>
7 #include <TTree.h>
8 #include <TF1.h>
9 #include <TFitResult.h>
10 #include <TFitResultPtr.h>
11 #include <TSystem.h>
12 #include <TTimeStamp.h>
13 #include "AliPHOSGeometry.h"
14
15 #include "TFile.h"
16
17 ClassImp(AliPHOSCpvGainCalibDA) ;
18
19 using namespace std;
20
21 using std::ifstream;
22 using std::ofstream;
23  
24
25
26 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
27 AliPHOSCpvGainCalibDA::AliPHOSCpvGainCalibDA():
28   TObject(),
29   fGeom(0)
30 {
31   //
32   //constructor
33   //
34   
35   fGeom=AliPHOSGeometry::GetInstance() ;
36   if(!fGeom) fGeom = AliPHOSGeometry::GetInstance("IHEP");
37
38   for(Int_t iDDL = 0; iDDL<2*AliPHOSCpvParam::kNDDL; iDDL++){
39     fDeadMap[iDDL] = 0x0;
40     fEntriesMap[iDDL] = 0x0;
41     for(Int_t iX=0; iX<AliPHOSCpvParam::kPadPcX; iX++) 
42       for(Int_t iY=1; iY<AliPHOSCpvParam::kPadPcY; iY++) 
43         fAmplA0Histo[iDDL][iX][iY] = 0x0;
44   }
45   CreateQAHistos();
46 }  //constructor
47 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
48 AliPHOSCpvGainCalibDA::~AliPHOSCpvGainCalibDA()
49 {
50   //
51   //destructor
52   //
53   for(Int_t iDDL=0; iDDL<2*AliPHOSCpvParam::kNDDL; iDDL++) {
54     for(Int_t iX=0; iX<AliPHOSCpvParam::kPadPcX; iX++) {
55       for(Int_t iY=1; iY<AliPHOSCpvParam::kPadPcY; iY++) {
56         if(fAmplA0Histo[iDDL][iX][iY]) delete fAmplA0Histo[iDDL][iX][iY];
57       }//iY
58     }//iX
59     if(fDeadMap[iDDL])delete fDeadMap[iDDL];
60     if(fEntriesMap[iDDL]) delete fEntriesMap[iDDL];
61   }//iDDL
62
63   //delete fhErrors;
64 }  //destructor
65 //***********************************************************************
66 void AliPHOSCpvGainCalibDA::InitCalibration(TFile *fCalibrSupplyRoot){
67   //tries to open CpvCalibrSupply.root for loading previously filled histograms;
68   //creates new histos otherwise.
69   //TFile *fCalibrSupplyRoot = TFile::Open("CpvCalibrSupply.root");
70   for(Int_t iDDL = 0;iDDL<2*AliPHOSCpvParam::kNDDL;iDDL++){
71     if(fCalibrSupplyRoot)
72       if(fCalibrSupplyRoot->Get(Form("hEntriesMap%d",iDDL))) 
73         fEntriesMap[iDDL]=new TH2I(*(TH2I*)(fCalibrSupplyRoot->Get(Form("hEntriesMap%d",iDDL))));
74       else fEntriesMap[iDDL]=0x0;
75     else fEntriesMap[iDDL]=0x0;
76     for(Int_t iX = 0;iX  <AliPHOSCpvParam::kPadPcX;iX++)
77       for(Int_t iY = 0;iY  <AliPHOSCpvParam::kPadPcY;iY++)
78         if(fCalibrSupplyRoot){
79           if(fCalibrSupplyRoot->Get(Form("hAmplA0_DDL%d_iX%d_iY%d",iDDL,iX,iY)))
80             fAmplA0Histo[iDDL][iX][iY]=new TH1F(*(TH1F*)(fCalibrSupplyRoot->Get(Form("hAmplA0_DDL%d_iX%d_iY%d",iDDL,iX,iY))));
81           else fAmplA0Histo[iDDL][iX][iY]=0x0;
82         }
83         else fAmplA0Histo[iDDL][iX][iY]=0x0;
84   }
85 }
86 //***********************************************************************
87 void AliPHOSCpvGainCalibDA::CreateA0Histos(Int_t iDDL){
88   //create 1D histos for particular DDL to fill them with raw amplitudes later
89   if(iDDL<0||iDDL>=2*AliPHOSCpvParam::kNDDL) return; //invalid ddl number
90   fEntriesMap[iDDL]=new TH2I(Form("hEntriesMap%d",iDDL),Form("A0 entries map, DDL = %d",iDDL),AliPHOSCpvParam::kPadPcX,0,AliPHOSCpvParam::kPadPcX,AliPHOSCpvParam::kPadPcY,0,AliPHOSCpvParam::kPadPcY);
91   fHistosList->Add(fEntriesMap[iDDL]);
92   for(Int_t iX = 0;iX  <AliPHOSCpvParam::kPadPcX;iX++)
93     for(Int_t iY = 0;iY  <AliPHOSCpvParam::kPadPcY;iY++)
94       fAmplA0Histo[iDDL][iX][iY]=new TH1F(Form("hAmplA0_DDL%d_iX%d_iY%d",iDDL,iX,iY),Form("Max amplitude in cluster DDL = %d X = %d Y = %d",iDDL,iX,iY),4096,0.,4096.);
95 }
96 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
97 Bool_t AliPHOSCpvGainCalibDA::SetDeadChannelMapFromFile(const char * filename = "CpvBadMap.root"){
98   //
99   //Set Dead Channel Map Cut from the file, if the input file is not present default value is set!
100   //Arguments: the name of the Dead Channel Map file 
101   //Returns: false if not possible to open file with provided filename
102
103   TFile *fDeadCh = TFile::Open(filename);
104   if(!fDeadCh)return 0;
105   for(Int_t iDDL = 0; iDDL < 2*AliPHOSCpvParam::kNDDL; iDDL+=2){
106     if(fDeadCh->Get(Form("hBadChMap%d",iDDL)))
107       fDeadMap[iDDL] = new TH2I(*(TH2I*)(fDeadCh->Get(Form("hBadChMap%d",iDDL))));
108   }
109   fDeadCh->Close();
110   return 1;
111 }//SetDeadChannelMapFromFile()
112 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
113 void AliPHOSCpvGainCalibDA::WriteA0HistosToFile() const
114 {
115   //write all A0 amplitude histos and A0 entries maps to CpvCalibrSupply.root
116   TFile * rootF = new TFile("CpvCalibrSupply.root","RECREATE");
117   rootF->cd();
118   
119   for(Int_t iDDL=0; iDDL<2*AliPHOSCpvParam::kNDDL; iDDL+=2){
120     if(fEntriesMap[iDDL]) fEntriesMap[iDDL]->Write();
121     for(Int_t iX=0; iX<AliPHOSCpvParam::kPadPcX; iX++) 
122       for(Int_t iY=0; iY<AliPHOSCpvParam::kPadPcY; iY++) 
123         if(fAmplA0Histo[iDDL][iX][iY])fAmplA0Histo[iDDL][iX][iY]->Write();
124   }
125   rootF->Close();
126 } //WriteAllHistsToFile()
127 //*************************************************************
128 Bool_t AliPHOSCpvGainCalibDA::FillAmplA0Histos(TClonesArray *digits){
129   // do a clusterization then find cell with max amlitude (so called A0) in every cluster
130   // fill corresponding histos with A0 amplitude
131   // return true in case of success (found at least 1 cluster).
132   if(!digits) return kFALSE;
133   Int_t nDig = digits->GetEntriesFast();
134   if(nDig < 2) return kFALSE;//why do we need at least 2 digits?
135   Bool_t stop = kFALSE;
136   Int_t nExcludedPoints = 0;
137   Bool_t *excludedPoints = new Bool_t[nDig];//points which already belongs to other clusters
138   for(int i=0;i<nDig;i++)excludedPoints[i]=kFALSE;
139   Int_t clusterIndex[10][5][5];//10 clusters max; this array contains digit numbers which belongs to particular cluster
140   Int_t clusterDDL[10];// DDL number of particular cluster
141   Int_t clusterX[10]; //X coordinate of cluster 
142   Int_t clusterY[10]; //Y coordinate of cluster
143   Float_t clusterAmplitude[10][5][5];
144   //=============================================================================
145   //========================= C L U S T E R I S A T I O N =======================
146   //=============================================================================
147   //here we define cluster as bunch of cells with at least 1 common verticies
148 // x= |_0_|_1_|_2_|_3_|_4_|
149 // y=0|   |   |   |   |   |
150 //    |___|___|___|___|___|
151 // y=1|   |   |   |   |   |
152 //    |___|___|___|___|___|
153 // y=2|   |   |   |   |   |
154 //    |___|___|___|___|___|
155 // y=3|   |   |   |   |   |
156 //    |___|___|___|___|___|
157 // y=4|   |   |   |   |   |
158 //    |___|___|___|___|___|
159   // initialize clusters array
160   for(Int_t iClus=0;iClus<10;iClus++)
161     for(Int_t ix=0;ix<5;ix++)
162       for(Int_t iy=0;iy<5;iy++)
163         clusterIndex[iClus][ix][iy]=-1;
164   
165   Int_t relId[4];
166   Int_t cluNumber = 0;
167   while(!stop){//we are going to find 10 or less clusters
168     Float_t qMax = 0.;//local maximum value
169     Int_t indMax = -1;//local maximum index
170     for(Int_t iDig = 0; iDig<nDig;iDig++){//find a local maximum
171       if(excludedPoints[iDig]==kTRUE)continue;//is this point already excluded?
172       AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( digits->At(iDig) ) ;
173       fGeom->AbsToRelNumbering(digit->GetId(),relId);
174       if(relId[1]!=-1){//exclude this digit
175         nExcludedPoints++; 
176         excludedPoints[iDig]=kTRUE;
177         continue;//this is not a CPV digit
178       }
179       int DDL = AliPHOSCpvParam::Mod2DDL(relId[0]);
180       if(IsBad(DDL,relId[2],relId[3])){//let's see if it is a bad pad
181         nExcludedPoints++; 
182         excludedPoints[iDig]=kTRUE;
183         continue;
184       }
185       if( digit->GetEnergy()>qMax) {qMax = digit->GetEnergy(); indMax = iDig;}
186     }
187     if(indMax<0){//did we find something?
188       stop=kTRUE;
189       continue;//no new local maximum 
190     }
191     //set found local maximum as center of new cluster
192     nExcludedPoints++; excludedPoints[indMax]=kTRUE; //do not forget to exclude found maximum from consideration
193     AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( digits->At(indMax) ) ;
194     fGeom->AbsToRelNumbering(digit->GetId(),relId);
195     clusterIndex[cluNumber][2][2] = indMax;
196     clusterAmplitude[cluNumber][2][2] = qMax;
197     clusterDDL[cluNumber] = AliPHOSCpvParam::Mod2DDL(relId[0]);
198     clusterX[cluNumber]=relId[2];
199     clusterY[cluNumber]=relId[3];
200     //let us try to find all other digits which belongs to the same cluster
201     Int_t pointsFound = 0;
202     do{
203       pointsFound=0;
204       for(Int_t iDig = 0; iDig<nDig;iDig++){
205         if(excludedPoints[iDig]==kTRUE)continue;//is this point already excluded?
206         AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( digits->At(iDig) ) ;
207         fGeom->AbsToRelNumbering(digit->GetId(),relId);
208         if(AliPHOSCpvParam::Mod2DDL(relId[0])!=clusterDDL[cluNumber]) continue;
209         //see if this current digit has common vertex with 
210         for(int ix = 0;ix<5;ix++)
211           for(int iy = 0;iy<5;iy++)
212             if(clusterIndex[cluNumber][ix][iy]>=0&&(TMath::Abs(relId[2]-clusterX[cluNumber]+ix-2)<=1&&TMath::Abs(relId[3]-clusterY[cluNumber]+iy-2)<=1)){//we found a cell!
213               pointsFound++;
214               clusterAmplitude[cluNumber][ix][iy] = digit->GetEnergy();
215               clusterIndex[cluNumber][ix][iy] = iDig;
216               //finally, exclude this point from consideration
217               nExcludedPoints++;
218               excludedPoints[iDig]=kTRUE;
219             }
220       }
221     }while (pointsFound!=0);
222     //OK, we have finished with this cluster
223     cluNumber++;
224     if(cluNumber>=10) stop=kTRUE; //we found enough clusters
225     if(nExcludedPoints>=nDig) stop=kTRUE;//we assigned all the digits
226   }
227   //cout<<"I found " <<cluNumber<<" clusters"<<endl;
228
229   //now we can operate with clusters
230   //=====================================================================================
231   //===================== F I L L I N G == O F == H I S T O G R A M S====================
232   //=====================================================================================
233   
234   fhClusterMult->Fill(cluNumber);
235   for (Int_t iClu = 0; iClu < cluNumber; iClu++){
236     if(!fEntriesMap[clusterDDL[iClu]]) CreateA0Histos(clusterDDL[iClu]);
237     // cout<<"iClu = "<<iClu<<
238     fAmplA0Histo[clusterDDL[iClu]][clusterX[iClu]-1][clusterY[iClu]-1]->Fill(clusterAmplitude[iClu][2][2]);
239     fEntriesMap[clusterDDL[iClu]]->Fill(clusterX[iClu]-1,clusterY[iClu]-1);
240     fhA0Value->Fill(clusterAmplitude[iClu][2][2]);
241     Double_t totAmpl = 0.;
242     for(int ix = 0; ix<5; ix++)
243       for(int iy = 0; iy<5; iy++)
244         if(clusterIndex[iClu][ix][iy]>=0){
245           fhAmplInClust->Fill(clusterAmplitude[iClu][ix][iy],ix*5+iy);
246           fhClusterShape->Fill(ix-2,iy-2);
247           totAmpl+=clusterAmplitude[iClu][ix][iy];
248         }
249   }
250 }
251 //*************************************************************
252 void AliPHOSCpvGainCalibDA::CreateQAHistos(){
253   fHistosList=new TList();
254   
255   fhClusterMult = new TH1F("fhClusterMult","Cluster Multiplicity in event",100,0,100);
256   fHistosList->Add(fhClusterMult);
257   
258   fhClusterShape=new TH2F("fhClusterShape","Shape of cluster", 5,-2.5,2.5, 5,-2.5,2.5 );
259   fHistosList->Add(fhClusterShape);
260
261   fhA0Value = new TH1F("fhA0Value","Max Amplitude in Cluster ",4096,0.,4096);
262   fHistosList->Add(fhA0Value);
263   
264   fhAmplInClust=new TH2F("fhAmplInClust", "amplitude distribution in cluster", 1000,0.,1000., 25,0.,25.);
265   fHistosList->Add(fhAmplInClust);
266
267   fhTotalClusterAmplitude = new TH1F("fhTotalClusterAmplitude", "Total Amplitude in Cluster",4096,0.,4096);
268   fHistosList->Add(fhTotalClusterAmplitude);
269 }