Adding some further mother volumes to speed-up the overlap checking and particle...
[u/mrichter/AliRoot.git] / ITS / StoreCalibSDD.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include "AliCDBManager.h"
3 #include "AliITSCalibrationSDD.h"
4 #include "AliITSgeomTGeo.h"
5 #include "AliCDBMetaData.h"
6 #include "AliCDBStorage.h"
7 #include "AliCDBId.h"
8 #include "AliCDBPath.h"
9 #include "AliCDBEntry.h"
10 #include <TObjArray.h>
11 #include <TRandom3.h>
12 #endif
13
14 void StoreCalibSDD(Int_t firstRun=0,Int_t lastRun=9999999){
15   ///////////////////////////////////////////////////////////////////////
16   // Macro to generate and store the calibration files for SDD         //
17   // Generates:                                                        //
18   //  1 file with 260 AliITSCalibrationSDD objects with                //
19   //    baselines, noise, gain, drift speed for each module (CalibSDD) //
20   ///////////////////////////////////////////////////////////////////////
21   
22   if(!AliCDBManager::Instance()->IsDefaultStorageSet()) {
23     AliCDBManager::Instance()->SetDefaultStorage("local://OCDB");
24   }
25   
26
27   AliCDBMetaData *md= new AliCDBMetaData(); // metaData describing the object
28   md->SetObjectClassName("TObjArray");
29   md->SetResponsible("Francesco Prino");
30   md->SetBeamPeriod(0);
31   md->SetComment("Simulated data");
32
33   AliCDBId idCalSDD("ITS/Calib/CalibSDD",firstRun, lastRun);
34   TObjArray calSDD(260);
35   calSDD.SetOwner(kFALSE);
36
37   // BAD modules data 
38   const Int_t nbadmod=6;
39   Int_t idBadMod[nbadmod];
40   idBadMod[0]=AliITSgeomTGeo::GetModuleIndex(3,4,1);
41   idBadMod[1]=AliITSgeomTGeo::GetModuleIndex(3,4,2);
42   idBadMod[2]=AliITSgeomTGeo::GetModuleIndex(3,4,3);
43   idBadMod[3]=AliITSgeomTGeo::GetModuleIndex(4,1,1);
44   idBadMod[4]=AliITSgeomTGeo::GetModuleIndex(4,12,1);
45   idBadMod[5]=AliITSgeomTGeo::GetModuleIndex(4,21,3);
46   
47   // Modules with bad left side 
48   const Int_t nbadleft=5;
49   Int_t idBadLeft[nbadleft];
50   idBadLeft[0]=AliITSgeomTGeo::GetModuleIndex(3,13,6);
51   idBadLeft[1]=AliITSgeomTGeo::GetModuleIndex(4,4,8);
52   idBadLeft[2]=AliITSgeomTGeo::GetModuleIndex(4,5,3);
53   idBadLeft[3]=AliITSgeomTGeo::GetModuleIndex(4,18,4);
54   idBadLeft[4]=AliITSgeomTGeo::GetModuleIndex(4,21,5);
55
56
57  // BAD anodes data
58   const Int_t nData = 209;
59   Int_t anodeUp[209] = {0,36,0,12,20,32,0,0,12,76,28,8,16,0,0,0,8,0,0,0,20,4,0,0,0,0,0,0
60                         ,0,0,8,0,0,0,0,0,0,0,0,0,0,0,12,0,8,0,4,4,0,160,0,0,0,252,16,0,8,8
61                         ,0,0,12,0,0,0,12,0,15,20,0,0,0,0,0,0,0,12,12,0,0,0,0,0,8,8,3,240,212
62                         ,12,0,8,0,12,0,0,12,0,8,24,0,0,12,0,0,0,0,40,0,0,40,12,28,0,0,12,12
63                         ,0,0,0,0,20,0,0,0,12,0,24,0,0,0,0,0,0,0,8,16,0,0,0,0,0,256,0,0,0,0,0,20
64                         ,0,12,0,0,0,0,24,0,0,0,0,0,0,0,20,0,0,16,0,0,0,0,24,0,0,0,8,0,16,40,0
65                         ,0,0,0,0,0,0,0,0,4,0,32,12,8,28,0,76,0,0,0,12,60,144,0,0,0,0,16,0,16,0,3 };
66
67   Int_t anodeDown[209] = {0,8,0,48,0,12,0,8,0,80,12,0,4,4,0,0,24,0,0,20,0,0,0,0,20,0,0,0,0,0,0,0
68                           ,0,0,0,0,4,0,24,4,0,0,8,0,0,36,20,0,12,236,0,0,12,256,16,8,32,4,0,0,24
69                           ,24,10,0,16,8,0,2,40,0,0,0,24,0,0,0,8,0,0,0,0,0,32,240,0,92,60,28,8,0,0
70                           ,2,0,0,0,0,12,48,0,0,0,0,0,36,11,12,0,0,0,12,12,11,0,20,0,0,12,20,0,0,4
71                           ,0,8,12,0,0,0,16,16,0,32,72,12,0,88,20,16,112,8,0,244,28,256,28,0,24,236
72                           ,56,0,68,0,4,20,208,20,12,4,28,12,0,0,20,12,0,100,0,16,8,8,0,24,16,0,12,12
73                           ,16,0,16,20,0,28,0,8,24,0,12,8,4,40,0,104,96,32,140,20,12,8,20,24,16,16,20
74                           ,8,140,96,0,32,20,44};
75
76   TRandom3 *gran = new TRandom3();
77   
78   for(Int_t mod=0;mod<260;mod++){
79     AliITSCalibrationSDD* resd = new AliITSCalibrationSDD("simulated");
80
81     Int_t nBadUp = 0;
82     Int_t nBadDown = 0;
83       
84
85     // gain
86     for(Int_t iWing=0; iWing<2;iWing++){
87       for(Int_t iChip=0; iChip<4;iChip++){
88         Float_t chipgain=gran->Gaus(1.,0.1);
89         if(chipgain<0.1) chipgain=0.1;
90         for(Int_t iChan=0; iChan<64;iChan++){
91           Float_t gain=gran->Gaus(chipgain,0.01);
92           if(gain<0.1) gain=0.1;
93           Int_t ian=resd->GetAnodeNumber(iWing,iChip,iChan);
94           resd->SetGain(ian,gain);
95         }
96       }
97     }
98
99     // bad channels
100     Int_t count=0;
101     do {
102       Int_t ranMod = Int_t(nData*gran->Uniform());
103       nBadUp   = anodeUp[ranMod];
104       nBadDown = anodeDown[ranMod];
105     }while (nBadUp+nBadDown>25);
106
107     resd->SetDeadChannels(nBadDown+nBadUp);
108     Int_t remainingBad = nBadDown;
109     while (remainingBad>0) {
110       Int_t nBadChannel;
111       if (remainingBad<5) {
112         nBadChannel = remainingBad;
113       } else {
114         Int_t width = remainingBad-(5-1);
115         if (width>4) width = 4;
116         nBadChannel = 5 + Int_t(width*gran->Uniform());
117       }
118         
119
120       Int_t iChannelPos = Int_t( (4*64-nBadChannel)*gran->Uniform() );
121       //      Int_t iChip = iChannelPos/64;
122       //      Int_t iChan = iChannelPos - iChip*64;
123       if(resd->IsBadChannel(iChannelPos)) continue;     
124       Int_t *clus = new Int_t[nBadChannel];
125       Int_t ich = iChannelPos;
126       for(Int_t i=0;i<nBadChannel;i++){
127         clus[i]=ich;
128         ich++;
129       }
130
131       for(Int_t i=0;i<nBadChannel;i++){
132         if(resd->IsBadChannel(clus[i])) break;
133         resd->SetBadChannel(count,clus[i]);
134         count++;          
135       }
136       remainingBad -= nBadChannel;
137       delete [] clus;
138     }
139
140     // May happen that, due to overlapping clusters, we
141     // have less bad channels than requested
142     // Let's put the remaining one per one ...
143     Int_t nSeToBad = 0;
144     for (Int_t i=0; i<4; i++){
145       for(Int_t j=0;j<64;j++){
146         Int_t ian=resd->GetAnodeNumber(0,i,j);
147         if (resd->GetChannelGain(ian)<0.0001) nSeToBad++;
148       }
149     }
150     while (nSeToBad<nBadDown) {
151       Int_t i = Int_t(4*64*gran->Uniform());
152       if(resd->IsBadChannel(i)==kFALSE){
153         resd->SetBadChannel(count,i);
154         count++;
155         nSeToBad++;
156       }
157     }
158       
159     remainingBad = nBadUp;
160     while (remainingBad>0) {
161       Int_t nBadChannel;
162       if (remainingBad<5) {
163         nBadChannel = remainingBad;
164       } else {
165         Int_t width = remainingBad-(5-1);
166         if (width>4) width = 4;
167         nBadChannel = 5 + Int_t(width*gran->Uniform());
168       }
169
170       Int_t iChannelPos = Int_t( (4*64-nBadChannel)*gran->Uniform() );
171       //      Int_t iChip = iChannelPos/64;
172       //      Int_t iChan = iChannelPos - iChip*64;
173       if(resd->IsBadChannel(iChannelPos)) continue;
174         
175       Int_t *clus = new Int_t[nBadChannel];
176       Int_t ich = iChannelPos;
177       for(Int_t i=0;i<nBadChannel;i++){
178         clus[i]=ich+256;
179         ich++;
180       }
181       for(Int_t i=0;i<nBadChannel;i++){
182         
183         if(resd->IsBadChannel(clus[i])) break;  
184         resd->SetBadChannel(count,clus[i]);
185         count++;
186           
187       }
188       remainingBad -= nBadChannel;
189       delete [] clus;
190     }
191       
192     nSeToBad = 0;
193     for (Int_t i=0; i<4; i++){
194       for(Int_t j=0;j<64;j++){
195         Int_t ian=resd->GetAnodeNumber(1,i,j);
196         if (resd->GetChannelGain(ian)<0.0001) nSeToBad++;
197       }
198     }
199
200     while (nSeToBad<nBadUp) {
201       Int_t i = Int_t(4*64*gran->Uniform());
202       if(resd->IsBadChannel(i+256)==kFALSE){
203         resd->SetBadChannel(count,i+256);
204         count++;
205         nSeToBad++;
206       }
207     }
208       
209     //baselines
210     /*
211       for(Int_t nan=0;nan<512;nan++){
212         
213       Int_t baseline = resd->GetBaseline(0);
214       Double_t noise = resd->GetNoiseAfterElectronics(0);
215       resd->SetBaseline(nan,baseline+(Int_t)gran->Gaus(0,4)); //baseline di def. e' 20, ho messo variazione di 4
216       resd->SetNoiseAfterElectronics(nan,noise+(Double_t)gran->Gaus(0,0.6));              
217                           
218       }
219     */
220
221     Int_t modId=mod+240;
222     // Bad modules
223     for(Int_t ibadm=0; ibadm<nbadmod; ibadm++){
224       if(modId==idBadMod[ibadm]) resd->SetBad();
225     }
226     // Modules with bad left side
227     for(Int_t ibadl=0; ibadl<nbadleft; ibadl++){
228       if(modId==idBadLeft[ibadl]) for(Int_t ichip=0;ichip<4;ichip++) resd->SetChipBad(ichip);
229     }
230
231     // Modules with bad pascal chips
232     if( modId==AliITSgeomTGeo::GetModuleIndex(4,4,2) ){
233       resd->SetChipBad(0);
234       resd->SetChipBad(3);
235     }
236
237     calSDD.Add(resd);
238     printf("Added module %d\n",mod);
239         }
240     
241   FILE* out = fopen("deadchannels.dat","w");
242   for(Int_t i=0;i<260;i++){
243     fprintf(out,"MODULE=%d\n",i);
244     AliITSCalibrationSDD* cl = (AliITSCalibrationSDD*)calSDD.At(i);
245     Int_t ndead=cl->GetDeadChannels();
246     fprintf(out,"n %d\n",ndead);
247     for(Int_t n=0;n<ndead;n++){
248       fprintf(out,"%d\n",cl->GetBadChannel(n));
249     }   
250   }
251   fclose(out);
252   AliCDBManager::Instance()->GetDefaultStorage()->Put(&calSDD, idCalSDD, md);   
253 }