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