Extacting the OCDB in a separate module. The detectors have write permission in the...
[u/mrichter/AliRoot.git] / FMD / AliFMDBackgroundCorrection.cxx
CommitLineData
df0d5480 1/**************************************************************************
2 * Copyright(c) 2004, 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// Thil class computes background corrections for the FMD. The correction is computed
17// in eta,phi cells and the objects stored can be put into alien to use with analysis.
18// It is based on the AliFMDInput class that is used to loop over hits and primaries.
19//
20// Author: Hans Hjersing Dalsgaard, NBI, hans.dalsgaard@cern.ch
21//
22//
23
24#include "AliSimulation.h"
25#include "TStopwatch.h"
26#include "iostream"
27#include "TGrid.h"
28#include "AliRunLoader.h"
29#include "AliGeomManager.h"
30#include "AliFMDGeometry.h"
31#include "AliStack.h"
32#include "TParticle.h"
33#include "TDatabasePDG.h"
34#include "TParticlePDG.h"
35#include "AliRun.h"
36#include "AliFMDBackgroundCorrection.h"
37#include "TSystem.h"
38#include "AliCDBManager.h"
39#include "TTree.h"
40#include "TClonesArray.h"
41#include "TBranch.h"
42#include "AliFMDHit.h"
43#include "AliLoader.h"
44#include "AliFMD.h"
df0d5480 45#include "AliGenEventHeader.h"
46#include "AliHeader.h"
47#include "TFile.h"
48#include "TAxis.h"
49#include "AliCDBId.h"
50#include "AliCDBMetaData.h"
51#include "TROOT.h"
52#include "AliFMDParameters.h"
53#include "AliLog.h"
26c8b5ff 54#include "TList.h"
8dc823cc 55#include "AliFMDAnaParameters.h"
b82e76e0 56#include "AliFMDAnaCalibBackgroundCorrection.h"
08d168d9 57#include "AliTrackReference.h"
58#include "AliFMDStripIndex.h"
df0d5480 59
60ClassImp(AliFMDBackgroundCorrection)
61//_____________________________________________________________________
62AliFMDBackgroundCorrection::AliFMDBackgroundCorrection() :
63 TNamed(),
ca3af5d9 64 fCorrectionArray(),
65 fPrimaryList()
df0d5480 66{}
67
68//_____________________________________________________________________
08d168d9 69AliFMDBackgroundCorrection::AliFMDInputBG::AliFMDInputBG(Bool_t hits_not_trackref) :
5d5bc6be 70 AliFMDInput("galice.root"),
8dc823cc 71 fPrimaryArray(),
72 fHitArray(),
627ba74b 73 fHitMap(),
df0d5480 74 fPrim(0),
75 fHits(0),
76 fZvtxCut(0),
77 fNvtxBins(0),
78 fPrevTrack(-1),
79 fPrevDetector(-1),
80 fPrevRing('Q'),
5d5bc6be 81 fPrevSec(-1),
df0d5480 82 fNbinsEta(100)
83{
84 AddLoad(kTracks);
08d168d9 85 if(hits_not_trackref)
86 AddLoad(kHits);
87 else
88 AddLoad(kTrackRefs);
df0d5480 89 AddLoad(kKinematics);
90 AddLoad(kHeader);
627ba74b 91
92
df0d5480 93}
94
95//_____________________________________________________________________
96
08d168d9 97void AliFMDBackgroundCorrection::GenerateBackgroundCorrection(Bool_t from_hits,
98 const Int_t nvtxbins,
df0d5480 99 Float_t zvtxcut,
ca3af5d9 100 const Int_t nBinsEta,
df0d5480 101 Bool_t storeInAlien,
102 Int_t runNo,
103 Int_t endRunNo,
104 const Char_t* filename,
105 Bool_t simulate,
109e97cd 106 Int_t nEvents,
107 Bool_t inFile,
108 const Char_t* infilename) {
df0d5480 109
ca3af5d9 110 //TGrid::Connect("alien:",0,0,"t");
df0d5480 111 if(simulate)
112 Simulate(nEvents);
113 else {
ca3af5d9 114 //AliCDBManager::Instance()->SetDefaultStorage("alien://Folder=/alice/data/2008/LHC08d/OCDB/");
162637e4 115 AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
df0d5480 116 AliCDBManager::Instance()->SetRun(runNo);
117
118#if defined(__CINT__)
119 gSystem->Load("liblhapdf");
120 gSystem->Load("libEGPythia6");
121 gSystem->Load("libpythia6");
122 gSystem->Load("libAliPythia6");
123 gSystem->Load("libgeant321");
124#endif
ca3af5d9 125 //
df0d5480 126 }
127
128 //Setting up the geometry
129 //-----------------------------------------------
130 if (AliGeomManager::GetGeometry() == NULL)
131 AliGeomManager::LoadGeometry();
132
133 AliFMDGeometry* geo = AliFMDGeometry::Instance();
134 geo->Init();
135 geo->InitTransformations();
ca3af5d9 136
df0d5480 137 AliInfo("Processing hits and primaries ");
138
08d168d9 139 AliFMDInputBG input(from_hits);
df0d5480 140
109e97cd 141 if(!inFile) {
142
143 input.SetVtxCutZ(zvtxcut);
144 input.SetNvtxBins(nvtxbins);
145 input.SetNbinsEta(nBinsEta);
146 input.Run();
147 }
df0d5480 148
109e97cd 149 AliInfo(Form("Found %d primaries and %d hits.", input.GetNprim(),input.GetNhits()));
150 TObjArray* hitArray ;
151 TObjArray* primaryArray;
152 if(inFile) {
153 TFile* infile = TFile::Open(infilename);
154 hitArray = new TObjArray();
155 primaryArray = new TObjArray();
156
157 for(Int_t det =1; det<=3;det++)
158 {
159 TObjArray* detArrayHits = new TObjArray();
160 detArrayHits->SetName(Form("FMD%d",det));
161 hitArray->AddAtAndExpand(detArrayHits,det);
162 Int_t nRings = (det==1 ? 1 : 2);
163 for(Int_t ring = 0;ring<nRings;ring++)
164 {
165
166 Char_t ringChar = (ring == 0 ? 'I' : 'O');
167 TObjArray* vtxArrayHits = new TObjArray();
168 vtxArrayHits->SetName(Form("FMD%d%c",det,ringChar));
169 detArrayHits->AddAtAndExpand(vtxArrayHits,ring);
170 for(Int_t v=0; v<nvtxbins;v++)
171 {
172
173 TH2F* hHits = (TH2F*)infile->Get(Form("hHits_FMD%d%c_vtx%d",det,ringChar,v));
174
175
176 vtxArrayHits->AddAtAndExpand(hHits,v);
177
178 }
179 }
180
181 }
182
183 for(Int_t iring = 0; iring<2;iring++) {
184 Char_t ringChar = (iring == 0 ? 'I' : 'O');
185 TObjArray* ringArray = new TObjArray();
186 ringArray->SetName(Form("FMD_%c",ringChar));
187 primaryArray->AddAtAndExpand(ringArray,iring);
188 for(Int_t v=0; v<nvtxbins;v++) {
189
190 TH2F* hPrimary = (TH2F*)infile->Get(Form("hPrimary_FMD_%c_vtx%d",ringChar,v));
191 ringArray->AddAtAndExpand(hPrimary,v);
192 }
193 }
194
195
196 }
197 else {
198 hitArray = input.GetHits();
199 primaryArray = input.GetPrimaries();
200 }
df0d5480 201 fCorrectionArray.SetName("FMD_bg_correction");
202 fCorrectionArray.SetOwner();
ca3af5d9 203
204 TList* primaryList = new TList();
26c8b5ff 205 primaryList->SetName("primaries");
ca3af5d9 206
26c8b5ff 207 TList* hitList = new TList();
208 hitList->SetName("hits");
209 TList* corrList = new TList();
210 corrList->SetName("corrections");
109e97cd 211
b82e76e0 212 AliFMDAnaCalibBackgroundCorrection* background = new AliFMDAnaCalibBackgroundCorrection();
213
df0d5480 214 for(Int_t det= 1; det <=3; det++) {
215 Int_t nRings = (det==1 ? 1 : 2);
216
217 TObjArray* detArrayCorrection = new TObjArray();
218 detArrayCorrection->SetName(Form("FMD%d",det));
219 fCorrectionArray.AddAtAndExpand(detArrayCorrection,det);
220
221
222 for(Int_t iring = 0; iring<nRings; iring++) {
223 TObjArray* primRingArray = (TObjArray*)primaryArray->At(iring);
224 Char_t ringChar = (iring == 0 ? 'I' : 'O');
225 TObjArray* vtxArrayCorrection = new TObjArray();
226 vtxArrayCorrection->SetName(Form("FMD%d%c",det,ringChar));
227 detArrayCorrection->AddAtAndExpand(vtxArrayCorrection,iring);
228
229 for(Int_t vertexBin=0;vertexBin<nvtxbins;vertexBin++) {
230 TObjArray* detArray = (TObjArray*)hitArray->At(det);
231 TObjArray* vtxArray = (TObjArray*)detArray->At(iring);
232 TH2F* hHits = (TH2F*)vtxArray->At(vertexBin);
26c8b5ff 233 hitList->Add(hHits);
df0d5480 234 TH2F* hPrimary = (TH2F*)primRingArray->At(vertexBin);
26c8b5ff 235 primaryList->Add(hPrimary);
df0d5480 236 TH2F* hCorrection = (TH2F*)hHits->Clone(Form("FMD%d%c_vtxbin_%d_correction",det,ringChar,vertexBin));
ca3af5d9 237 hCorrection->Divide(hPrimary);
26c8b5ff 238 corrList->Add(hCorrection);
8dc823cc 239 hCorrection->SetTitle(hCorrection->GetName());
df0d5480 240 vtxArrayCorrection->AddAtAndExpand(hCorrection,vertexBin);
b82e76e0 241 background->SetBgCorrection(det,ringChar,vertexBin,hCorrection);
df0d5480 242 }
243
244 }
245 }
246
247 TAxis refAxis(nvtxbins,-1*zvtxcut,zvtxcut);
109e97cd 248 if(inFile) {
249 TFile* infile = TFile::Open(infilename);
250 TAxis* refaxis = (TAxis*)infile->Get("vertexbins");
251 background->SetRefAxis(refaxis);
252
253 }
254 else background->SetRefAxis(&refAxis);
df0d5480 255
8dc823cc 256 TFile* fout = new TFile(filename,"RECREATE");
df0d5480 257 refAxis.Write("vertexbins");
26c8b5ff 258
259 hitList->Write();
260 primaryList->Write();
261 corrList->Write();
ca3af5d9 262
8dc823cc 263 TObjArray* container = new TObjArray();
264 container->SetOwner();
265 container->AddAtAndExpand(&refAxis,0);
266 container->AddAtAndExpand(&fCorrectionArray,1);
267 container->AddAtAndExpand(hitArray,2);
268 container->AddAtAndExpand(primaryArray,3);
b82e76e0 269
270
df0d5480 271 if(storeInAlien) {
272 AliCDBManager* cdb = AliCDBManager::Instance();
162637e4 273 cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
8dc823cc 274 AliCDBId id(AliFMDAnaParameters::GetBackgroundPath(),runNo,endRunNo);
df0d5480 275
276 AliCDBMetaData* meta = new AliCDBMetaData;
277 meta->SetResponsible(gSystem->GetUserInfo()->fRealName.Data());
278 meta->SetAliRootVersion(gROOT->GetVersion());
279 meta->SetBeamPeriod(1);
280 meta->SetComment("Background Correction for FMD");
b82e76e0 281 meta->SetProperty("key1", background );
282 cdb->Put(background, id, meta);
df0d5480 283
284 }
285
8dc823cc 286 fout->Close();
287
df0d5480 288
ca3af5d9 289 }
290
df0d5480 291//_____________________________________________________________________
292void AliFMDBackgroundCorrection::Simulate(Int_t nEvents) {
293
294 AliSimulation sim ;
295 sim.SetRunNumber(0);
296 TGrid::Connect("alien:",0,0,"t");
297 sim.SetDefaultStorage("alien://Folder=/alice/data/2008/LHC08d/OCDB/");
298 sim.SetConfigFile("Config.C");
299 sim.SetRunQA("FMD:");
300 TStopwatch timer;
301 timer.Start();
302 sim.RunSimulation(nEvents);
303 timer.Stop();
304 timer.Print();
305
306}
307
308//_____________________________________________________________________
faf80567 309Bool_t
310AliFMDBackgroundCorrection::AliFMDInputBG::ProcessHit(AliFMDHit* h,
311 TParticle* /*p*/)
312{
df0d5480 313
08d168d9 314 if(!h)
315 return kTRUE;
faf80567 316 Bool_t retval = ProcessEvent(h->Detector(),
317 h->Ring(),
98accb9e 318 h->Sector(),
319 h->Strip(),
faf80567 320 h->Track(),
321 h->Q());
08d168d9 322
323 return retval;
324}
325//_____________________________________________________________________
faf80567 326Bool_t
327AliFMDBackgroundCorrection::AliFMDInputBG::ProcessTrackRef(AliTrackReference* tr,
328 TParticle* p)
329{
08d168d9 330 if(!tr)
331 return kTRUE;
332 UShort_t det,sec,strip;
333 Char_t ring;
334 AliFMDStripIndex::Unpack(tr->UserId(),det,ring,sec,strip);
faf80567 335 Int_t nTrack = tr->GetTrack();
336 TDatabasePDG* pdgDB = TDatabasePDG::Instance();
08d168d9 337 TParticlePDG* pdgPart = pdgDB->GetParticle(p->GetPdgCode());
faf80567 338 Float_t charge = (pdgPart ? pdgPart->Charge() : 0);
339 Bool_t retval = ProcessEvent(det,ring,sec,strip,nTrack,charge);
08d168d9 340 return retval;
341
342}
343//_____________________________________________________________________
faf80567 344Bool_t
345AliFMDBackgroundCorrection::AliFMDInputBG::ProcessEvent(UShort_t det,
346 Char_t ring,
347 UShort_t sec,
348 UShort_t strip,
349 Int_t nTrack,
350 Float_t charge)
08d168d9 351{
352
df0d5480 353
df0d5480 354
faf80567 355 if(charge != 0 &&
627ba74b 356 ((nTrack != fPrevTrack) ||
357 (det != fPrevDetector) ||
358 (ring != fPrevRing) ||
359 (sec != fPrevSec))) {
360 fHitMap.operator()(det,ring,sec,strip) = 1;
08d168d9 361 fHits++;
5d5bc6be 362
08d168d9 363 }
ca3af5d9 364
5d5bc6be 365 fPrevDetector = det;
366 fPrevRing = ring;
367 fPrevSec = sec;
368 fPrevTrack = nTrack;
369
df0d5480 370 return kTRUE;
08d168d9 371
df0d5480 372}
373
374//_____________________________________________________________________
faf80567 375Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Init()
376{
df0d5480 377 fPrimaryArray.SetOwner();
378 fPrimaryArray.SetName("FMD_primary");
379
627ba74b 380 fPrimaryMapInner.SetBins(fNbinsEta, -6,6, 20, 0, 2*TMath::Pi());
381 fPrimaryMapOuter.SetBins(fNbinsEta, -6,6, 40, 0, 2*TMath::Pi());
382 fPrimaryMapInner.SetName("fPrimaryMapInner");
383 fPrimaryMapInner.SetName("fPrimaryMapOuter");
384
385 fPrimaryMapInner.Sumw2();
386 fPrimaryMapOuter.Sumw2();
df0d5480 387 for(Int_t iring = 0; iring<2;iring++) {
388 Char_t ringChar = (iring == 0 ? 'I' : 'O');
389 TObjArray* ringArray = new TObjArray();
390 ringArray->SetName(Form("FMD_%c",ringChar));
391 fPrimaryArray.AddAtAndExpand(ringArray,iring);
392 Int_t nSec = (iring == 1 ? 40 : 20);
393 for(Int_t v=0; v<fNvtxBins;v++) {
394
395 TH2F* hPrimary = new TH2F(Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
396 Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
397 fNbinsEta, -6,6, nSec, 0,2*TMath::Pi());
b82e76e0 398 hPrimary->Sumw2();
df0d5480 399 ringArray->AddAtAndExpand(hPrimary,v);
400 }
401 }
402
403
404 fHitArray.SetOwner();
405 fHitArray.SetName("FMD_hits");
406
faf80567 407 for(Int_t det =1; det<=3;det++) {
408 TObjArray* detArrayHits = new TObjArray();
409 detArrayHits->SetName(Form("FMD%d",det));
410 fHitArray.AddAtAndExpand(detArrayHits,det);
411 Int_t nRings = (det==1 ? 1 : 2);
412 for(Int_t ring = 0;ring<nRings;ring++) {
413 Int_t nSec = (ring == 1 ? 40 : 20);
414 Char_t ringChar = (ring == 0 ? 'I' : 'O');
415 TObjArray* vtxArrayHits = new TObjArray();
416 vtxArrayHits->SetName(Form("FMD%d%c",det,ringChar));
417 detArrayHits->AddAtAndExpand(vtxArrayHits,ring);
418 for(Int_t v=0; v<fNvtxBins;v++) {
419 TH2F* hHits = new TH2F(Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
420 Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
421 fNbinsEta, -6,6, nSec, 0, 2*TMath::Pi());
422 hHits->Sumw2();
423 vtxArrayHits->AddAtAndExpand(hHits,v);
424
425 }
df0d5480 426 }
faf80567 427 }
df0d5480 428
df0d5480 429 AliFMDInput::Init();
430
df0d5480 431 return kTRUE;
432}
ca3af5d9 433//_____________________________________________________________________
434
faf80567 435Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Begin(Int_t event )
436{
df0d5480 437
faf80567 438 Bool_t retVal = AliFMDInput::Begin(event);
439 AliStack* partStack = fLoader->Stack();
440 Int_t nTracks = partStack->GetNtrack();
ca3af5d9 441
faf80567 442
443 for(Int_t j=0;j<nTracks;j++) {
5d5bc6be 444 TParticle* p = partStack->Particle(j);
445 TDatabasePDG* pdgDB = TDatabasePDG::Instance();
446 TParticlePDG* pdgPart = pdgDB->GetParticle(p->GetPdgCode());
447 Float_t charge = (pdgPart ? pdgPart->Charge() : 0);
448 Float_t phi = TMath::ATan2(p->Py(),p->Px());
faf80567 449
450 if(phi<0) phi = phi+2*TMath::Pi();
627ba74b 451 // if(p->Theta() == 0) continue;
452 Float_t eta = p->Eta();
453
454 // std::cout<<-1*TMath::Log(TMath::Tan(0.5*p->Theta()))<<std::endl;
faf80567 455
456 Bool_t primary = partStack->IsPhysicalPrimary(j);
457 //(charge!=0)&&(TMath::Abs(p->Vx() - vertex.At(0))<0.01)&&(TMath::Abs(p->Vy() - vertex.At(1))<0.01)&&(TMath::Abs(p->Vz() - vertex.At(2))<0.01);
458 if(primary && charge !=0) {
459
627ba74b 460
461
faf80567 462 fPrim++;
faf80567 463
627ba74b 464 fPrimaryMapInner.Fill(eta,phi);
465 fPrimaryMapOuter.Fill(eta,phi);
faf80567 466 }
467 }
109e97cd 468
ca3af5d9 469 return retVal;
470}
471//_____________________________________________________________________
627ba74b 472Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::End() {
473
474 Bool_t retval = AliFMDInput::End();
475
476 AliGenEventHeader* genHeader = fLoader->GetHeader()->GenEventHeader();
477 TArrayF vertex;
478 genHeader->PrimaryVertex(vertex);
479
480 if(TMath::Abs(vertex.At(2)) > fZvtxCut)
481 return kTRUE;
482
483 Double_t delta = 2*fZvtxCut/fNvtxBins;
484 Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
485 Int_t vertexBin = (Int_t)vertexBinDouble;
486 //Primaries
487 TObjArray* innerArray = (TObjArray*)fPrimaryArray.At(0);
488 TObjArray* outerArray = (TObjArray*)fPrimaryArray.At(1);
489
490 TH2F* hPrimaryInner = (TH2F*)innerArray->At(vertexBin);
491 TH2F* hPrimaryOuter = (TH2F*)outerArray->At(vertexBin);
492
493 hPrimaryInner->Add(&fPrimaryMapInner);
494 hPrimaryOuter->Add(&fPrimaryMapOuter);
495
496 //Hits
497 for(UShort_t det=1;det<=3;det++) {
498 Int_t nRings = (det==1 ? 1 : 2);
499 for (UShort_t ir = 0; ir < nRings; ir++) {
500 Char_t ring = (ir == 0 ? 'I' : 'O');
501 UShort_t nsec = (ir == 0 ? 20 : 40);
502 UShort_t nstr = (ir == 0 ? 512 : 256);
503
504 for(UShort_t sec =0; sec < nsec; sec++) {
505
506 for(UShort_t strip = 0; strip < nstr; strip++) {
507
508 if(fHitMap.operator()(det,ring,sec,strip) > 0) {
509
510 Double_t x,y,z;
511 AliFMDGeometry* fmdgeo = AliFMDGeometry::Instance();
512 fmdgeo->Detector2XYZ(det,ring,sec,strip,x,y,z);
513
514 Int_t iring = (ring == 'I' ? 0 : 1);
515
516 TObjArray* detArray = (TObjArray*)fHitArray.At(det);
517 TObjArray* vtxArray = (TObjArray*)detArray->At(iring);
518 TH2F* hHits = (TH2F*)vtxArray->At(vertexBin);
519
520 Float_t phi = TMath::ATan2(y,x);
521 if(phi<0) phi = phi+2*TMath::Pi();
522 Float_t r = TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2));
523 Float_t theta = TMath::ATan2(r,z-vertex.At(2));
524 Float_t eta = -1*TMath::Log(TMath::Tan(0.5*theta));
525 hHits->Fill(eta,phi);
526 }
527
528 }
529 }
530 }
531 }
532
533 fPrimaryMapInner.Reset();
534 fPrimaryMapOuter.Reset();
535 fHitMap.Reset(0);
536
537 return retval;
538}
539
ca3af5d9 540//
541// EOF
542//
df0d5480 543
544