Dump object functionality
[u/mrichter/AliRoot.git] / PWGPP / 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 }