PWG2/FORWARD/analysis/AliFMDAnalysisTaskCollector.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 <TMath.h>
   9 #include <iostream>
  10
  11 #include "AliFMDAnalysisTaskCollector.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 "AliESDVertex.h"
  20 #include "AliFMDAnaParameters.h"
  21 //#include "AliFMDGeometry.h"
  22
  23
  24 ClassImp(AliFMDAnalysisTaskCollector)
  25
  26 //____________________________________________________________________
  27 Double_t  AliFMDAnalysisTaskCollector::TripleLandau(Double_t *x, Double_t *par) {
  28
  29   Double_t energy        = x[0];
  30   Double_t constant = par[0];
  31   Double_t mpv      = par[1];
  32   Double_t sigma    = par[2];
  33   Double_t alpha    = par[3];
  34   Double_t beta     = par[4];
  35
  36   Double_t f = constant*(TMath::Landau(energy,mpv,sigma,kTRUE)+
  37                          alpha*TMath::Landau(energy,2*mpv+2*sigma*TMath::Log(2),2*sigma,kTRUE)+
  38                          beta*TMath::Landau(energy,3*mpv+3*sigma*TMath::Log(3),3*sigma,kTRUE) );
  39
  40   return f;
  41 }
  42 //____________________________________________________________________
  43
  44 AliFMDAnalysisTaskCollector::AliFMDAnalysisTaskCollector()
  45 : fDebug(0),
  46   fOutputList(0),
  47   fArray(0),
  48   fZvtxDist(0)
  49 {
  50   // Default constructor
  51
  52
  53 }
  54 //____________________________________________________________________
  55 AliFMDAnalysisTaskCollector::AliFMDAnalysisTaskCollector(const char* name):
  56     AliAnalysisTaskSE(name),
  57     fDebug(0),
  58     fOutputList(0),
  59     fArray(0),
  60     fZvtxDist(0)
  61 {
  62   // Default constructor
  63
  64   DefineOutput(1, TList::Class());
  65 }
  66 //____________________________________________________________________
  67 void AliFMDAnalysisTaskCollector::UserCreateOutputObjects()
  68 {
  69   // Create the output container
  70   printf("AnalysisTaskFMD::CreateOutPutData() \n");
  71
  72   fOutputList = new TList();//(TList*)GetOutputData(0);
  73   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  74
  75   fArray     = new TObjArray();
  76   fArray->SetName("FMD");
  77   fArray->SetOwner();
  78   TH1F* hEdist = 0;
  79
  80   for(Int_t nEta = 0; nEta <= pars->GetNetaBins()+1; nEta++) {
  81     TObjArray* etaArray = new TObjArray();
  82     fArray->AddAtAndExpand(etaArray,nEta);
  83     for(Int_t det =1; det<=3;det++)
  84       {
  85         TObjArray* detArray = new TObjArray();
  86         detArray->SetName(Form("FMD%d",det));
  87         etaArray->AddAtAndExpand(detArray,det);
  88         Int_t nRings = (det==1 ? 1 : 2);
  89         for(Int_t ring = 0;ring<nRings;ring++)
  90           {
  91             Char_t ringChar = (ring == 0 ? 'I' : 'O');
  92             hEdist = new TH1F(Form("FMD%d%c_etabin%d",det,ringChar,nEta),Form("FMD%d%c_etabin%d",det,ringChar,nEta),200,0,6);
  93             hEdist->SetXTitle("#Delta E / E_{MIP}");
  94             fOutputList->Add(hEdist);
  95             detArray->AddAtAndExpand(hEdist,ring);
  96           }
  97       }
  98
  99   }
 100
 101   fZvtxDist  = new TH1F("ZvtxDist","Vertex distribution",100,-30,30);
 102   fZvtxDist->SetXTitle("z vertex");
 103   fOutputList->Add(fZvtxDist);
 104 }
 105
 106 //____________________________________________________________________
 107 void AliFMDAnalysisTaskCollector::UserExec(Option_t */*option*/)
 108 {
 109
 110   AliESDEvent* esd = (AliESDEvent*)InputEvent();
 111   AliESD* old = esd->GetAliESDOld();
 112   if (old) {
 113     esd->CopyFromOldESD();
 114   }
 115   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 116
 117   Bool_t trigger = pars->IsEventTriggered(esd);
 118   if(!trigger)
 119     return;
 120   Double_t vertex[3];
 121
 122   pars->GetVertex(esd,vertex);
 123   if(vertex[0] == 0 && vertex[1] == 0 && vertex[2] == 0)
 124     return;
 125
 126   fZvtxDist->Fill(vertex[2]);
 127
 128   if(TMath::Abs(vertex[2]) > pars->GetVtxCutZ())
 129     return;
 130
 131   AliESDFMD* fmd = esd->GetFMDData();
 132   if (!fmd) return;
 133
 134   for(UShort_t det=1;det<=3;det++) {
 135
 136     Int_t nRings = (det==1 ? 1 : 2);
 137     for (UShort_t ir = 0; ir < nRings; ir++) {
 138
 139       Char_t   ring = (ir == 0 ? 'I' : 'O');
 140       UShort_t nsec = (ir == 0 ? 20  : 40);
 141       UShort_t nstr = (ir == 0 ? 512 : 256);
 142       TH2F* hBg = pars->GetBackgroundCorrection(det,ring,0);
 143
 144       for(UShort_t sec =0; sec < nsec;  sec++)  {
 145         for(UShort_t strip = 0; strip < nstr; strip++) {
 146
 147
 148           Float_t mult = fmd->Multiplicity(det,ring,sec,strip);
 149           if(mult == AliESDFMD::kInvalidMult || mult == 0) continue;
 150
 151           Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex[2]);
 152
 153           Int_t nEta = hBg->GetXaxis()->FindBin(eta);
 154
 155           TObjArray* etaArray = (TObjArray*)fArray->At(nEta);
 156           TObjArray* detArray = (TObjArray*)etaArray->At(det);
 157           TH1F* Edist = (TH1F*)detArray->At(ir);
 158
 159           Edist->Fill(mult);
 160
 161         }
 162       }
 163     }
 164   }
 165
 166   PostData(1, fOutputList);
 167
 168 }
 169
 170 //____________________________________________________________________
 171 void AliFMDAnalysisTaskCollector::ReadFromFile(const Char_t* filename, Bool_t store, Int_t speciesOption)  {
 172
 173   //speciesOption:
 174   //0: p+p Landau fit
 175   //1: Pb+Pb triple landau convolution fit
 176
 177   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 178   pars->Init(kTRUE,AliFMDAnaParameters::kBackgroundCorrection);
 179
 180   TFile fin(filename,"UPDATE");
 181
 182   TList* list = (TList*)fin.Get("energyDist");
 183
 184   AliFMDAnaCalibEnergyDistribution* EnergyDist = new AliFMDAnaCalibEnergyDistribution();
 185
 186   EnergyDist->SetNetaBins(pars->GetNetaBins());
 187   EnergyDist->SetEtaLimits(pars->GetEtaMin(),pars->GetEtaMax());
 188
 189   for(Int_t nEta = 1; nEta <= pars->GetNetaBins(); nEta++) {
 190
 191     for(Int_t det = 1; det<=3; det++) {
 192       Int_t nRings  =  (det == 1 ? 1 : 2);
 193       for(Int_t ring = 0;ring<nRings; ring++) {
 194         Char_t ringChar = (ring == 0 ? 'I' : 'O');
 195
 196         TH1F* hEdist = (TH1F*)list->FindObject(Form("FMD%d%c_etabin%d",det,ringChar,nEta));
 197         TF1* fitFunc = 0 ;
 198
 199         if(hEdist->GetEntries() != 0) {
 200
 201           hEdist->GetXaxis()->SetRangeUser(0.2,hEdist->GetXaxis()->GetXmax());
 202
 203           if(speciesOption == 0)
 204             fitFunc =  new TF1("FMDfitFunc","landau",hEdist->GetBinCenter(hEdist->GetMaximumBin())-0.2,3);
 205           if(speciesOption == 1) {
 206             fitFunc = new TF1("FMDfitFunc",TripleLandau,hEdist->GetBinCenter(hEdist->GetMaximumBin())-0.2,5,5);
 207             fitFunc->SetParNames("constant","MPV","sigma","2-Mip weight","3-Mip weight");
 208             fitFunc->SetParameters(10,0.8,0.1,0.05,0.01);
 209             fitFunc->SetParLimits(1,0.6,1.2);
 210             fitFunc->SetParLimits(3,0,1);
 211             fitFunc->SetParLimits(4,0,1);
 212
 213           }
 214
 215
 216           hEdist->Fit(fitFunc,"","",hEdist->GetBinCenter(hEdist->GetMaximumBin())-0.2,3);
 217           fitFunc->Write(Form("FMD%d%c_etabin%d_fitfunc",det,ringChar,nEta),TObject::kWriteDelete);
 218
 219         }
 220
 221         TH2F* hBg = pars->GetBackgroundCorrection(det,ringChar,0);
 222         EnergyDist->SetEnergyDistribution(det,ringChar,hBg->GetXaxis()->GetBinCenter(nEta),hEdist);
 223       }
 224
 225     }
 226
 227   }
 228
 229   fin.Write();
 230   fin.Close();
 231
 232   if(store) {
 233     TFile fcalib(pars->GetPath(pars->GetEdistID() ),"RECREATE");
 234     EnergyDist->Write(AliFMDAnaParameters::GetEdistID());
 235     fcalib.Close();
 236   }
 237
 238
 239
 240 }
 241
 242 //____________________________________________________________________
 243 //
 244 // EOF
 245 //