Generation of dictionaries and rootmaps with Root6
[u/mrichter/AliRoot.git] / PWGPP / AliAnalysisTaskIPInfo.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2009, 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 // Class AliAnalysisTaskIPInfo
18 // AliAnalysisTask to extract from the ESD the IP position and sigma
19 // as well as to estimate the primary vertex and tracks DCA resolution.
20 // Uses external class AliIntSpotEstimator
21 //
22 // Author: ruben.shahoyan@cern.ch
23 //*************************************************************************
24
25 #include <TChain.h>
26 #include <TTree.h>
27 #include <TNtuple.h>
28 #include <TBranch.h>
29 #include <TClonesArray.h>
30 #include <TObjArray.h>
31 #include <TH1F.h>
32 #include <TH2F.h>  
33 #include <TNtuple.h>  
34 #include <TCanvas.h>
35
36 #include "AliAnalysisTask.h"
37 #include "AliAnalysisManager.h"
38
39 #include "AliESDtrack.h"
40 #include "AliExternalTrackParam.h"
41 #include "AliESDVertex.h"
42 #include "AliESDEvent.h"
43 #include "AliESDfriend.h"
44 #include "AliVertexerTracks.h"
45 #include "AliESDInputHandler.h"
46 #include "AliAnalysisTaskIPInfo.h"
47 #include "AliIntSpotEstimator.h"
48 #include "AliMultiplicity.h"
49
50
51 ClassImp(AliAnalysisTaskIPInfo)
52
53 const Char_t* AliAnalysisTaskIPInfo::fEstNames[kNEst] = {"ITSTPC","TPC","SPD"};
54
55 //________________________________________________________________________
56 AliAnalysisTaskIPInfo::AliAnalysisTaskIPInfo(const char *name) 
57 : AliAnalysisTask(name, "IP analysis"),
58   fESD(0),fESDfriend(0),fOutput(0),fTracks(50)
59 {
60
61   // Define input and output slots here
62   // Input slot #0 works with a TChain
63   DefineInput(0, TChain::Class());
64   // Output slot #0 writes into a TList container
65   DefineOutput(0, TList::Class());  //My private output
66   //
67   for (int i=0;i<kNEst;i++) fIPEst[i] = 0;
68   //
69 }
70
71 //________________________________________________________________________
72 AliAnalysisTaskIPInfo::~AliAnalysisTaskIPInfo()
73 {
74   // Destructor
75   if (fOutput) {
76     delete fOutput;
77     fOutput = 0;
78   }
79 }
80 //________________________________________________________________________
81 void AliAnalysisTaskIPInfo::SetIPCenIni(Int_t estID, Double_t x,Double_t y,Double_t z)
82 {
83   // set initial estimate of the IP center
84   if (estID<0 || estID>= kNEst) return;
85   fIPCenIni[estID][0] = x;
86   fIPCenIni[estID][1] = y;
87   fIPCenIni[estID][2] = z;
88 }
89
90 //________________________________________________________________________
91 void AliAnalysisTaskIPInfo::SetOptions(Int_t estID, Bool_t recoVtx,
92                                        Double_t outcut,Int_t ntrIP,Int_t nPhiBins,Int_t nestb,
93                                        Double_t estmin,Double_t estmax,
94                                        Int_t ntrBins,Int_t ntMn,Int_t ntMx,
95                                        Int_t nPBins,Double_t pmn,Double_t pmx,Bool_t fillNt)
96 {
97   // set options for estimators
98   if (estID<0 || estID>= kNEst) return;
99   fNTrMinIP[estID] = ntrIP;
100   fRecoVtx[estID]  = recoVtx;
101   fNPhiBins[estID] = nPhiBins;
102   fNEstb[estID]    = nestb;
103   fNTrBins[estID]  = ntrBins;
104   fNPBins[estID]   = nPBins;
105   fNTrMin[estID]   = ntMn;
106   fNTrMax[estID]   = ntMx;
107   fOutCut[estID]   = outcut;
108   fEstMin[estID]   = estmin;
109   fEstMax[estID]   = estmax;
110   fPMin[estID]     = pmn;
111   fPMax[estID]     = pmx;
112   fFillNt[estID]   = fillNt;
113 }
114
115 //________________________________________________________________________
116 void AliAnalysisTaskIPInfo::ConnectInputData(Option_t *) 
117 {
118   // Connect ESD or AOD here
119   // Called once
120   //
121   AliInfo("HERE");
122   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
123   if (!tree) {
124     Printf("ERROR: Could not read chain from input slot 0");
125   } 
126   else {
127     tree->SetBranchAddress("ESDfriend.",&fESDfriend);
128     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());    
129     if (!esdH) Printf("ERROR: Could not get ESDInputHandler");
130     else fESD = esdH->GetEvent();
131   }
132   //
133   return;
134 }
135
136 //________________________________________________________________________
137 void AliAnalysisTaskIPInfo::CreateOutputObjects()
138 {
139   // Create estimators
140   // Several histograms are more conveniently managed in a TList
141   fOutput = new TList;
142   fOutput->SetOwner();
143   //
144   Double_t err[3]={0.0400,0.0400,7.5};
145   //
146   for (int i=0;i<kNEst;i++) if (fNEstb[i]>1) {
147       TString nm   = GetName();
148       TString nmes = fEstNames[i];
149       nm += nmes;
150       fIPEst[i] = new AliIntSpotEstimator(nm.Data(),fOutCut[i],fNTrMinIP[i],fNPhiBins[i],
151                                           fNEstb[i],fEstMin[i],fEstMax[i],
152                                          fNTrBins[i],fNTrMin[i],fNTrMax[i],
153                                          fNPBins[i],fPMin[i],fPMax[i],fFillNt[i]);
154       AliESDVertex *initVertex = new AliESDVertex(fIPCenIni[i],err);
155       fIPEst[i]->GetVertexer()->SetVtxStart(initVertex);
156       delete initVertex;
157       fIPEst[i]->GetVertexer()->SetConstraintOff();
158       fIPEst[i]->GetVertexer()->SetMinClusters(2);
159       fIPEst[i]->SetIPCenIni(fIPCenIni[i]);
160       if (nmes == "TPC") fIPEst[i]->GetVertexer()->SetTPCMode();
161       else               fIPEst[i]->GetVertexer()->SetITSMode();
162       //
163       fOutput->Add(fIPEst[i]);
164       if (fIPEst[i]->GetNtuple()) fOutput->Add(fIPEst[i]->GetNtuple());
165     }
166   //
167   return;
168 }
169
170 //________________________________________________________________________
171 void AliAnalysisTaskIPInfo::Exec(Option_t *) 
172 {
173   static TClonesArray tracks("AliExternalTrackParam",50);
174   //
175   // Main loop
176   // Called for each event
177   if (!fESD) {
178     Printf("ERROR: fESD not available");
179     return;
180   }
181   fESD->SetESDfriend(fESDfriend);
182   //
183   const AliESDVertex *vtx;
184   UShort_t *trackID;
185   Int_t ntracks;
186   //
187   for (int ie=0;ie<kNEst;ie++) {
188     //
189     if (!fIPEst[ie]) continue;
190     if (ie==kTPC) {
191       fIPEst[kTPC]->GetVertexer()->SetFieldkG( fESD->GetMagneticField() );
192       vtx = fRecoVtx[kTPC] ? fIPEst[kTPC]->GetVertexer()->FindPrimaryVertex(fESD) : fESD->GetPrimaryVertexTPC();
193       if (vtx) {
194         ntracks = vtx->GetNIndices();
195         trackID = (UShort_t*)vtx->GetIndices();
196         for (int i=ntracks;i--;) fTracks.Add((TObject*)fESD->GetTrack(trackID[i])->GetTPCInnerParam());
197         fIPEst[kTPC]->ProcessEvent(&fTracks);
198         fTracks.Clear();
199       }
200     }
201     else if (ie==kITSTPC) {
202       fIPEst[kITSTPC]->GetVertexer()->SetFieldkG( fESD->GetMagneticField() );
203       vtx = fRecoVtx[kITSTPC] ? fIPEst[kITSTPC]->GetVertexer()->FindPrimaryVertex(fESD) : fESD->GetPrimaryVertex();
204       if (vtx) {
205         ntracks = vtx->GetNIndices();
206         trackID = (UShort_t*)vtx->GetIndices();
207         for (int i=ntracks;i--;) fTracks.Add((TObject*)fESD->GetTrack(trackID[i]));
208         fIPEst[kITSTPC]->ProcessEvent(&fTracks);
209         fTracks.Clear();
210       }
211     }
212     else if (ie==kSPD) {
213       fIPEst[kSPD]->GetVertexer()->SetFieldkG( fESD->GetMagneticField() );
214       ntracks = CreateSPDTracklets(tracks);
215       for (int i=ntracks;i--;) fTracks.Add((TObject*)tracks[i]);
216       fIPEst[kSPD]->ProcessEvent(&fTracks);
217       fTracks.Clear();
218     }
219   }
220   //
221   PostData(0, fOutput);
222   //
223   return;
224 }      
225
226 //________________________________________________________________________
227 Int_t AliAnalysisTaskIPInfo::CreateSPDTracklets(TClonesArray& tracks)
228 {
229   // create traclets from multiplicity class
230   double cv[21] = {
231     25e-4,
232     0   ,  25e-4,
233     0   ,     0,  40e-2,
234     0   ,     0,      0,   1e-2,
235     0   ,     0,      0,      0,   1e-2,
236     0   ,     0,      0,      0,      0,   1e-2
237   };
238   //
239   double xyzt[3],pxyz[3];
240   int nSPDtracks = 0;
241   tracks.Delete();
242   const AliMultiplicity *mult = fESD->GetMultiplicity();
243   const AliESDVertex *spdVtx  = fESD->GetPrimaryVertexSPD();
244   int nTracklets = 0;
245   if (mult && spdVtx && (nTracklets=mult->GetNumberOfTracklets())>2 ) {
246     const Double_t kRLay1=3.9, kRLay2=7.6;
247     double xv = spdVtx->GetX();
248     double yv = spdVtx->GetY();
249     double zv = spdVtx->GetZ();
250     for (int i=0;i<nTracklets;i++) { // get cluster coordinates from tracklet
251       double phi1 = mult->GetPhi(i);
252       double tht1 = mult->GetTheta(i);
253       double phi2 = phi1 - mult->GetDeltaPhi(i);
254       double tht2 = tht1 - mult->GetDeltaTheta(i);
255       double cs = TMath::Cos(phi1);
256       double sn = TMath::Sin(tht1);
257       double det = xv*sn+yv*sn;
258       det = det*det + kRLay1*kRLay1;
259       double t  = -(xv*cs+yv*sn) + TMath::Sqrt(det);
260       double x1 = cs*t;
261       double y1 = sn*t;
262       double z1 = zv + TMath::Sqrt(x1*x1+y1*y1)/TMath::Tan(tht1);
263       x1 += xv;
264       y1 += yv;
265       //
266       cs = TMath::Cos(phi2);
267       sn = TMath::Sin(tht2);
268       det = xv*sn+yv*sn;
269       det = det*det + kRLay2*kRLay2;
270       t  = -(xv*cs+yv*sn) + TMath::Sqrt(det);
271       double dx = cs*t;
272       double dy = sn*t;
273       double dz = zv + TMath::Sqrt(dx*dx+dy*dy)/TMath::Tan(tht2);
274       dx += xv-x1;
275       dy += yv-y1;
276       dz += -z1;
277       double dr = TMath::Sqrt(dx*dx+dy*dy+dz*dz);
278       pxyz[0] = dx/dr; // direction cosines
279       pxyz[1] = dy/dr;
280       pxyz[2] = dz/dr;
281       t = (xv-x1)*pxyz[0] + (yv-y1)*pxyz[1] + (zv-z1)*pxyz[2];
282       xyzt[0] = x1 + t*pxyz[0];  // PCA to vertex
283       xyzt[1] = y1 + t*pxyz[1];
284       xyzt[2] = z1 + t*pxyz[2];
285       //
286       new(tracks[nSPDtracks++]) AliExternalTrackParam(xyzt,pxyz,cv,0);
287     }
288   }
289   return nSPDtracks;
290 }
291
292 //________________________________________________________________________
293 void AliAnalysisTaskIPInfo::Terminate(Option_t *) 
294 {
295   // Draw result to the screen
296   // Called once at the end of the query
297   fOutput = dynamic_cast<TList*> (GetOutputData(0));
298   if (!fOutput) {     
299     Printf("ERROR: fOutput not available");
300     return;
301   }
302   //
303   return;
304 }