PWGPP/ITS/AliMeanVertexCalibTask.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 AliMeanVertexCalibTask
  18 // AliAnalysisTask to extract from ESD the information on primary vertex
  19 // reconstruction in order to compute the MeanVertex object
  20 //
  21 // Author:  D.Caffarri, davide.caffarri@pd.infn.it
  22 //          A.Dainese, andrea.dainese@pd.infn.it
  23 //*************************************************************************
  24
  25
  26 #include <TH1F.h>
  27 #include <TH2F.h>
  28
  29 #include "AliMultiplicity.h"
  30 #include "AliESDVertex.h"
  31 #include "AliESDEvent.h"
  32 #include "AliVertexerTracks.h"
  33
  34 #include "AliMeanVertexCalibTask.h"
  35
  36
  37 ClassImp(AliMeanVertexCalibTask)
  38
  39 //_____________________________________________________________________
  40 AliMeanVertexCalibTask::AliMeanVertexCalibTask(const char *name):
  41   AliAnalysisTaskSE(name),
  42   fESD(0),
  43   fOutput(0),
  44   fOnlyITSTPCTracks(kFALSE),
  45   fOnlyITSSATracks(kTRUE)
  46 {
  47
  48   // Constructor
  49
  50   // Define input and output slots here
  51   // Output slot #0 writes into a TList container
  52   DefineOutput(1, TList::Class());  //My private output
  53 }
  54
  55 //________________________________________________________________________
  56 AliMeanVertexCalibTask::~AliMeanVertexCalibTask()
  57 {
  58   // Destructor
  59
  60   if (fOutput) {
  61     delete fOutput;
  62     fOutput = 0;
  63   }
  64
  65 }
  66 //________________________________________________________________________
  67 void AliMeanVertexCalibTask::UserCreateOutputObjects()
  68 {
  69
  70   // Several histograms are more conveniently managed in a TList
  71   fOutput = new TList;
  72   fOutput->SetOwner();
  73
  74
  75   TH1F* hSPDVertexX = new TH1F("hSPDVertexX","SPDVertex x; x vertex [cm]; events",200,-1,1);
  76   fOutput->Add(hSPDVertexX);
  77   TH1F* hSPDVertexY = new TH1F("hSPDVertexY","SPDVertex y; y vertex [cm]; events",200,-1,1);
  78   fOutput->Add(hSPDVertexY);
  79   TH1F* hSPDVertexZ = new TH1F("hSPDVertexZ","SPDVertex z; z vertex [cm]; events",200,-20,20);
  80   fOutput->Add(hSPDVertexZ);
  81   TH1F* hTRKVertexX = new TH1F("hTRKVertexX","TRKVertex x; x vertex [cm]; events",200,-1,1);
  82   fOutput->Add(hTRKVertexX);
  83   TH1F* hTRKVertexY = new TH1F("hTRKVertexY","TRKVertex y; y vertex [cm]; events",200,-1,1);
  84   fOutput->Add(hTRKVertexY);
  85   TH1F* hTRKVertexZ = new TH1F("hTRKVertexZ","TRKVertex z; z vertex [cm]; events",200,-20,20);
  86   fOutput->Add(hTRKVertexZ);
  87
  88   TH2F *hTRKVertexXZ = new TH2F("hTRKVertexXZ", "TRKVertex XZ corr", 200, -1, 1, 200, -20, 20);
  89   fOutput->Add(hTRKVertexXZ);
  90
  91   TH2F *hTRKVertexYZ = new TH2F("hTRKVertexYZ", "TRKVertex YZ corr", 200, -1, 1, 200, -20, 20);
  92   fOutput->Add(hTRKVertexYZ);
  93
  94   TH1F* hTRKVertexXdefMult = new TH1F("hTRKVertexXdefMult","TRKVertex x Mult; x vertex [cm] 30<Mult<45; events",500,-1,1);
  95   fOutput->Add(hTRKVertexXdefMult);
  96   TH1F* hTRKVertexYdefMult = new TH1F("hTRKVertexYdefMult","TRKVertex y Mult; y vertex [cm] 30<Mult<45; events",500,-1,1);
  97   fOutput->Add(hTRKVertexYdefMult);
  98
  99         TH1F* hTRKVertexXHighMult = new TH1F("hTRKVertexXHighMult","TRKVertex x High Mult; x vertex [cm] Mult>1500; events",500,-0.5,0.5);
 100         fOutput->Add(hTRKVertexXHighMult);
 101         TH1F* hTRKVertexYHighMult = new TH1F("hTRKVertexYHighMult","TRKVertex y High Mult; y vertex [cm] Mult>1500; events",500,-0.5,0.5);
 102         fOutput->Add(hTRKVertexYHighMult);
 103
 104
 105   TH1F* hITSSAVertexX = new TH1F("hITSSAVertexX","ITSSAVertex x; x vertex [cm]; events",200,-1,1);
 106   fOutput->Add(hITSSAVertexX);
 107   TH1F* hITSSAVertexY = new TH1F("hITSSAVertexY","ITSSAVertex y; y vertex [cm]; events",200,-1,1);
 108   fOutput->Add(hITSSAVertexY);
 109   TH1F* hITSSAVertexZ = new TH1F("hITSSAVertexZ","ITSSAVertex z; z vertex [cm]; events",200,-20,20);
 110   fOutput->Add(hITSSAVertexZ);
 111
 112   TH2F *hITSSAVertexXZ = new TH2F("hITSSAVertexXZ", "ITSSAVertex XZ corr", 200, -1, 1, 200, -20, 20);
 113   fOutput->Add(hITSSAVertexXZ);
 114
 115   TH2F *hITSSAVertexYZ = new TH2F("hITSSAVertexYZ", "ITSSAVertex YZ corr", 200, -1, 1, 200, -20, 20);
 116   fOutput->Add(hITSSAVertexYZ);
 117
 118   TH1F* hITSSAVertexXdefMult = new TH1F("hITSSAVertexXdefMult","ITSSAVertex x Mult; x vertex [cm] 30<Mult<45; events",500,-1,1);
 119   fOutput->Add(hITSSAVertexXdefMult);
 120   TH1F* hITSSAVertexYdefMult = new TH1F("hITSSAVertexYdefMult","ITSSAVertex y Mult; y vertex [cm] 30<Mult<45; events",500,-1,1);
 121   fOutput->Add(hITSSAVertexYdefMult);
 122
 123
 124    TH1F* hITSSAVertexXHighMult = new TH1F("hITSSAVertexXHighMult","ITSSAVertex x High Mult; x vertex [cm] Mult>1500; events",500,-0.5,0.5);
 125    fOutput->Add(hITSSAVertexXHighMult);
 126    TH1F* hITSSAVertexYHighMult = new TH1F("hITSSAVertexYHighMult","ITSSAVertex y High Mult; y vertex [cm] Mult>1500; events",500,-0.5,0.5);
 127    fOutput->Add(hITSSAVertexYHighMult);
 128
 129
 130
 131   PostData(1, fOutput);
 132
 133   return;
 134 }
 135
 136
 137 //________________________________________________________________________
 138 void AliMeanVertexCalibTask::UserExec(Option_t *)
 139 {
 140   // Main loop
 141   // Called for each event
 142
 143   if (!InputEvent()) {
 144     Printf("ERROR: fESD not available");
 145     return;
 146   }
 147
 148   AliESDEvent* esdE = (AliESDEvent*) InputEvent();
 149
 150   const AliMultiplicity *alimult = esdE->GetMultiplicity();
 151   Int_t ntrklets=0;
 152   if(alimult) ntrklets = alimult->GetNumberOfTracklets();
 153
 154
 155   const AliESDVertex *trkv = esdE->GetPrimaryVertexTracks();
 156   const AliESDVertex *spdv = esdE->GetPrimaryVertexSPD();
 157
 158   AliESDVertex *itsSAv = ReconstructPrimaryVertex(kFALSE);
 159
 160   if(spdv) {
 161     if(spdv->GetNContributors()>0) {
 162       TString title=spdv->GetTitle();
 163       if(title.Contains("3D")) {
 164         ((TH1F*)fOutput->FindObject("hSPDVertexX"))->Fill(spdv->GetXv());
 165         ((TH1F*)fOutput->FindObject("hSPDVertexY"))->Fill(spdv->GetYv());
 166       }
 167       ((TH1F*)fOutput->FindObject("hSPDVertexZ"))->Fill(spdv->GetZv());
 168     }
 169   }
 170
 171
 172   if(trkv) {
 173     if(trkv->GetNContributors()>0) {
 174       ((TH1F*)fOutput->FindObject("hTRKVertexX"))->Fill(trkv->GetXv());
 175       ((TH1F*)fOutput->FindObject("hTRKVertexY"))->Fill(trkv->GetYv());
 176       ((TH1F*)fOutput->FindObject("hTRKVertexZ"))->Fill(trkv->GetZv());
 177
 178       if (ntrklets>30 && ntrklets<45){
 179         ((TH1F*)fOutput->FindObject("hTRKVertexXdefMult"))->Fill(trkv->GetXv());
 180         ((TH1F*)fOutput->FindObject("hTRKVertexYdefMult"))->Fill(trkv->GetYv());
 181       }
 182
 183           if (ntrklets>1500){
 184                         ((TH1F*)fOutput->FindObject("hTRKVertexXHighMult"))->Fill(trkv->GetXv());
 185                         ((TH1F*)fOutput->FindObject("hTRKVertexYHighMult"))->Fill(trkv->GetYv());
 186                 }
 187
 188       ((TH2F*)fOutput->FindObject("hTRKVertexXZ"))->Fill(trkv->GetXv(),trkv->GetZv());
 189       ((TH2F*)fOutput->FindObject("hTRKVertexYZ"))->Fill(trkv->GetYv(),trkv->GetZv());
 190
 191     }
 192   }
 193
 194   if (itsSAv){
 195     if (itsSAv->GetNContributors()>0){
 196
 197       ((TH1F*)fOutput->FindObject("hITSSAVertexX"))->Fill(itsSAv->GetXv());
 198       ((TH1F*)fOutput->FindObject("hITSSAVertexY"))->Fill(itsSAv->GetYv());
 199       ((TH1F*)fOutput->FindObject("hITSSAVertexZ"))->Fill(itsSAv->GetZv());
 200
 201                 if (ntrklets>30 && ntrklets<45){
 202                         ((TH1F*)fOutput->FindObject("hITSSAVertexXdefMult"))->Fill(itsSAv->GetXv());
 203                     ((TH1F*)fOutput->FindObject("hITSSAVertexYdefMult"))->Fill(itsSAv->GetYv());
 204                 }
 205
 206                 if (ntrklets>1500){
 207                         ((TH1F*)fOutput->FindObject("hITSSAVertexXHighMult"))->Fill(itsSAv->GetXv());
 208                         ((TH1F*)fOutput->FindObject("hITSSAVertexYHighMult"))->Fill(itsSAv->GetYv());
 209                 }
 210
 211       ((TH2F*)fOutput->FindObject("hITSSAVertexXZ"))->Fill(itsSAv->GetXv(),itsSAv->GetZv());
 212       ((TH2F*)fOutput->FindObject("hITSSAVertexYZ"))->Fill(itsSAv->GetYv(),itsSAv->GetZv());
 213
 214     }
 215   }
 216
 217   delete itsSAv;
 218
 219   PostData(1, fOutput);
 220
 221   return;
 222
 223 }
 224
 225
 226 //________________________________________________________________________
 227 void AliMeanVertexCalibTask::Terminate(Option_t *)
 228 {
 229   // Draw result to the screen
 230   // Called once at the end of the query
 231   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 232   if (!fOutput) {
 233     Printf("ERROR: fOutput not available");
 234     return;
 235   }
 236
 237
 238   return;
 239
 240 }
 241
 242
 243 //__________________________________________________________________________
 244 AliESDVertex* AliMeanVertexCalibTask::ReconstructPrimaryVertex(Bool_t constr,Int_t mode) const {
 245   // On the fly reco of ITS+TPC vertex from ESD
 246   // mode=0 use all tracks
 247   // mode=1 use odd-number tracks
 248   // mode=2 use even-number tracks
 249
 250   AliESDEvent* evt = (AliESDEvent*) fInputEvent;
 251   AliVertexerTracks vertexer(evt->GetMagneticField());
 252   if(evt->GetNumberOfTracks()<500) {
 253     vertexer.SetITSMode(); // defaults
 254     vertexer.SetMinClusters(4); // default is 5
 255   } else {
 256     vertexer.SetITSMode(0.1,0.1,0.5,5,1,3.,100.,1000.,3.,30.,1,1);// PbPb
 257   }
 258   if (fOnlyITSSATracks) vertexer.SetITSpureSA();
 259
 260   Float_t diamondcovxy[3]; evt->GetDiamondCovXY(diamondcovxy);
 261   Double_t pos[3]={evt->GetDiamondX(),evt->GetDiamondY(),0};
 262   Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
 263   AliESDVertex *initVertex = new AliESDVertex(pos,cov,1.,1);
 264   vertexer.SetVtxStart(initVertex);
 265   delete initVertex;
 266   if(!constr) vertexer.SetConstraintOff();
 267
 268   if(fOnlyITSTPCTracks || fOnlyITSSATracks || mode>0) {
 269     Int_t iskip=0;
 270     Int_t ntracks = evt->GetNumberOfTracks();
 271     Int_t *skip = new Int_t[ntracks];
 272     for(Int_t i=0;i<ntracks;i++) skip[i]=-1;
 273     for(Int_t itr=0;itr<ntracks; itr++) {
 274       AliESDtrack* track = evt->GetTrack(itr);
 275       if(fOnlyITSTPCTracks && track->GetNcls(1)==0) { // skip ITSSA
 276         skip[iskip++]=itr;
 277         continue;
 278       }
 279       if(fOnlyITSSATracks && track->GetNcls(1)>0) { // skip ITSTPC
 280         skip[iskip++]=itr;
 281         continue;
 282       }
 283       if(mode==1 && itr%2==0) skip[iskip++]=itr;
 284       if(mode==2 && itr%2==1) skip[iskip++]=itr;
 285     }
 286     vertexer.SetSkipTracks(iskip,skip);
 287     delete [] skip; skip=NULL;
 288   }
 289
 290   return vertexer.FindPrimaryVertex(evt);
 291 }