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->GetXv();
 248     double yv = spdVtx->GetYv();
 249     double zv = spdVtx->GetZv();
 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 }