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