TRD/AliTRDAnalysisTaskTP.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 /* $Id: AliTRDAnalysisTaskTP.cxx 42548 2010-07-27 08:10:51Z cblume $ */
  17
  18 ////////////////////////////////////////////////////////////////////////////
  19 //                                                                        //
  20 //  Track point maker for the alignment of TRD                            //
  21 //                                                                        //
  22 //  Author:                                                               //
  23 //     Sebastian Huber (S.Huber@gsi.de)                                   //
  24 //                                                                        //
  25 ////////////////////////////////////////////////////////////////////////////
  26
  27 #include <iostream>
  28
  29 #include <TROOT.h>
  30 #include "TH1D.h"
  31 #include "TH2D.h"
  32 #include "TChain.h"
  33 #include <TLinearFitter.h>
  34
  35 #include "AliVEvent.h"
  36 #include "AliESDEvent.h"
  37 #include "AliInputEventHandler.h"
  38 #include "AliESDInputHandler.h"
  39 #include "AliAnalysisTaskSE.h"
  40 #include "AliAnalysisManager.h"
  41 #include "AliAlignObjParams.h"
  42 #include "AliTrackPointArray.h"
  43 #include "AliESDfriend.h"
  44 #include "AliESDtrack.h"
  45 #include "AliESDfriendTrack.h"
  46 #include "TFile.h"
  47 #include "TTree.h"
  48 #include "TSystem.h"
  49 #include "TTimeStamp.h"
  50 #include "TCollection.h"
  51 #include "AliLog.h"
  52 #include "AliGeomManager.h"
  53
  54 #include "AliTRDAnalysisTaskTP.h"
  55
  56 ClassImp(AliTRDAnalysisTaskTP)
  57
  58 AliTRDAnalysisTaskTP::AliTRDAnalysisTaskTP()
  59   :AliAnalysisTaskSE(),
  60   fArrHists(0x0),
  61   fArrTTree(0x0),
  62   fTree(NULL),
  63   fESD(0x0),
  64   fModpop(0x0),
  65   fBug(0x0),
  66   fNevents(0x0),
  67   fNtracks(0x0),
  68   fNAcceptedTracks(0),
  69   fArray(0x0),
  70   fFile(0x0)
  71 {
  72   // default constructor
  73
  74 }
  75 //____________________________________________________________
  76 AliTRDAnalysisTaskTP::AliTRDAnalysisTaskTP(const char *name) :
  77   AliAnalysisTaskSE(name),
  78   fArrHists(0x0),
  79   fArrTTree(0x0),
  80   fTree(0x0),
  81   fESD(0x0),
  82   fModpop(0x0),
  83   fBug(0x0),
  84   fNevents(0),
  85   fNtracks(0),
  86   fNAcceptedTracks(0),
  87   fArray(0x0),
  88   fFile(0x0)
  89 {
  90   // Constructor
  91   DefineOutput(1, TTree::Class());
  92   DefineOutput(2, TObjArray::Class());
  93
  94 }
  95 //____________________________________________________________
  96 AliTRDAnalysisTaskTP::~AliTRDAnalysisTaskTP() {
  97   // destructor
  98
  99 }
 100 //____________________________________________________________
 101 void AliTRDAnalysisTaskTP::UserCreateOutputObjects() {
 102
 103   AliAlignObjParams alobj;  // initialize align obj.
 104   TString option = GetOption();
 105
 106 if (!fArrHists) fArrHists=new TObjArray;
 107
 108   fModpop = new TH2D("modpop","modpop",90,-0.5,89.5,30,-0.5,29.5);
 109   fModpop->SetXTitle("module nr");
 110   fModpop->SetYTitle("layer nr");
 111   fArrHists->Add(fModpop);
 112
 113   OpenFile(1);
 114   fTree = new TTree("spTree", "Tree with track space point arrays");
 115   fTree->Branch("SP","AliTrackPointArray", &fArray);
 116
 117 }
 118 //____________________________________________________________
 119 void AliTRDAnalysisTaskTP::UserExec(Option_t *) {
 120
 121   //AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
 122   fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
 123   if(!fESD){
 124   //cout << "ERROR: fESDs not available " << endl;
 125   return;
 126   }
 127
 128   fESDfriend = dynamic_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
 129   if (!fESDfriend) {
 130   //cout << "ERROR: fESDfriends not available " << endl;
 131   return;
 132   }
 133
 134   TLinearFitter fitter(2, "pol1");
 135   TLinearFitter fitterz(2, "pol1");
 136
 137   // track cuts
 138   int tpc = 0; // require tpc
 139   int ptu = 0; // require certain pt's (magnetic field and tpc presumably on)
 140
 141   const Float_t kMaxDelta       = 1;
 142   const Float_t kMinNcl         = 60;
 143   const Float_t kMinPtLow       = 0.2;
 144   const Float_t kMinNclLow      = 100;
 145   const Float_t kMinPt0         = 2;
 146   const Float_t kMinPt          = 0;
 147   UInt_t status = AliESDtrack::kTRDrefit;
 148   if (tpc) status |= AliESDtrack::kTPCrefit;
 149
 150   const Float_t kMinRadius2  = 2*2;
 151   const Float_t kMaxRadius2  = 400*400;
 152   const Float_t kDeadSpace   = 4;
 153   const Float_t kTan = TMath::Tan(10*TMath::DegToRad());
 154   Int_t ntracks = fESD->GetNumberOfTracks();
 155   const AliTrackPointArray *array=0;
 156   AliTrackPointArray *tmpArray = 0;
 157   // trdarray contains all trd points in this event, for duplication detection
 158   AliTrackPointArray *trdarray = new AliTrackPointArray(1000);
 159   int ntrdarray = 0;
 160
 161   for (Int_t itrack=0; itrack < ntracks; itrack++) { //track loop
 162
 163     AliESDtrack * track = fESD->GetTrack(itrack);
 164     fNtracks++;
 165     if (!track) continue;
 166
 167     if (track->GetP() < kMinPt) continue;
 168     if (track->GetKinkIndex(0)!=0) continue;
 169     if (tpc) if (track->GetTPCNcls()<kMinNcl) continue;
 170     if (ptu) if (track->GetP() < kMinPtLow) continue;
 171     if (ptu) if (track->GetP() < kMinPt0 && track->GetTPCNcls()<kMinNclLow) continue;
 172
 173     AliESDfriendTrack *friendtrack = fESDfriend->GetTrack(itrack);
 174     if (!friendtrack){
 175     continue;
 176     }
 177
 178     array = friendtrack->GetTrackPointArray();
 179     if (!array) continue;
 180     Int_t npoints = array->GetNPoints();
 181     if (tmpArray) delete tmpArray;
 182     tmpArray = new AliTrackPointArray(npoints);
 183     Int_t current = 0;
 184     int ntpc = 0;
 185     int ntrd = 0;
 186     for (Int_t ipoint=0; ipoint<npoints; ipoint++){
 187
 188       AliTrackPoint p;
 189       array->GetPoint(p, ipoint);
 190
 191       UShort_t volid = array->GetVolumeID()[ipoint];
 192       Int_t iModule;
 193       AliGeomManager::ELayerID layer = AliGeomManager::VolUIDToLayer(volid,iModule);
 194       if ((layer < AliGeomManager::kFirstLayer) || (layer >= AliGeomManager::kLastLayer)) continue;
 195
 196       if ((iModule >= AliGeomManager::LayerSize(layer)) || (iModule < 0)) continue;
 197
 198       Float_t r2 = p.GetX()*p.GetX()+p.GetY()*p.GetY();
 199       if ( r2<kMinRadius2 || r2 > kMaxRadius2 ) continue;
 200
 201
 202
 203       if (layer>=AliGeomManager::kTPC1 && layer<=AliGeomManager::kTPC2){
 204         if (p.GetCov()[0]<0 || p.GetCov()[3]<0 ||  p.GetCov()[5]<0) continue;
 205
 206         AliTrackPoint& plocal = p.MasterToLocal();
 207         Double_t ylocal  = plocal.GetY();
 208         Double_t zlocal  = plocal.GetZ();
 209         Double_t xlocal  = plocal.GetX();
 210         Float_t edgey = TMath::Abs(plocal.GetX()*kTan);
 211         Int_t nclose=0;
 212         fitter.ClearPoints();
 213         fitterz.ClearPoints();
 214         for (Int_t jpoint=ipoint-7; jpoint<=ipoint+7; jpoint++){
 215           if (jpoint<0 || jpoint>=npoints) continue;
 216           if (ipoint==jpoint) continue;
 217           UShort_t volidL = array->GetVolumeID()[jpoint];
 218           if (volidL!=volid) continue;
 219           AliTrackPoint pc;
 220           array->GetPoint(pc, jpoint);
 221           AliTrackPoint &pcl=  pc.MasterToLocal();
 222           Double_t dx = pcl.GetX()-xlocal;
 223           fitter.AddPoint(&dx,pcl.GetY(),1);
 224           fitterz.AddPoint(&dx,pcl.GetZ(),1);
 225           nclose++;
 226         }
 227         if (nclose<6) continue;
 228
 229         fitter.Eval();
 230         fitterz.Eval();
 231         Double_t fity =fitter.GetParameter(0);
 232         Double_t fitz =fitterz.GetParameter(0);
 233         if (TMath::Abs(ylocal-fity)>kMaxDelta) continue;
 234         if (TMath::Abs(zlocal-fitz)>kMaxDelta) continue;
 235         if (TMath::Abs(fity)>edgey-kDeadSpace) continue;
 236         ntpc++;
 237       }
 238
 239       if (layer>=AliGeomManager::kTRD1 && layer<=AliGeomManager::kTRD6){
 240
 241         trdarray->AddPoint(ntrdarray++,&p);
 242         fModpop->Fill(iModule,layer);
 243         ntrd++;
 244       }
 245       tmpArray->AddPoint(current,&p);
 246       current++;
 247     }
 248     if (ntpc < 100) continue;
 249     if (ntrd < 4) continue;
 250     if (fArray) delete fArray;
 251     fArray = new AliTrackPointArray(current);
 252     for (Int_t ipoint=0; ipoint<current; ipoint++){
 253       AliTrackPoint p;
 254       tmpArray->GetPoint(p, ipoint);
 255       fArray->AddPoint(ipoint,&p);
 256     }
 257     fNAcceptedTracks++;
 258     fTree->Fill();
 259   }
 260   delete trdarray;
 261   fNevents++;
 262   PostData(1,fTree);
 263   PostData(2,fArrHists);
 264 }
 265 //____________________________________________________________
 266 void AliTRDAnalysisTaskTP::Terminate(Option_t */*option*/) {
 267   //
 268   // Terminate
 269   //
 270 }
 271
 272