FMD/AliFMDBackgroundCorrection.cxx

   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"
  45 #include "TH2F.h"
  46 #include "AliGenEventHeader.h"
  47 #include "AliHeader.h"
  48 #include "TFile.h"
  49 #include "TAxis.h"
  50 #include "AliCDBId.h"
  51 #include "AliCDBMetaData.h"
  52 #include "TROOT.h"
  53 #include "AliFMDParameters.h"
  54 #include "AliLog.h"
  55 #include "TList.h"
  56 #include "AliFMDAnaParameters.h"
  57 #include "AliFMDAnaCalibBackgroundCorrection.h"
  58
  59 ClassImp(AliFMDBackgroundCorrection)
  60 //_____________________________________________________________________
  61 AliFMDBackgroundCorrection::AliFMDBackgroundCorrection() :
  62   TNamed(),
  63   fCorrectionArray(),
  64   fPrimaryList()
  65 {}
  66
  67 //_____________________________________________________________________
  68 AliFMDBackgroundCorrection::AliFMDInputBG::AliFMDInputBG() :
  69   AliFMDInput(),
  70   fPrimaryArray(),
  71   fHitArray(),
  72   fPrim(0),
  73   fHits(0),
  74   fZvtxCut(0),
  75   fNvtxBins(0),
  76   fPrevTrack(-1),
  77   fPrevDetector(-1),
  78   fPrevRing('Q'),
  79   fNbinsEta(100)
  80 {
  81   AddLoad(kTracks);
  82   AddLoad(kHits);
  83   AddLoad(kKinematics);
  84   AddLoad(kHeader);
  85 }
  86
  87 //_____________________________________________________________________
  88
  89 void AliFMDBackgroundCorrection::GenerateBackgroundCorrection(const Int_t nvtxbins,
  90                                                               Float_t zvtxcut,
  91                                                               const Int_t nBinsEta,
  92                                                               Bool_t storeInAlien,
  93                                                               Int_t runNo,
  94                                                               Int_t endRunNo,
  95                                                               const Char_t* filename,
  96                                                               Bool_t simulate,
  97                                                               Int_t nEvents) {
  98
  99   //TGrid::Connect("alien:",0,0,"t");
 100   if(simulate)
 101     Simulate(nEvents);
 102   else {
 103     //AliCDBManager::Instance()->SetDefaultStorage("alien://Folder=/alice/data/2008/LHC08d/OCDB/");
 104     AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
 105     AliCDBManager::Instance()->SetRun(runNo);
 106
 107 #if defined(__CINT__)
 108     gSystem->Load("liblhapdf");
 109     gSystem->Load("libEGPythia6");
 110     gSystem->Load("libpythia6");
 111     gSystem->Load("libAliPythia6");
 112     gSystem->Load("libgeant321");
 113 #endif
 114       //
 115   }
 116
 117   //Setting up the geometry
 118   //-----------------------------------------------
 119   if (AliGeomManager::GetGeometry() == NULL)
 120     AliGeomManager::LoadGeometry();
 121
 122   AliFMDGeometry* geo = AliFMDGeometry::Instance();
 123   geo->Init();
 124   geo->InitTransformations();
 125
 126   AliInfo("Processing hits and primaries ");
 127
 128   AliFMDInputBG input;
 129   input.SetVtxCutZ(zvtxcut);
 130   input.SetNvtxBins(nvtxbins);
 131   input.SetNbinsEta(nBinsEta);
 132   input.Run();
 133
 134   AliInfo(Form("Found %d primaries and %d hits.", input.GetNprim(),input.GetNhits()));
 135
 136   TObjArray* hitArray = input.GetHits();
 137   TObjArray* primaryArray = input.GetPrimaries();
 138   fCorrectionArray.SetName("FMD_bg_correction");
 139   fCorrectionArray.SetOwner();
 140
 141   TList* primaryList     = new TList();
 142   primaryList->SetName("primaries");
 143
 144   TList* hitList     = new TList();
 145   hitList->SetName("hits");
 146   TList* corrList    = new TList();
 147   corrList->SetName("corrections");
 148
 149   AliFMDAnaCalibBackgroundCorrection* background = new AliFMDAnaCalibBackgroundCorrection();
 150
 151   for(Int_t det= 1; det <=3; det++) {
 152     Int_t nRings = (det==1 ? 1 : 2);
 153
 154     TObjArray* detArrayCorrection = new TObjArray();
 155     detArrayCorrection->SetName(Form("FMD%d",det));
 156     fCorrectionArray.AddAtAndExpand(detArrayCorrection,det);
 157
 158
 159     for(Int_t iring = 0; iring<nRings; iring++) {
 160       TObjArray* primRingArray = (TObjArray*)primaryArray->At(iring);
 161       Char_t ringChar = (iring == 0 ? 'I' : 'O');
 162       TObjArray* vtxArrayCorrection = new TObjArray();
 163       vtxArrayCorrection->SetName(Form("FMD%d%c",det,ringChar));
 164       detArrayCorrection->AddAtAndExpand(vtxArrayCorrection,iring);
 165
 166       for(Int_t vertexBin=0;vertexBin<nvtxbins;vertexBin++) {
 167         TObjArray* detArray  = (TObjArray*)hitArray->At(det);
 168         TObjArray* vtxArray  = (TObjArray*)detArray->At(iring);
 169         TH2F* hHits          = (TH2F*)vtxArray->At(vertexBin);
 170         hitList->Add(hHits);
 171         TH2F* hPrimary  = (TH2F*)primRingArray->At(vertexBin);
 172         primaryList->Add(hPrimary);
 173         TH2F* hCorrection = (TH2F*)hHits->Clone(Form("FMD%d%c_vtxbin_%d_correction",det,ringChar,vertexBin));
 174         hCorrection->Divide(hPrimary);
 175         corrList->Add(hCorrection);
 176         hCorrection->SetTitle(hCorrection->GetName());
 177         vtxArrayCorrection->AddAtAndExpand(hCorrection,vertexBin);
 178         background->SetBgCorrection(det,ringChar,vertexBin,hCorrection);
 179       }
 180
 181     }
 182   }
 183
 184   TAxis refAxis(nvtxbins,-1*zvtxcut,zvtxcut);
 185
 186   background->SetRefAxis(&refAxis);
 187
 188   TFile*  fout = new TFile(filename,"RECREATE");
 189   refAxis.Write("vertexbins");
 190
 191   hitList->Write();
 192   primaryList->Write();
 193   corrList->Write();
 194
 195   TObjArray* container = new TObjArray();
 196   container->SetOwner();
 197   container->AddAtAndExpand(&refAxis,0);
 198   container->AddAtAndExpand(&fCorrectionArray,1);
 199   container->AddAtAndExpand(hitArray,2);
 200   container->AddAtAndExpand(primaryArray,3);
 201
 202
 203   if(storeInAlien) {
 204     AliCDBManager* cdb = AliCDBManager::Instance();
 205     cdb->SetDefaultStorage("local://$ALICE_ROOT");
 206     AliCDBId      id(AliFMDAnaParameters::GetBackgroundPath(),runNo,endRunNo);
 207
 208     AliCDBMetaData* meta = new AliCDBMetaData;
 209     meta->SetResponsible(gSystem->GetUserInfo()->fRealName.Data());
 210     meta->SetAliRootVersion(gROOT->GetVersion());
 211     meta->SetBeamPeriod(1);
 212     meta->SetComment("Background Correction for FMD");
 213     meta->SetProperty("key1", background );
 214     cdb->Put(background, id, meta);
 215
 216   }
 217
 218   fout->Close();
 219
 220
 221   }
 222
 223 //_____________________________________________________________________
 224 void AliFMDBackgroundCorrection::Simulate(Int_t nEvents) {
 225
 226   AliSimulation sim ;
 227   sim.SetRunNumber(0);
 228   TGrid::Connect("alien:",0,0,"t");
 229   sim.SetDefaultStorage("alien://Folder=/alice/data/2008/LHC08d/OCDB/");
 230   sim.SetConfigFile("Config.C");
 231   sim.SetRunQA("FMD:");
 232   TStopwatch timer;
 233   timer.Start();
 234   sim.RunSimulation(nEvents);
 235   timer.Stop();
 236   timer.Print();
 237
 238 }
 239
 240 //_____________________________________________________________________
 241 Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::ProcessHit( AliFMDHit* h, TParticle* p) {
 242
 243   AliGenEventHeader* genHeader = fLoader->GetHeader()->GenEventHeader();
 244   TArrayF vertex;
 245   genHeader->PrimaryVertex(vertex);
 246
 247   if(TMath::Abs(vertex.At(2)) > fZvtxCut)
 248     return kTRUE;
 249
 250   Double_t delta = 2*fZvtxCut/fNvtxBins;
 251   Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
 252   Int_t vertexBin = (Int_t)vertexBinDouble;
 253
 254   Int_t i = h->Track();
 255
 256   if(h)
 257     if(h->Q() !=  0 && (i != fPrevTrack || h->Detector() != fPrevDetector || h->Ring() != fPrevRing)) {
 258
 259       Int_t det = h->Detector();
 260       Char_t ring = h->Ring();
 261       Int_t iring = (ring == 'I' ? 0 : 1);
 262
 263       TObjArray* detArray  = (TObjArray*)fHitArray.At(det);
 264       TObjArray* vtxArray  = (TObjArray*)detArray->At(iring);
 265       TH2F* hHits          = (TH2F*)vtxArray->At(vertexBin);
 266
 267       Float_t phi   = TMath::ATan2(h->Py(),h->Px());
 268       if(phi<0)
 269         phi = phi+2*TMath::Pi();
 270       Float_t theta = TMath::ATan2(TMath::Sqrt(TMath::Power(h->X(),2)+TMath::Power(h->Y(),2)),h->Z()+vertex.At(2));
 271       Float_t eta   = -1*TMath::Log(TMath::Tan(0.5*theta));
 272       hHits->Fill(eta,phi);
 273       fHits++;
 274       fPrevDetector = det;
 275       fPrevRing     = ring;
 276     }
 277
 278   fPrevTrack = i;
 279   return kTRUE;
 280 }
 281
 282 //_____________________________________________________________________
 283 Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Init() {
 284
 285   fPrimaryArray.SetOwner();
 286   fPrimaryArray.SetName("FMD_primary");
 287
 288   for(Int_t iring = 0; iring<2;iring++) {
 289     Char_t ringChar = (iring == 0 ? 'I' : 'O');
 290     TObjArray* ringArray = new TObjArray();
 291     ringArray->SetName(Form("FMD_%c",ringChar));
 292     fPrimaryArray.AddAtAndExpand(ringArray,iring);
 293     Int_t nSec = (iring == 1 ? 40 : 20);
 294     for(Int_t v=0; v<fNvtxBins;v++) {
 295
 296       TH2F* hPrimary       = new TH2F(Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
 297                                       Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
 298                                       fNbinsEta, -6,6, nSec, 0,2*TMath::Pi());
 299       hPrimary->Sumw2();
 300       ringArray->AddAtAndExpand(hPrimary,v);
 301     }
 302   }
 303
 304
 305   fHitArray.SetOwner();
 306   fHitArray.SetName("FMD_hits");
 307
 308   for(Int_t det =1; det<=3;det++)
 309     {
 310       TObjArray* detArrayHits = new TObjArray();
 311       detArrayHits->SetName(Form("FMD%d",det));
 312       fHitArray.AddAtAndExpand(detArrayHits,det);
 313       Int_t nRings = (det==1 ? 1 : 2);
 314       for(Int_t ring = 0;ring<nRings;ring++)
 315         {
 316           Int_t nSec = (ring == 1 ? 40 : 20);
 317           Char_t ringChar = (ring == 0 ? 'I' : 'O');
 318           TObjArray* vtxArrayHits = new TObjArray();
 319           vtxArrayHits->SetName(Form("FMD%d%c",det,ringChar));
 320           detArrayHits->AddAtAndExpand(vtxArrayHits,ring);
 321           for(Int_t v=0; v<fNvtxBins;v++)
 322             {
 323
 324               TH2F* hHits          = new TH2F(Form("hHits_FMD%d%c_vtx%d",det,ringChar,v),
 325                                               Form("hHits_FMD%d%c_vtx%d",det,ringChar,v),
 326                                               fNbinsEta, -6,6, nSec, 0, 2*TMath::Pi());
 327               hHits->Sumw2();
 328               vtxArrayHits->AddAtAndExpand(hHits,v);
 329
 330             }
 331         }
 332
 333     }
 334
 335
 336
 337   AliFMDInput::Init();
 338
 339
 340   return kTRUE;
 341 }
 342 //_____________________________________________________________________
 343
 344 Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Begin(Int_t event ) {
 345
 346   Bool_t retVal = AliFMDInput::Begin(event);
 347
 348   AliStack* partStack=fLoader->Stack();
 349
 350   Int_t nTracks=partStack->GetNtrack();
 351   AliGenEventHeader* genHeader = fLoader->GetHeader()->GenEventHeader();
 352   TArrayF vertex;
 353   genHeader->PrimaryVertex(vertex);
 354
 355   if(TMath::Abs(vertex.At(2)) > fZvtxCut)
 356     return kTRUE;
 357
 358   Double_t delta = 2*fZvtxCut/fNvtxBins;
 359   Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
 360   Int_t vertexBin = (Int_t)vertexBinDouble;
 361
 362   for(Int_t j=0;j<nTracks;j++)
 363         {
 364           TParticle* p  = partStack->Particle(j);
 365           TDatabasePDG* pdgDB = TDatabasePDG::Instance();
 366           TParticlePDG* pdgPart = pdgDB->GetParticle(p->GetPdgCode());
 367           Float_t charge = (pdgPart ? pdgPart->Charge() : 0);
 368           Float_t phi = TMath::ATan2(p->Py(),p->Px());
 369
 370           if(phi<0)
 371             phi = phi+2*TMath::Pi();
 372           if(p->Theta() == 0) continue;
 373               Float_t eta   = -1*TMath::Log(TMath::Tan(0.5*p->Theta()));
 374
 375           Bool_t primary = (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);
 376           if(primary) {
 377             fPrim++;
 378             TObjArray* innerArray = (TObjArray*)fPrimaryArray.At(0);
 379             TObjArray* outerArray = (TObjArray*)fPrimaryArray.At(1);
 380
 381             TH2F* hPrimaryInner  = (TH2F*)innerArray->At(vertexBin);
 382             TH2F* hPrimaryOuter  = (TH2F*)outerArray->At(vertexBin);
 383             hPrimaryInner->Fill(eta,phi);
 384             hPrimaryOuter->Fill(eta,phi);
 385           }
 386         }
 387
 388   return retVal;
 389 }
 390 //_____________________________________________________________________
 391 //
 392 // EOF
 393 //
 394
 395