]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/Upgrade/macros/spaceChargeFluctuation.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
o test
[u/mrichter/AliRoot.git] / TPC / Upgrade / macros / spaceChargeFluctuation.C
CommitLineData
c0172f82 1/*
2 .x ~/rootlogon.C
3 .L $ALICE_ROOT/TPC/Upgrade/macros/spaceChargeFluctuation.C+
4// GenerateMapRaw("/hera/alice/local/filtered/alice/data/2010/LHC10e/000129527/raw/merged_highpt_1.root","output.root", 50);
5 GenerateMapRawIons( "/hera/alice/local/filtered/alice/data/2010/LHC10e/000129527/raw/merged_highpt_1.root","output.root", 25);
6
7 */
8#include "TMath.h"
9#include "TRandom.h"
10#include "TTreeStream.h"
11#include "TVectorD.h"
12#include "TCanvas.h"
13#include "TStopwatch.h"
14#include "AliTPCParam.h"
15#include "AliTPCcalibDB.h"
16#include "AliTPCAltroMapping.h"
17#include "AliAltroRawStream.h"
18#include "AliSysInfo.h"
19#include "AliTPCRawStreamV3.h"
20#include "AliCDBManager.h"
21#include "TGeoGlobalMagField.h"
22#include "AliMagF.h"
23#include "AliRawReaderRoot.h"
24#include "AliRawReader.h"
25#include "TH3.h"
26#include "TH2.h"
27
28
29TH1 * GenerateMapRawIons(const char *fileName="raw.root", const char *outputName="histo.root", Int_t maxEvents=25);
30
31
32void spaceChargeFluctuation(Int_t mode=0){
33 //
34 //
35 if (mode==0) GenerateMapRawIons();
36
37}
38
39void spaceChargeFluctuationMC(Int_t nframes){
40 //
41 // Toy MC to generate space charge fluctuation
42 //
43 //
44 TTreeSRedirector *pcstream = new TTreeSRedirector("spaceChargeFluctuation.root","recreate");
45 Double_t interactionRate=100000;
46 Double_t driftTime=0.1;
47
48 for (Int_t iframe=0; iframe<nframes; iframe++){
49 Int_t nevents=gRandom->Poisson(interactionRate*driftTime);
50 Int_t ntracksAll=0;
51 TVectorD vecAllPhi128(128);
52 TVectorD vecAllPhi36(36);
53 for (Int_t ievent=0; ievent<nevents; ievent++){
54 //
55 //
56 Int_t ntracks=gRandom->Exp(500)+gRandom->Exp(500); // to be taken form the distribution
57 ntracksAll+=ntracks;
58 for (Int_t isec=0; isec<128; isec++){
59 vecAllPhi128(isec)+=ntracks*gRandom->Rndm()/128;
60 }
61 for (Int_t isec=0; isec<36; isec++){
62 vecAllPhi36(isec)+=ntracks*gRandom->Rndm()/36;
63 }
64 }
65 (*pcstream)<<"ntracks"<<
66 "nevents="<<nevents<< // number of events withing time frame
67 "ntracksAll="<<ntracksAll<< // number of tracks within time frame
68 "vecAllPhi36="<<&vecAllPhi36<< // number of tracks in phi bin (36 bins) within time frame
69 "vecAllPhi128="<<&vecAllPhi128<< // number of tracks in phi bin (128 bins) within time frame
70 "\n";
71 }
72 delete pcstream;
73}
74
75
76// void spaceChargeFluctuationDraw(){
77// //
78// //
79// //
80// TFile *f = TFile::Open("spaceChargeFluctuation.root");
81// TTree * tree = (TTree*)f->Get("ntracks");
82// TCanvas * canvasFluc = new TCanvas("cavasFluc","canvasFluc");
83// }
84
85
86TH1 * GenerateMapRaw(const char *fileName, const char *outputName, Int_t maxEvents){
87 //
88 // Generate 3D maps of electrons
89 // different threshold considered
90 // Paramaters:
91 // fileName - name of input raw file
92 // outputName - name of output file with the space charge histograms
93 // maxEvents - grouping of the events
94 //
95 //
96 TTreeSRedirector * pcstream = new TTreeSRedirector(outputName, "recreate");
97
98 const char * ocdbpath =gSystem->Getenv("OCDB_PATH") ? gSystem->Getenv("OCDB_PATH"):"local://$ALICE_ROOT/OCDB/";
99 AliCDBManager * man = AliCDBManager::Instance();
100 man->SetDefaultStorage(ocdbpath);
101 man->SetRun(0);
102 TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG, AliMagF::kBeamTypepp, 2.76/2.));
103 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
104 AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters();
105 //
106 TStopwatch timer;
107 timer.Start();
108 TH3F * hisQ3D[3]={0};
109 TH2F * hisQ2D[3]={0};
110 TH1F * hisQ1D[3]={0};
111 //
112 for (Int_t ith=0; ith<3; ith++){
113 char chname[100];
114 snprintf(chname,100,"hisQ3D_Th%d",2*ith+2);
115 hisQ3D[ith] = new TH3F(chname, chname,180, 0,2*TMath::TwoPi(),param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) ,100,0,1100);
116 snprintf(chname,100,"hisQ2D_Th%d",2*ith+2);
117 hisQ2D[ith] = new TH2F(chname,chname,180, 0,2*TMath::TwoPi(), param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
118 //
119 snprintf(chname,100,"hisQ1D_Th%d",2*ith+2);
120 hisQ1D[ith] = new TH1F(chname,chname, param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
121 hisQ3D[ith]->SetDirectory(0);
122 hisQ1D[ith]->SetDirectory(0);
123 hisQ1D[ith]->SetDirectory(0);
124 }
125 //
126 //
127 AliRawReader *reader = new AliRawReaderRoot(fileName);
128 reader->Reset();
129 AliAltroRawStream* stream = new AliAltroRawStream(reader);
130 stream->SelectRawData("TPC");
131 Int_t evtnr=0;
132 Int_t chunkNr=0;
133 //
134 while (reader->NextEvent()) {
135 if(evtnr>maxEvents) {
136 chunkNr++;
137 evtnr=0;
138 pcstream->GetFile()->mkdir(Form("Chunk%d",chunkNr));
139 pcstream->GetFile()->cd(Form("Chunk%d",chunkNr));
140 for (Int_t ith=0; ith<3; ith++){
141 hisQ1D[ith]->Write(Form("His1D_%d",ith));
142 hisQ2D[ith]->Write(Form("His2D_%d",ith));
143 hisQ3D[ith]->Write(Form("His3D_%d",ith));
5c70ec78 144 (*pcstream)<<"histo"<<
145 Form("hist1D_%d.=",ith*2+2)<<hisQ1D[ith]<<
146 Form("hist2D_%d.=",ith*2+2)<<hisQ2D[ith]<<
147 Form("hist3D_%d.=",ith*2+2)<<hisQ2D[ith]<<
148 "\n";
c0172f82 149 hisQ1D[ith]->Reset();
150 hisQ2D[ith]->Reset();
151 hisQ3D[ith]->Reset();
152 }
153 }
154 cout<<"Chunk=\t"<<chunkNr<<"\tEvt=\t"<<evtnr<<endl;
155 evtnr++;
156 AliSysInfo::AddStamp(Form("Event%d",evtnr),evtnr);
157 AliTPCRawStreamV3 input(reader,(AliAltroMapping**)mapping);
158 //
159 while (input.NextDDL()){
160 Int_t sector = input.GetSector();
161 while ( input.NextChannel() ) {
162 Int_t row = input.GetRow();
163 Int_t rowAbs = row+ ((sector<36) ? 0: param->GetNRow(0));
164 Int_t pad = input.GetPad();
165 Int_t nPads = param->GetNPads(sector,row);
166 Double_t localX = param->GetPadRowRadii(sector,row);
167 Double_t localY = (pad-nPads)*param->GetPadPitchWidth(sector);
168 Double_t localPhi= TMath::ATan2(localY,localX);
169 Double_t phi = TMath::Pi()*((sector%18)+0.5)/9+localPhi;
170 if (sector%36>=18) phi+=TMath::TwoPi();
171 Double_t padLength=param->GetPadPitchLength(sector,row);
172 // Double_t padWidth=param->GetPadPitchWidth(sector);
173 // Double_t zLength=param->GetZLength();
174 Double_t deltaPhi=hisQ3D[0]->GetXaxis()->GetBinWidth(0);
175 Double_t deltaZ= hisQ3D[0]->GetZaxis()->GetBinWidth(0);
176 //
177 while ( input.NextBunch() ){
178 Int_t startTbin = (Int_t)input.GetStartTimeBin();
179 Int_t bunchlength = (Int_t)input.GetBunchLength();
180 const UShort_t *sig = input.GetSignals();
181 Int_t aboveTh[3]={0};
182 for (Int_t i=0; i<bunchlength; i++){
183 for (Int_t ith=0; ith<3; ith++){
184 if (sig[i]>(ith*2)+2) aboveTh[ith]++;
185 }
186 }
187 for (Int_t ith=0; ith<3; ith++){
188 if (aboveTh[ith%3]>1){
189 for (Int_t i=0; i<bunchlength; i++){
190 //
191 // normalization
192 //
193 Double_t volume3D = padLength*(localX*deltaPhi)*deltaZ;
194 if (sector%36<18) hisQ1D[ith]->Fill(rowAbs, sig[i]/volume3D);
195 hisQ2D[ith]->Fill(phi,rowAbs, sig[i]/volume3D);
196 hisQ3D[ith]->Fill(phi,rowAbs,startTbin+i,sig[i]/volume3D);
197 }
198 }
199 }
200 }
201 }
202 }
203 }
204 timer.Print();
205 delete pcstream;
206 return 0;
207}
208
209
210
211
212
213TH1 * GenerateMapRawIons(const char *fileName, const char *outputName, Int_t maxEvents){
214 //
215 // Generate 3D maps of the space charge for the rad data maps
216 // different threshold considered
217 // Paramaters:
218 // fileName - name of input raw file
219 // outputName - name of output file with the space charge histograms
220 // maxEvents - grouping of the events
221 //
222 //
223 TTreeSRedirector * pcstream = new TTreeSRedirector(outputName, "recreate");
224 const char * ocdbpath =gSystem->Getenv("OCDB_PATH") ? gSystem->Getenv("OCDB_PATH"):"local://$ALICE_ROOT/OCDB/";
225 AliCDBManager * man = AliCDBManager::Instance();
226 man->SetDefaultStorage(ocdbpath);
227 man->SetRun(0);
228 TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG, AliMagF::kBeamTypepp, 2.76/2.));
229 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
230 AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters();
231 TStopwatch timer;
232 timer.Start();
233 // arrays of space charges - different elements corresponds to different threshold to accumulate charge
234 TH3F * hisQ3D[3]={0}; // 3D maps space charge from drift volume
235 TH2F * hisQ2D[3]={0};
236 TH1F * hisQ1D[3]={0};
237 TH3F * hisQ3DROC[3]={0}; // 3D maps space charge from ROC
238 TH2F * hisQ2DROC[3]={0};
239 TH1F * hisQ1DROC[3]={0};
5c70ec78 240 Int_t nbinsRow=param->GetNRowLow()+param->GetNRowUp();
241 Double_t *xbins = new Double_t[nbinsRow+1];
242 xbins[0]=param->GetPadRowRadiiLow(0)-1; //underflow bin
243 for (Int_t ibin=0; ibin<param->GetNRowLow();ibin++) xbins[1+ibin]=param->GetPadRowRadiiLow(ibin);
244 for (Int_t ibin=0; ibin<param->GetNRowUp();ibin++) xbins[1+ibin+param->GetNRowLow()]=param->GetPadRowRadiiUp(ibin);
245
246
c0172f82 247 //
248 for (Int_t ith=0; ith<3; ith++){
249 char chname[100];
250 snprintf(chname,100,"hisQ3D_Th%d",2*ith+2);
251 hisQ3D[ith] = new TH3F(chname, chname,180, 0,2*TMath::TwoPi(),param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) ,125,-250,250);
252 snprintf(chname,100,"hisQ2D_Th%d",2*ith+2);
253 hisQ2D[ith] = new TH2F(chname,chname,180, 0,2*TMath::TwoPi(), param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
254 snprintf(chname,100,"hisQ1D_Th%d",2*ith+2);
255 hisQ1D[ith] = new TH1F(chname,chname, param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
256 //
257 snprintf(chname,100,"hisQ3DROC_Th%d",2*ith+2);
258 hisQ3DROC[ith] = new TH3F(chname, chname,180, 0,2*TMath::TwoPi(),param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) ,125,-250,250);
259 snprintf(chname,100,"hisQ2DROC_Th%d",2*ith+2);
260 hisQ2DROC[ith] = new TH2F(chname,chname,180, 0,2*TMath::TwoPi(), param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
261 snprintf(chname,100,"hisQ1DROC_Th%d",2*ith+2);
262 hisQ1DROC[ith] = new TH1F(chname,chname, param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
5c70ec78 263 hisQ1D[ith]->GetXaxis()->Set(nbinsRow,xbins);
264 hisQ2D[ith]->GetYaxis()->Set(nbinsRow,xbins);
265 hisQ3D[ith]->GetYaxis()->Set(nbinsRow,xbins);
266 hisQ1DROC[ith]->GetXaxis()->Set(nbinsRow,xbins);
267 hisQ2DROC[ith]->GetYaxis()->Set(nbinsRow,xbins);
268 hisQ3DROC[ith]->GetYaxis()->Set(nbinsRow,xbins);
c0172f82 269 //
270 hisQ3D[ith]->SetDirectory(0);
271 hisQ1D[ith]->SetDirectory(0);
272 hisQ1D[ith]->SetDirectory(0);
273 hisQ3DROC[ith]->SetDirectory(0);
274 hisQ1DROC[ith]->SetDirectory(0);
275 hisQ1DROC[ith]->SetDirectory(0);
276 }
277 //
278 //
279 AliRawReader *reader = new AliRawReaderRoot(fileName);
280 reader->Reset();
281 AliAltroRawStream* stream = new AliAltroRawStream(reader);
282 stream->SelectRawData("TPC");
283 Int_t evtnr=0;
284 Int_t chunkNr=0;
285 //
286
287 while (reader->NextEvent()) {
288 Double_t shiftZ= gRandom->Rndm()*250;
289 //
290 if(evtnr>=maxEvents) {
291 chunkNr++;
292 evtnr=0;
293 pcstream->GetFile()->mkdir(Form("Chunk%d",chunkNr));
294 pcstream->GetFile()->cd(Form("Chunk%d",chunkNr));
295 for (Int_t ith=0; ith<3; ith++){
296 hisQ1D[ith]->Write(Form("His1DDrift_%d",ith));
297 hisQ2D[ith]->Write(Form("His2DDrift_%d",ith));
298 hisQ3D[ith]->Write(Form("His3DDrift_%d",ith));
299 hisQ1DROC[ith]->Write(Form("His1DROC_%d",ith));
300 hisQ2DROC[ith]->Write(Form("His2DROC_%d",ith));
301 hisQ3DROC[ith]->Write(Form("His3DROC_%d",ith));
5c70ec78 302 (*pcstream)<<"histo"<<
303 Form("hist1D_%d.=",ith*2+2)<<hisQ1D[ith]<<
304 Form("hist2D_%d.=",ith*2+2)<<hisQ2D[ith]<<
305 Form("hist3D_%d.=",ith*2+2)<<hisQ2D[ith]<<
306 Form("hist1DROC_%d.=",ith*2+2)<<hisQ1DROC[ith]<<
307 Form("hist2DROC_%d.=",ith*2+2)<<hisQ2DROC[ith]<<
308 Form("hist3DROC_%d.=",ith*2+2)<<hisQ2DROC[ith];
309 }
310 (*pcstream)<<"histo"<<"\n";
311 for (Int_t ith=0; ith<3; ith++){
c0172f82 312 hisQ1D[ith]->Reset();
313 hisQ2D[ith]->Reset();
314 hisQ3D[ith]->Reset();
315 hisQ1DROC[ith]->Reset();
316 hisQ2DROC[ith]->Reset();
317 hisQ3DROC[ith]->Reset();
318 }
319 }
320 cout<<"Chunk=\t"<<chunkNr<<"\tEvt=\t"<<evtnr<<endl;
321 evtnr++;
322 AliSysInfo::AddStamp(Form("Event%d",evtnr),evtnr);
323 AliTPCRawStreamV3 input(reader,(AliAltroMapping**)mapping);
324 //
325 while (input.NextDDL()){
326 Int_t sector = input.GetSector();
327 while ( input.NextChannel() ) {
328 Int_t row = input.GetRow();
329 Int_t rowAbs = row+ ((sector<36) ? 0: param->GetNRow(0));
330 Int_t pad = input.GetPad();
331 Int_t nPads = param->GetNPads(sector,row);
332 Double_t localX = param->GetPadRowRadii(sector,row);
333 Double_t localY = (pad-nPads)*param->GetPadPitchWidth(sector);
334 Double_t localPhi= TMath::ATan2(localY,localX);
335 Double_t phi = TMath::Pi()*((sector%18)+0.5)/9+localPhi;
336 if (sector%36>=18) phi+=TMath::TwoPi();
337 Double_t padLength=param->GetPadPitchLength(sector,row);
338 // Double_t padWidth=param->GetPadPitchWidth(sector);
339 // Double_t zLength=param->GetZLength();
340 Double_t deltaPhi=hisQ3D[0]->GetXaxis()->GetBinWidth(0);
341 Double_t deltaZ= hisQ3D[0]->GetZaxis()->GetBinWidth(0);
342 //
343 while ( input.NextBunch() ){
344 Int_t startTbin = (Int_t)input.GetStartTimeBin();
345 Int_t bunchlength = (Int_t)input.GetBunchLength();
346 const UShort_t *sig = input.GetSignals();
347 Int_t aboveTh[3]={0};
348 for (Int_t i=0; i<bunchlength; i++){
349 for (Int_t ith=0; ith<3; ith++){
350 if (sig[i]>(ith*2)+2) aboveTh[ith]++;
351 }
352 }
353 for (Int_t ith=0; ith<3; ith++){
354 if (aboveTh[ith%3]>1){
355 for (Int_t i=0; i<bunchlength; i++){
356 //
357 // normalization
358 //
359 Double_t volume3D = padLength*(localX*deltaPhi)*deltaZ;
360 Double_t zIonDrift = (param->GetZLength()-startTbin*param->GetZWidth());
361 zIonDrift+=shiftZ;
362 if (TMath::Abs(zIonDrift)<param->GetZLength()){
363 if ((sector%36)>=18) zIonDrift*=-1; // c side has opposite sign
5c70ec78 364 if (sector%36<18) hisQ1D[ith]->Fill(localX, sig[i]/volume3D);
365 hisQ2D[ith]->Fill(phi,localX, sig[i]/volume3D);
366 hisQ3D[ith]->Fill(phi,localX,zIonDrift,sig[i]/volume3D);
c0172f82 367 }
368 //
369 Double_t zIonROC = ((sector%36)<18)? shiftZ: -shiftZ; // z position of the "ion disc" - A side C side opposite sign
5c70ec78 370 if (sector%36<18) hisQ1DROC[ith]->Fill(localX, sig[i]/volume3D);
371 hisQ2DROC[ith]->Fill(phi,localX, sig[i]/volume3D);
372 hisQ3DROC[ith]->Fill(phi,localX,zIonROC,sig[i]/volume3D);
c0172f82 373 }
374 }
375 }
376 }
377 }
378 }
379 }
380 timer.Print();
381 delete pcstream;
382 return 0;
383}