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