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