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