FMD/analysis/AliFMDAnalysisTaskDndeta.cxx

   1
   2 #include <TROOT.h>
   3 #include <TSystem.h>
   4 #include <TInterpreter.h>
   5 #include <TChain.h>
   6 #include <TFile.h>
   7 #include <TList.h>
   8 #include <iostream>
   9 #include "TH1F.h"
  10 #include "TH2F.h"
  11 #include "AliFMDAnalysisTaskDndeta.h"
  12 #include "AliAnalysisManager.h"
  13 #include "AliESDFMD.h"
  14 #include "AliESDEvent.h"
  15 #include "AliAODEvent.h"
  16 #include "AliAODHandler.h"
  17 #include "AliMCEventHandler.h"
  18 #include "AliStack.h"
  19 #include "AliLog.h"
  20 #include "AliESDVertex.h"
  21 #include "TMath.h"
  22 #include "AliFMDAnaParameters.h"
  23 //#include "AliFMDGeometry.h"
  24 #include "AliGenEventHeader.h"
  25 #include "AliHeader.h"
  26 //#include "TDatabasePDG.h"
  27 //#include "TParticlePDG.h"
  28 #include "AliFMDStripIndex.h"
  29 ClassImp(AliFMDAnalysisTaskDndeta)
  30
  31
  32 AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta()
  33 : fDebug(0),
  34   fOutputList(0),
  35   fInputList(0),
  36 // fArray(0),
  37   fInputArray(0),
  38   fVertexString(0x0),
  39   fNevents(),
  40   fNMCevents(),
  41   fStandalone(kTRUE),
  42 //fMCevent(0),
  43   fLastTrackByStrip(0)
  44 //  fPrimary(kTRUE),
  45 //  fRecordHits(kFALSE)
  46 {
  47   // Default constructor
  48   DefineInput (0, TList::Class());
  49   DefineOutput(0, TList::Class());
  50 }
  51 //_____________________________________________________________________
  52 AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta(const char* name, Bool_t SE):
  53     AliAnalysisTask(name, "Density"),
  54     fDebug(0),
  55     fOutputList(0),
  56     fInputList(0),
  57     //  fArray(),
  58     fInputArray(0),
  59     fVertexString(0x0),
  60     fNevents(),
  61     fNMCevents(),
  62     fStandalone(kTRUE),
  63     //  fMCevent(0),
  64     fLastTrackByStrip(0)
  65     //   fPrimary(kTRUE),
  66     //  fRecordHits(kFALSE)
  67 {
  68   fStandalone = SE;
  69   if(fStandalone) {
  70     DefineInput (0, TList::Class());
  71     DefineInput(1, TObjString::Class());
  72     DefineOutput(0, TList::Class());
  73
  74   }
  75 }
  76 //_____________________________________________________________________
  77 void AliFMDAnalysisTaskDndeta::CreateOutputObjects()
  78 {
  79   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  80
  81   // fArray.SetName("FMD");
  82   //  fArray.SetOwner();
  83
  84   if(!fOutputList)
  85     fOutputList = new TList();
  86   fOutputList->SetName("BackgroundCorrected");
  87
  88
  89   TH2F* hMult = 0;
  90   TH2F* hMultTrVtx = 0;
  91   TH1F* hHits = 0;
  92   TH1F* hPrimVertexBin = 0;
  93
  94
  95   TH2F* hBgTmp   = pars->GetBackgroundCorrection(1, 'I', 0);
  96   TH1F* hPrimary = new TH1F("hMultvsEta","hMultvsEta",
  97                             hBgTmp->GetNbinsX(),
  98                             hBgTmp->GetXaxis()->GetXmin(),
  99                             hBgTmp->GetXaxis()->GetXmax());
 100   hPrimary->Sumw2();
 101   fOutputList->Add(hPrimary);
 102   Int_t nVtxbins = pars->GetNvtxBins();
 103
 104
 105   for(Int_t det =1; det<=3;det++)
 106     {
 107       // TObjArray* detArray = new TObjArray();
 108       // detArray->SetName(Form("FMD%d",det));
 109       // fArray.AddAtAndExpand(detArray,det);
 110       Int_t nRings = (det==1 ? 1 : 2);
 111       for(Int_t ring = 0;ring<nRings;ring++)
 112         {
 113           Char_t ringChar = (ring == 0 ? 'I' : 'O');
 114           Int_t  nSec     = (ring == 0 ? 20 : 40);
 115
 116           //  TObjArray* vtxArray = new TObjArray();
 117           // vtxArray->SetName(Form("FMD%d%c",det,ringChar));
 118           // detArray->AddAtAndExpand(vtxArray,ring);
 119           for(Int_t i = 0; i< nVtxbins; i++) {
 120             TH2F* hBg = pars->GetBackgroundCorrection(det, ringChar, i);
 121             hMult  = new TH2F(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),
 122                               hBg->GetNbinsX(),
 123                               hBg->GetXaxis()->GetXmin(),
 124                               hBg->GetXaxis()->GetXmax(),
 125                               nSec, 0, 2*TMath::Pi());
 126
 127             hHits  = new TH1F(Form("hHits_FMD%d%c_vtxbin%d",det,ringChar,i),Form("hHits_FMD%d%c_vtxbin%d",det,ringChar,i),
 128                               hBg->GetNbinsX(),
 129                               hBg->GetXaxis()->GetXmin(),
 130                               hBg->GetXaxis()->GetXmax());
 131
 132
 133
 134             hMult->Sumw2();
 135             hHits->Sumw2();
 136             fOutputList->Add(hMult);
 137             // fOutputList->Add(hMultTrVtx);
 138             fOutputList->Add(hHits);
 139             //vtxArray->AddAtAndExpand(hMult,i);
 140
 141           }
 142         }
 143     }
 144
 145   for(Int_t i = 0; i< nVtxbins; i++) {
 146
 147     hPrimVertexBin = new TH1F(Form("primmult_vtxbin%d",i),
 148                               Form("primmult_vtxbin%d",i),
 149                               hBgTmp->GetNbinsX(),
 150                               hBgTmp->GetXaxis()->GetXmin(),
 151                               hBgTmp->GetXaxis()->GetXmax());
 152     hPrimVertexBin->Sumw2();
 153     fOutputList->Add(hPrimVertexBin);
 154
 155   }
 156
 157   fNevents.SetBins(nVtxbins,0,nVtxbins);
 158   fNevents.SetName("nEvents");
 159   fNMCevents.SetBins(nVtxbins,0,nVtxbins);
 160   fNMCevents.SetName("nMCEvents");
 161   fOutputList->Add(&fNevents);
 162   fOutputList->Add(&fNMCevents);
 163
 164 }
 165 //_____________________________________________________________________
 166 void AliFMDAnalysisTaskDndeta::ConnectInputData(Option_t */*option*/)
 167 {
 168   if(fStandalone) {
 169     fInputList   = (TList*)GetInputData(0);
 170     fVertexString = (TObjString*)GetInputData(1);
 171   }
 172 }
 173 //_____________________________________________________________________
 174 void AliFMDAnalysisTaskDndeta::Exec(Option_t */*option*/)
 175 {
 176   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 177   Int_t vtxbin   = fVertexString->GetString().Atoi();
 178   fNevents.Fill(vtxbin);
 179   for(UShort_t det=1;det<=3;det++) {
 180     //TObjArray* detInputArray = (TObjArray*)fInputArray->At(det);
 181     // TObjArray* detArray = (TObjArray*)fArray.At(det);
 182     Int_t nRings = (det==1 ? 1 : 2);
 183     for (UShort_t ir = 0; ir < nRings; ir++) {
 184       Char_t ringChar = (ir == 0 ? 'I' : 'O');
 185       //TObjArray* vtxInputArray = (TObjArray*)detInputArray->At(ir);
 186       //TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
 187       // TH2F* hMultTotal = (TH2F*)vtxArray->At(vtxbin);
 188
 189       TH2F* hMultTotal = (TH2F*)fOutputList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
 190
 191
 192       TH2F* hMultInput = (TH2F*)fInputList->FindObject(Form("mult_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
 193
 194       hMultTotal->Add(hMultInput);
 195       //hMultTrVtx->Add(hMultInput);
 196
 197     }
 198   }
 199
 200   if(pars->GetProcessPrimary())
 201     ProcessPrimary();
 202
 203   if(fStandalone) {
 204     PostData(0, fOutputList);
 205   }
 206
 207 }
 208 //_____________________________________________________________________
 209 void AliFMDAnalysisTaskDndeta::Terminate(Option_t */*option*/) {
 210
 211
 212   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 213
 214   Int_t nVtxbins = pars->GetNvtxBins();
 215
 216   for(UShort_t det=1;det<=3;det++) {
 217     // TObjArray* detArray = (TObjArray*)fArray.At(det);
 218     Int_t nRings = (det==1 ? 1 : 2);
 219     for (UShort_t ir = 0; ir < nRings; ir++) {
 220       // TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
 221       Char_t ringChar = (ir == 0 ? 'I' : 'O');
 222       for(Int_t i =0; i<nVtxbins; i++) {
 223
 224         //TH2F* hMultTotal = (TH2F*)vtxArray->At(i);
 225         //hMultTotal->Scale(pars->GetTriggerEfficiency(i));
 226         TH2F* hMultTotal = (TH2F*)fOutputList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i));
 227         TH2F* hMultTrVtx = (TH2F*)hMultTotal->Clone(Form("dNdeta_FMD_TrVtx_%d%c_vtxbin%d",det,ringChar,i));
 228
 229         hMultTotal->Scale(pars->GetEventSelectionEfficiency(i));
 230         //std::cout<<pars->GetTriggerEfficiency(i)<<std::endl;
 231         std::cout<<pars->GetEventSelectionEfficiency(i)<<std::endl;
 232         TH1D* hMultProj   = hMultTotal->ProjectionX(Form("dNdeta_FMD%d%c_vtxbin%d_proj",det,ringChar,i),1,hMultTotal->GetNbinsY());
 233         TH1D* hMultProjTrVtx   = hMultTrVtx->ProjectionX(Form("dNdeta_FMD_TrVtx_%d%c_vtxbin%d_proj",det,ringChar,i),1,hMultTotal->GetNbinsY());
 234         fOutputList->Add(hMultTrVtx);
 235         fOutputList->Add(hMultProj);
 236         fOutputList->Add(hMultProjTrVtx);
 237       }
 238     }
 239   }
 240
 241 }
 242 //_____________________________________________________________________
 243 void AliFMDAnalysisTaskDndeta::ProcessPrimary() {
 244
 245   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 246   AliMCEvent* mcEvent = eventHandler->MCEvent();
 247   if(!mcEvent)
 248     return;
 249
 250   fLastTrackByStrip.Reset(-1);
 251
 252   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 253
 254   AliMCParticle* particle = 0;
 255   AliStack* stack = mcEvent->Stack();
 256
 257   TH1F* hPrimary = (TH1F*)fOutputList->FindObject("hMultvsEta");
 258   AliHeader* header            = mcEvent->Header();
 259   AliGenEventHeader* genHeader = header->GenEventHeader();
 260
 261   TArrayF vertex;
 262   genHeader->PrimaryVertex(vertex);
 263   if(TMath::Abs(vertex.At(2)) > pars->GetVtxCutZ())
 264     return;
 265   Double_t delta           = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
 266   Double_t vertexBinDouble = (vertex.At(2) + pars->GetVtxCutZ()) / delta;
 267   Int_t    vertexBin       = (Int_t)vertexBinDouble;
 268
 269   Bool_t firstTrack = kTRUE;
 270
 271   // we loop over the primaries only unless we need the hits (diagnostics running slowly)
 272   Int_t nTracks = stack->GetNprimary();
 273   if(pars->GetProcessHits())
 274     nTracks = stack->GetNtrack();
 275
 276   for(Int_t i = 0 ;i<nTracks;i++) {
 277     particle = mcEvent->GetTrack(i);
 278     if(!particle)
 279       continue;
 280
 281     if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
 282       hPrimary->Fill(particle->Eta());
 283
 284
 285       TH1F* hPrimVtxBin = (TH1F*)fOutputList->FindObject(Form("primmult_vtxbin%d",vertexBin));
 286       hPrimVtxBin->Fill(particle->Eta());
 287       if(firstTrack) {
 288         fNMCevents.Fill(vertexBin);
 289         firstTrack = kFALSE;
 290       }
 291
 292     }
 293     if(pars->GetProcessPrimary()) {
 294       for(Int_t j=0; j<particle->GetNumberOfTrackReferences();j++) {
 295
 296       AliTrackReference* ref = particle->GetTrackReference(j);
 297       UShort_t det,sec,strip;
 298       Char_t   ring;
 299       if(ref->DetectorId() != AliTrackReference::kFMD)
 300         continue;
 301       AliFMDStripIndex::Unpack(ref->UserId(),det,ring,sec,strip);
 302       Float_t thisStripTrack = fLastTrackByStrip(det,ring,sec,strip);
 303       if(particle->Charge() != 0 && i != thisStripTrack ) {
 304         //Double_t x,y,z;
 305
 306         Float_t   eta   = pars->GetEtaFromStrip(det,ring,sec,strip,vertex.At(2));//-1*TMath::Log(TMath::Tan(0.5*theta));
 307         TH1F* hHits = (TH1F*)fOutputList->FindObject(Form("hHits_FMD%d%c_vtxbin%d",det,ring,vertexBin));
 308         hHits->Fill(eta);
 309         Float_t nstrips = (ring =='O' ? 256 : 512);
 310
 311         fLastTrackByStrip(det,ring,sec,strip) = (Float_t)i;
 312
 313         if(strip >0)
 314           fLastTrackByStrip(det,ring,sec,strip-1) = (Float_t)i;
 315         if(strip < (nstrips - 1))
 316           fLastTrackByStrip(det,ring,sec,strip+1) = (Float_t)i;
 317
 318
 319       }
 320       }
 321     }
 322
 323   }
 324
 325 }
 326 //_____________________________________________________________________
 327 //
 328 //
 329 // EOF