PWG2/FORWARD/analysis/AliFMDAnalysisTaskGenerateBackground.cxx

   1 #include "AliFMDAnalysisTaskGenerateBackground.h"
   2 #include "AliESDEvent.h"
   3 #include "iostream"
   4 #include "AliESDFMD.h"
   5 #include "TH2F.h"
   6 #include "AliTrackReference.h"
   7 #include "AliStack.h"
   8 #include "AliFMDAnaParameters.h"
   9 #include "AliFMDStripIndex.h"
  10 #include "AliStack.h"
  11 #include "AliMCParticle.h"
  12 #include "AliMCEvent.h"
  13 //#include "AliFMDGeometry.h"
  14 #include "TArray.h"
  15 #include "AliGenEventHeader.h"
  16 #include "AliHeader.h"
  17 #include "AliFMDAnaCalibBackgroundCorrection.h"
  18 //#include "AliCDBManager.h"
  19 //#include "AliCDBId.h"
  20 //#include "AliCDBMetaData.h"
  21 #include "TSystem.h"
  22 #include "TROOT.h"
  23 #include "TAxis.h"
  24 ClassImp(AliFMDAnalysisTaskGenerateBackground)
  25
  26 //_____________________________________________________________________
  27 AliFMDAnalysisTaskGenerateBackground::AliFMDAnalysisTaskGenerateBackground():
  28 AliAnalysisTaskSE(),
  29   fListOfHits(),
  30   fListOfPrimaries(),
  31   fListOfCorrection(),
  32   fVertexBins(),
  33   fLastTrackByStrip(0),
  34   fHitsByStrip(0),
  35   fZvtxCut(10),
  36   fNvtxBins(10),
  37   fNbinsEta(200),
  38   fBackground(0)
  39 {
  40   // Default constructor
  41 }
  42 //_____________________________________________________________________
  43 AliFMDAnalysisTaskGenerateBackground::AliFMDAnalysisTaskGenerateBackground(const char* name):
  44   AliAnalysisTaskSE(name),
  45   fListOfHits(),
  46   fListOfPrimaries(),
  47   fListOfCorrection(),
  48   fVertexBins(),
  49   fLastTrackByStrip(0),
  50   fHitsByStrip(0),
  51   fZvtxCut(10),
  52   fNvtxBins(10),
  53   fNbinsEta(200),
  54   fBackground(0)
  55 {
  56
  57   DefineOutput(1, TList::Class());
  58   DefineOutput(2, TList::Class());
  59   DefineOutput(3, TH1F::Class());
  60   DefineOutput(4, TList::Class());
  61 }
  62 //_____________________________________________________________________
  63 void AliFMDAnalysisTaskGenerateBackground::UserCreateOutputObjects()
  64 {
  65 // Create the output containers
  66 //
  67
  68   std::cout<<"Creating output objects"<<std::endl;
  69   for(Int_t iring = 0; iring<2;iring++) {
  70     Char_t ringChar = (iring == 0 ? 'I' : 'O');
  71     Int_t nSec = (iring == 1 ? 40 : 20);
  72     for(Int_t v=0; v<fNvtxBins;v++) {
  73
  74       TH2F* hPrimary       = new TH2F(Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
  75                                       Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
  76                                       fNbinsEta, -6,6, nSec, 0,2*TMath::Pi());
  77       hPrimary->Sumw2();
  78       fListOfPrimaries.Add(hPrimary);
  79     }
  80   }
  81
  82
  83   for(Int_t det =1; det<=3;det++) {
  84     Int_t nRings = (det==1 ? 1 : 2);
  85     for(Int_t ring = 0;ring<nRings;ring++) {
  86       Int_t nSec = (ring == 1 ? 40 : 20);
  87       Char_t ringChar = (ring == 0 ? 'I' : 'O');
  88       TH1F* doubleHits = new TH1F(Form("DoubleHits_FMD%d%c",det,ringChar),
  89                                   Form("DoubleHits_FMD%d%c",det,ringChar),fNbinsEta, -6,6);
  90       TH1F* allHits = new TH1F(Form("allHits_FMD%d%c",det,ringChar),
  91                                Form("allHits_FMD%d%c",det,ringChar), fNbinsEta, -6,6);
  92
  93       doubleHits->Sumw2();
  94       allHits->Sumw2();
  95       fListOfHits.Add(allHits);
  96       fListOfHits.Add(doubleHits);
  97
  98       for(Int_t v=0; v<fNvtxBins;v++) {
  99         TH2F* hHits = new TH2F(Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
 100                                Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
 101                                fNbinsEta, -6,6, nSec, 0, 2*TMath::Pi());
 102         hHits->Sumw2();
 103         fListOfHits.Add(hHits);
 104
 105       }
 106     }
 107   }
 108
 109   fVertexBins.SetName("VertexBins");
 110   fVertexBins.GetXaxis()->Set(fNvtxBins,-1*fZvtxCut,fZvtxCut);
 111
 112 }
 113 //_____________________________________________________________________
 114 void AliFMDAnalysisTaskGenerateBackground::Init()
 115 {
 116   fLastTrackByStrip.Reset(-1);
 117
 118
 119 }
 120 //_____________________________________________________________________
 121 void AliFMDAnalysisTaskGenerateBackground::UserExec(Option_t */*option*/)
 122 {
 123
 124   fLastTrackByStrip.Reset(-1);
 125   fHitsByStrip.Reset(0);
 126   AliMCEvent* mcevent = MCEvent();
 127   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 128
 129   AliMCParticle* particle = 0;
 130   AliStack* stack = mcevent->Stack();
 131
 132   UShort_t det,sec,strip;
 133   Char_t   ring;
 134
 135   Int_t nTracks = mcevent->GetNumberOfTracks();
 136   AliHeader* header            = mcevent->Header();
 137   AliGenEventHeader* genHeader = header->GenEventHeader();
 138
 139   TArrayF vertex;
 140   genHeader->PrimaryVertex(vertex);
 141
 142   if(TMath::Abs(vertex.At(2)) > fZvtxCut)
 143     return;
 144   for(Int_t i = 0 ;i<nTracks;i++) {
 145     particle = mcevent->GetTrack(i);
 146
 147     if(!particle)
 148       continue;
 149
 150     Double_t delta           = 2*fZvtxCut/fNvtxBins;
 151     Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
 152     Int_t    vertexBin       = (Int_t)vertexBinDouble;
 153
 154     if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
 155
 156
 157       TH2F* hPrimaryInner = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",'I',vertexBin));
 158       TH2F* hPrimaryOuter = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",'O',vertexBin));
 159       hPrimaryInner->Fill(particle->Eta(),particle->Phi());
 160       hPrimaryOuter->Fill(particle->Eta(),particle->Phi());
 161     }
 162
 163     for(Int_t j=0; j<particle->GetNumberOfTrackReferences();j++) {
 164
 165       AliTrackReference* ref = particle->GetTrackReference(j);
 166
 167       if(ref->DetectorId() != AliTrackReference::kFMD)
 168         continue;
 169       AliFMDStripIndex::Unpack(ref->UserId(),det,ring,sec,strip);
 170       Float_t thisStripTrack = fLastTrackByStrip(det,ring,sec,strip);
 171       if(particle->Charge() != 0 && i != thisStripTrack ) {
 172         //Double_t x,y,z;
 173         //AliFMDGeometry* fmdgeo = AliFMDGeometry::Instance();
 174         //fmdgeo->Detector2XYZ(det,ring,sec,strip,x,y,z);
 175         Float_t phi = pars->GetPhiFromSector(det,ring,sec);
 176         Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex.At(2));
 177         //Float_t   phi   = TMath::ATan2(y,x);
 178         //if(phi<0) phi   = phi+2*TMath::Pi();
 179         //      Float_t   r     = TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2));
 180         //Float_t   theta = TMath::ATan2(r,z-vertex.At(2));
 181         //Float_t   eta   = -1*TMath::Log(TMath::Tan(0.5*theta));
 182         TH2F* hHits = (TH2F*)fListOfHits.FindObject(Form("hHits_FMD%d%c_vtx%d", det,ring,vertexBin));
 183         hHits->Fill(eta,phi);
 184         Float_t nstrips = (ring =='O' ? 256 : 512);
 185         fHitsByStrip(det,ring,sec,strip) +=1;
 186         TH1F* allHits = (TH1F*)fListOfHits.FindObject(Form("allHits_FMD%d%c",det,ring));
 187         TH1F* doubleHits = (TH1F*)fListOfHits.FindObject(Form("DoubleHits_FMD%d%c",det,ring));
 188
 189         if(fHitsByStrip(det,ring,sec,strip) == 1)
 190           allHits->Fill(eta);
 191
 192         doubleHits->Fill(eta);
 193         /*if(fHitsByStrip.operator()(det,ring,sec,strip) == 2){
 194           TH1F* doubleHits = (TH1F*)fListOfHits.FindObject("DoubleHits");
 195           doubleHits->Fill(eta,2);
 196           }*/
 197         //if(fHitsByStrip.operator()(det,ring,sec,strip) > 1){
 198         //  doubleHits->Fill(eta);
 199         //      }
 200
 201         fLastTrackByStrip(det,ring,sec,strip) = (Float_t)i;
 202         if(strip >0)
 203           fLastTrackByStrip(det,ring,sec,strip-1) = (Float_t)i;
 204         if(strip < (nstrips - 1))
 205           fLastTrackByStrip(det,ring,sec,strip+1) = (Float_t)i;
 206       }
 207     }
 208
 209   }
 210
 211
 212   PostData(1, &fListOfHits);
 213   PostData(2, &fListOfPrimaries);
 214   PostData(3, &fVertexBins);
 215 }
 216 //_____________________________________________________________________
 217 void AliFMDAnalysisTaskGenerateBackground::Terminate(Option_t */*option*/)
 218 {
 219   /*  TH1F* allHits = (TH1F*)fListOfHits.FindObject("allHits");
 220   TH1F* doubleHits = (TH1F*)fListOfHits.FindObject("DoubleHits");
 221
 222   doubleHits->Divide(allHits);
 223   GenerateCorrection();
 224   PostData(1, &fListOfHits);
 225   PostData(4, &fListOfCorrection);*/
 226
 227 }
 228 //_____________________________________________________________________
 229 void AliFMDAnalysisTaskGenerateBackground::GenerateCorrection() {
 230
 231   fBackground = new AliFMDAnaCalibBackgroundCorrection();
 232
 233   for(Int_t det= 1; det <=3; det++) {
 234     Int_t nRings = (det==1 ? 1 : 2);
 235
 236     for(Int_t iring = 0; iring<nRings; iring++) {
 237       Char_t ring = (iring == 0 ? 'I' : 'O');
 238       TH1F* allHits = (TH1F*)fListOfHits.FindObject(Form("allHits_FMD%d%c",det,ring));
 239       TH1F* doubleHits = (TH1F*)fListOfHits.FindObject(Form("DoubleHits_FMD%d%c",det,ring));
 240       fBackground->SetDoubleHitCorrection(det,ring,doubleHits);
 241       doubleHits->Divide(allHits);
 242       for(Int_t vertexBin=0;vertexBin<fNvtxBins  ;vertexBin++) {
 243         TH2F* hHits          = (TH2F*)fListOfHits.FindObject(Form("hHits_FMD%d%c_vtx%d", det,ring,vertexBin));
 244         TH2F* hPrimary  = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",ring,vertexBin));
 245         TH2F* hCorrection = (TH2F*)hHits->Clone(Form("FMD%d%c_vtxbin_%d_correction",det,ring,vertexBin));
 246         hCorrection->Divide(hPrimary);
 247         hCorrection->SetTitle(hCorrection->GetName());
 248         fListOfCorrection.Add(hCorrection);
 249         fBackground->SetBgCorrection(det,ring,vertexBin,hCorrection);
 250
 251
 252       }
 253
 254     }
 255   }
 256   TAxis refAxis(fNvtxBins,-1*fZvtxCut,fZvtxCut);
 257   fBackground->SetRefAxis(&refAxis);
 258
 259 }
 260 //_____________________________________________________________________
 261 void AliFMDAnalysisTaskGenerateBackground::ReadFromFile(const Char_t* filename, Bool_t storeInOCDB, Int_t /*runNo*/) {
 262
 263   TFile infile(filename);
 264   TH1F* hVertex = (TH1F*)infile.Get("VertexBins");
 265   fZvtxCut = hVertex->GetXaxis()->GetXmax();
 266   fNvtxBins = hVertex->GetXaxis()->GetNbins();
 267   fVertexBins.SetName("VertexBins");
 268   fVertexBins.GetXaxis()->Set(fNvtxBins,-1*fZvtxCut,fZvtxCut);
 269
 270   TList* listOfHits = (TList*)infile.Get("Hits");
 271   TList* listOfPrim = (TList*)infile.Get("Primaries");
 272
 273   for(Int_t det =1; det<=3;det++)
 274       {
 275         Int_t nRings = (det==1 ? 1 : 2);
 276         for(Int_t ring = 0;ring<nRings;ring++)
 277           {
 278             Char_t ringChar = (ring == 0 ? 'I' : 'O');
 279             TH1F* allHits = (TH1F*)listOfHits->FindObject(Form("allHits_FMD%d%c",det,ringChar));
 280             TH1F* doubleHits = (TH1F*)listOfHits->FindObject(Form("DoubleHits_FMD%d%c",det,ringChar));
 281             fListOfHits.Add(allHits);
 282             fListOfHits.Add(doubleHits);
 283             for(Int_t v=0; v<fNvtxBins;v++)
 284               {
 285
 286                 TH2F* hHits          = (TH2F*)listOfHits->FindObject(Form("hHits_FMD%d%c_vtx%d", det,ringChar,v));
 287                 fListOfHits.Add(hHits);
 288               }
 289           }
 290       }
 291   for(Int_t iring = 0; iring<2;iring++) {
 292     Char_t ringChar = (iring == 0 ? 'I' : 'O');
 293     for(Int_t v=0; v<fNvtxBins;v++) {
 294
 295       TH2F* hPrimary       = (TH2F*)listOfPrim->FindObject( Form("hPrimary_FMD_%c_vtx%d",ringChar,v));
 296       fListOfPrimaries.Add(hPrimary);
 297
 298     }
 299   }
 300   GenerateCorrection();
 301
 302   TFile fout("backgroundFromFile.root","recreate");
 303   fListOfHits.Write();
 304   fListOfPrimaries.Write();
 305   fListOfCorrection.Write();
 306   fVertexBins.Write();
 307
 308   fout.Close();
 309
 310   if(storeInOCDB) {
 311     TFile fcalib("$ALICE_ROOT/FMD/Correction/Background/background.root","RECREATE");
 312     fBackground->Write(AliFMDAnaParameters::GetBackgroundID());
 313     fcalib.Close();
 314     /*  AliCDBManager* cdb = AliCDBManager::Instance();
 315     cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
 316     AliCDBId      id(AliFMDAnaParameters::GetBackgroundPath(),runNo,999999999);
 317
 318     AliCDBMetaData* meta = new AliCDBMetaData;
 319     meta->SetResponsible(gSystem->GetUserInfo()->fRealName.Data());
 320     meta->SetAliRootVersion(gROOT->GetVersion());
 321     meta->SetBeamPeriod(1);
 322     meta->SetComment("Background Correction for FMD");
 323     meta->SetProperty("key1", fBackground );
 324     cdb->Put(fBackground, id, meta);
 325     */
 326   }
 327
 328 }
 329 //_____________________________________________________________________
 330 //
 331 // EOF
 332 //