PWG1/AliAnalysisTaskVtXY.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //--------------------------------------------------------------------------
  17 //    Implementation of the beam interaction spot location and size estimate
  18 //    via VertexerTracksNoConstraint resolution extrapolation
  19 //    to the infinit multiplicity alias zero resolution
  20 //
  21 // Origin: AliAnalysisTaskVtXY
  22 //         A. Jacholkowski, Catania
  23 //         adam.jacholkowski@cern.ch
  24 //-----------------------------------------------------------------
  25 #include "TChain.h"
  26 #include "TTree.h"
  27 #include "TStyle.h"
  28 #include "TH2F.h"
  29 #include "TProfile.h"
  30 #include "TCanvas.h"
  31
  32 #include "AliAnalysisTask.h"
  33 #include "AliAnalysisManager.h"
  34
  35 #include "AliESDEvent.h"
  36 #include "AliESDInputHandler.h"
  37
  38 #include "AliAnalysisTaskVtXY.h"
  39
  40 #include "AliESDVertex.h"
  41 #include "AliVertexerTracks.h"
  42 ClassImp(AliAnalysisTaskVtXY)
  43
  44 //________________________________________________________________________
  45 AliAnalysisTaskVtXY::AliAnalysisTaskVtXY(const char *name)
  46   : AliAnalysisTask(name, ""),
  47   fESD(0),
  48   fList(0),
  49   fHistVtx(0),
  50   fHistVty(0)
  51 {
  52   // Constructor
  53
  54   // Define input and output slots here
  55   // Input slot #0 works with a TChain
  56   DefineInput(0, TChain::Class());
  57   DefineOutput(0, TList::Class());
  58   // Output slot #0 writes into a TList container
  59 }
  60
  61 //________________________________________________________________________
  62 void AliAnalysisTaskVtXY::ConnectInputData(Option_t *)
  63 {
  64   // Connect ESD or AOD here
  65   // Called once
  66
  67   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
  68   if (!tree) {
  69     Printf("ERROR: Could not read chain from input slot 0");
  70   } else {
  71     // Disable all branches and enable only the needed ones
  72     // The next two lines are different when data produced as AliESDEvent is read
  73     //tree->SetBranchStatus("*", kFALSE);
  74     //tree->SetBranchStatus("Tracks.*", kTRUE);
  75
  76     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
  77
  78     if (!esdH) {
  79       Printf("ERROR: Could not get ESDInputHandler");
  80     } else
  81       fESD = esdH->GetEvent();
  82   }
  83 }
  84
  85 //________________________________________________________________________
  86 void AliAnalysisTaskVtXY::CreateOutputObjects()
  87 {
  88   // Create histograms
  89   // Called once
  90
  91    fHistVtx = new TProfile("fHistVtx","Xvert-RMS vs sigmaX(NC)",100,0.,0.05,-1.,1.,"s");
  92    fHistVty = new TProfile("fHistVty","Yvert-RMS vs sigmaY(NC)",100,0.,0.05,-1.,1.,"s");
  93    fHistVtx->GetXaxis()->SetTitle("sigmaX(NContributors)");
  94    fHistVtx->GetYaxis()->SetTitle("Xv-av&spread");
  95    fHistVty->GetXaxis()->SetTitle("sigmaY(NContributors)");
  96    fHistVty->GetYaxis()->SetTitle("Yv-av&spread");
  97    fHistVtx->SetMarkerStyle(kFullCircle);
  98    fHistVty->SetMarkerStyle(kOpenCircle);
  99    //
 100    fList = new TList();
 101    fList->Add(fHistVtx);
 102    fList->Add(fHistVty);
 103 }
 104
 105 //________________________________________________________________________
 106 void AliAnalysisTaskVtXY::Exec(Option_t *)
 107 {
 108   // Main loop
 109   // Called for each event
 110   // 0 condition / existence of ESD
 111   if (!fESD) {
 112     Printf("ERROR: fESD not available");
 113     return;
 114   }
 115   Printf("There are %d tracks in this event", fESD->GetNumberOfTracks());
 116   Int_t Ngood = 0;
 117   // Track loop kept in case of need
 118   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 119     AliESDtrack* track = fESD->GetTrack(iTracks);
 120     if (!track) {
 121       Printf("ERROR: Could not receive track %d", iTracks);
 122       continue;
 123     }
 124     Ngood++;
 125   } //end track loop
 126   // redo the ESD vertex withour constraint
 127   AliVertexerTracks vertexer(fESD->GetMagneticField());
 128   vertexer.SetITSMode();
 129   AliESDVertex *vertex = vertexer.FindPrimaryVertex(fESD);
 130   // ***************************************************
 131   //  const AliESDVertex *vtxESD = 0;
 132   // if(fESD) {vtxESD = fESD->GetPrimaryVertexTracks();}
 133   Double_t ESDvtx[3] ={999.,999.,9999.};
 134   Double_t XRes,YRes ;
 135   Double_t sig[3];
 136   // Cond 1.vertex has to exist
 137    if (!vertex) {
 138      Printf("ERROR: vertex not available");
 139      return;
 140     }
 141    // Cond 2 - vertexer has not failed
 142    if(!vertex->GetStatus()){
 143      Printf("WARNING: vertexer has failed");
 144        return;
 145    }
 146      TString vtitle = vertex->GetTitle();
 147    // Cond 3 --> vertexer type check
 148      if(!vtitle.Contains("VertexerTracksNoConstraint")){
 149        Printf("WARNING: not VertexerTracksNoConstraint");
 150       return;
 151    }
 152    ULong64_t TrigMask=0;
 153    TrigMask = fESD->GetTriggerMask();
 154    // Cond.4 Trigger mask eventually - not activated
 155    if(TrigMask==0){
 156      Printf("Trigger Mask = 0");
 157      // return;
 158     }
 159     if(vertex){
 160     vertex->GetXYZ(ESDvtx);
 161     XRes = vertex->GetXRes();
 162     YRes = vertex->GetYRes();
 163     vertex->GetSigmaXYZ(sig);
 164     fHistVtx->Fill(XRes,ESDvtx[0]);
 165     fHistVty->Fill(YRes,ESDvtx[1]);
 166     PostData(0,fList);
 167     }
 168   //
 169 }
 170
 171 //________________________________________________________________________
 172 void AliAnalysisTaskVtXY::Terminate(Option_t *)
 173 {
 174   // Draw result to the screen
 175   // Called once at the end of the query
 176   fList = dynamic_cast<TList*>(GetOutputData(0));
 177   if(!fList){
 178     Printf("ERROR: fList not available");
 179     return;
 180   }
 181   fHistVtx = (TProfile *)fList->At(0);
 182   fHistVty = (TProfile *)fList->At(1);
 183   TCanvas *c1 = new TCanvas("AliAnalysisTaskVtXY","Vtx analysis",10,10,800,800);
 184   c1->SetFillColor(10); c1->SetHighLightColor(10);
 185   c1->Divide(1,2);
 186   c1->cd(1);
 187   fHistVtx->DrawCopy("");
 188   c1->cd(2);
 189   fHistVty->DrawCopy("");
 190   // derive histos with fits from the above profiles
 191   Double_t spread = 0.;
 192   Double_t content = 0.;
 193   Double_t entries = 0.;
 194   Double_t error = 0.;
 195   TH1D * hsx = new TH1D("hsx","Xvtx spreads",100,0.,0.050);
 196   //fill
 197   for(Int_t i=0; i<100; i++){
 198     spread = fHistVtx->GetBinError(i);
 199     content = fHistVtx->GetBinContent(i);
 200     entries = fHistVtx->GetBinEntries(i);
 201     error = 0.;
 202     if(entries<10){content = 0.; spread = 0.;}
 203     if(entries>0){error = 0.0005 + spread/TMath::Sqrt(entries);}
 204     hsx->SetBinContent(i,spread);
 205     hsx->SetBinError(i,error);
 206   }
 207   hsx->GetXaxis()->SetTitle("sigma-Xvert");
 208   hsx->GetYaxis()->SetTitle("RMS(Xv) [cm]");
 209   hsx->GetYaxis()->SetTitleOffset(-0.3);
 210   hsx->GetYaxis()->SetRangeUser(0.,0.15);
 211   hsx->SetMarkerColor(3);
 212   hsx->SetMarkerStyle(20);
 213   /*TCanvas *cx = */new TCanvas("cx","",50,50,800,800);
 214   gStyle->SetOptFit(111);
 215   TPad *px = new TPad("px","",0,0,1,1);
 216   px->Draw();
 217   px->cd();
 218   px->SetFillColor(42);
 219   px->SetFrameFillColor(10);
 220   hsx->Fit("pol3","R","",0.001,0.02);
 221   hsx->Draw();
 222   // VtY
 223   TH1D * hsy = new TH1D("hsy","Yvtx spreads",100,0.,0.050);
 224   //fill
 225   for(Int_t i=0; i<100; i++){
 226     spread = fHistVty->GetBinError(i);
 227     content = fHistVty->GetBinContent(i);
 228     entries = fHistVty->GetBinEntries(i);
 229     error = 0.;
 230     if(entries<10){content = 0.; spread = 0.;}
 231     if(entries>0){error = 0.0005 + spread/TMath::Sqrt(entries);}
 232     hsy->SetBinContent(i,spread);
 233     hsy->SetBinError(i,error);
 234   }
 235   hsy->GetXaxis()->SetTitle("sigma-Yvert");
 236   hsy->GetYaxis()->SetTitle("RMS(Yv) [cm]");
 237   hsy->GetYaxis()->SetTitleOffset(-0.3);
 238   hsy->GetYaxis()->SetRangeUser(0.,0.15);
 239   hsy->SetMarkerColor(2);
 240   hsy->SetMarkerStyle(20);
 241   /*TCanvas *cy = */new TCanvas("cy","",100,100,800,800);
 242   TPad *py = new TPad("py","",0,0,1,1);
 243   py->Draw();
 244   py->cd();
 245   py->SetFillColor(42);
 246   py->SetFrameFillColor(10);
 247   hsy->Fit("pol3","R","",0.001,0.02);
 248   hsy->Draw();
 249 }