PWGLF/FORWARD/trains/MyAnalysis.C

   1 // MyAnalysis.C
   2 #ifndef __CINT__
   3 # include <AliAnalysisManager.h>
   4 # include <AliESDEvent.h>
   5 # include <AliMultiplicity.h>
   6 # include <AliVEventHandler.h>
   7 # include <AliESDVertex.h>
   8 # include <AliProdInfo.h>
   9 # include <TH1D.h>
  10 # include <TH2D.h>
  11 #else
  12 class TH1D;
  13 class TH2D;
  14 class AliProdInfo;
  15 #endif
  16 #include <AliAnalysisTaskSE.h>
  17 class MyAnalysis : public AliAnalysisTaskSE
  18 {
  19 public:
  20   MyAnalysis()
  21     : AliAnalysisTaskSE(), fList(0), fMult(0), fVz(0), fProd(0)
  22   {}
  23   MyAnalysis(const char* name)
  24     : AliAnalysisTaskSE(name), fList(0), fMult(0), fVz(0), fProd(0)
  25   {
  26     DefineOutput(1, TList::Class());
  27     DefineOutput(2, TList::Class()); // For output from Terminate
  28     fBranchNames = "AliMultiplicity.,SPDVertex.,PrimaryVertex.";
  29   }
  30   MyAnalysis(const MyAnalysis& o)
  31   : AliAnalysisTaskSE(o), fList(o.fList), fMult(o.fMult), fVz(o.fVz)
  32   {}
  33   virtual ~MyAnalysis() {}
  34   MyAnalysis& operator=(const MyAnalysis&) { return *this; }
  35   virtual void UserCreateOutputObjects()
  36   {
  37     fList = new TList();
  38     fList->SetName("Sums");
  39     fList->SetOwner();
  40
  41     fMult = new TH2D("mult", "SPD tracklets", 80, -2, 2, 10, -10, 10);
  42     fMult->SetXTitle("#eta");
  43     fMult->SetYTitle("v_{z} [cm]");
  44     fMult->Sumw2();
  45     fMult->SetDirectory(0); // Disassociate from file
  46
  47     fVz = new TH1D("vz", "Interaction point", 10, -10, 10);
  48     fVz->SetXTitle("v_{z} [cm]");
  49     fVz->Sumw2();
  50     fVz->SetDirectory(0); // Disassociate from file
  51
  52     fList->Add(fMult);
  53     fList->Add(fVz);
  54
  55     PostData(1, fList);
  56   }
  57   virtual void UserExec(Option_t* )
  58   {
  59     if (!fProd) {
  60       AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
  61       AliVEventHandler *inputHandler=mgr->GetInputEventHandler();
  62       if (inputHandler) {
  63         Info("", "Got input handler");
  64         TList *uiList = inputHandler->GetUserInfo();
  65         if (uiList) {
  66           Info("", "Got user list:");
  67           uiList->ls();
  68
  69           fProd = new AliProdInfo(uiList);
  70           Info("", "Lising production information");
  71           fProd->List();
  72         }
  73       }
  74     }
  75
  76     AliESDEvent* event = dynamic_cast<AliESDEvent*>(InputEvent());
  77     if (!event) return;
  78     if (event->IsPileupFromSPD(3,0.8))   return;
  79
  80     const AliESDVertex* vtx = event->GetPrimaryVertexSPD();
  81     if (!vtx || !vtx->GetStatus()) return;
  82     if (vtx->IsFromVertexerZ() &&
  83         (vtx->GetDispersion() > 0.2 ||  vtx->GetZRes() > 1.25 * 0.2))
  84       return;
  85
  86     const AliMultiplicity* mult = event->GetMultiplicity();
  87     if (!mult) return;
  88
  89     Double_t vz = vtx->GetZ();
  90     fVz->Fill(vz);
  91
  92     Int_t nTracklets = mult->GetNumberOfTracklets();
  93     for (Int_t i = 0; i < nTracklets; i++)
  94       fMult->Fill(mult->GetEta(i), vz);
  95
  96     PostData(1, fList);
  97   }
  98   void  Terminate(Option_t *)
  99   {
 100     TList* l = dynamic_cast<TList*>(GetOutputData(1));
 101     if (!l) {
 102       Warning("Terminate", "No out data # 1 found");
 103       return;
 104     }
 105
 106     TH2D* mult = static_cast<TH2D*>(l->FindObject("mult"));
 107     TH1D* vz   = static_cast<TH1D*>(l->FindObject("vz"));
 108     if (!mult || !vz) {
 109       Warning("Terminate", "Either 'mult' (%p) or 'vz' (%p) or both not found",
 110               mult, vz);
 111       return;
 112     }
 113
 114     TList* output = new TList;   // Needed for new output from Terminate
 115     output->SetName("Results");  // 1st output re-opened read-only
 116     output->SetOwner();
 117
 118     TH2D* out = static_cast<TH2D*>(mult->Clone("dndeta"));
 119     out->SetTitle("dN_{ch}/d#eta from SPD tracklets per vertex bin");
 120     out->SetZTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
 121     out->SetDirectory(0); // Disassociate from file
 122     Int_t    nVz  = mult->GetNbinsY();
 123     Int_t    nEta = mult->GetNbinsX();
 124     for (Int_t iVz = 1; iVz <= nVz; iVz++) {
 125       Double_t nEv = vz->GetBinContent(iVz);
 126       Double_t e1  = vz->GetBinError(iVz);
 127       Double_t sca = (nEv == 0 ? 0 : 1. / nEv);
 128       for (Int_t iEta = 1; iEta <= nEta; iEta++) {
 129         Double_t c  = mult->GetBinContent(iEta,iVz);
 130         Double_t e  = mult->GetBinError(iEta,iVz);
 131         Double_t ee = TMath::Sqrt(c*c * e1*e1 + nEv*nEv * e*e) * sca*sca;
 132         out->SetBinContent(iEta, iVz, sca * c);
 133         out->SetBinError(iEta, iVz, ee);
 134       }
 135     }
 136     Double_t etaMin = mult->GetXaxis()->GetXmin();
 137     Double_t etaMax = mult->GetXaxis()->GetXmax();
 138     out->Scale(Double_t(nEta) / (etaMax-etaMin));
 139
 140     output->Add(out);
 141     PostData(2, output);
 142   }
 143 protected:
 144   TList*  fList;
 145   TH2D*   fMult;
 146   TH1D*   fVz;
 147   AliProdInfo* fProd;
 148   ClassDef(MyAnalysis, 1);
 149 };
 150 //
 151 // EOF
 152 //
 153