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   //
  73   // Default constructor
  74   //
  75
  76 }
  77
  78 //____________________________________________________________
  79 AliTRDAnalysisTaskTP::AliTRDAnalysisTaskTP(const char *name) :
  80   AliAnalysisTaskSE(name),
  81   fArrHists(0x0),
  82   fArrTTree(0x0),
  83   fTree(0x0),
  84   fESD(0x0),
  85   fModpop(0x0),
  86   fBug(0x0),
  87   fNevents(0),
  88   fNtracks(0),
  89   fNAcceptedTracks(0),
  90   fArray(0x0),
  91   fFile(0x0)
  92 {
  93   //
  94   // Constructor
  95   //
  96
  97   DefineOutput(1, TTree::Class());
  98   DefineOutput(2, TObjArray::Class());
  99
 100 }
 101
 102 //____________________________________________________________
 103 AliTRDAnalysisTaskTP::~AliTRDAnalysisTaskTP()
 104 {
 105   //
 106   // destructor
 107   //
 108
 109 }
 110
 111 //____________________________________________________________
 112 void AliTRDAnalysisTaskTP::UserCreateOutputObjects()
 113 {
 114   //
 115   // Create the output objects
 116   //
 117
 118   AliAlignObjParams alobj;  // initialize align obj.
 119   TString option = GetOption();
 120
 121   if (!fArrHists) fArrHists=new TObjArray;
 122
 123   fModpop = new TH2D("modpop","modpop",90,-0.5,89.5,30,-0.5,29.5);
 124   fModpop->SetXTitle("module nr");
 125   fModpop->SetYTitle("layer nr");
 126   fArrHists->Add(fModpop);
 127
 128   OpenFile(1);
 129   fTree = new TTree("spTree", "Tree with track space point arrays");
 130   fTree->Branch("SP","AliTrackPointArray", &fArray);
 131
 132 }
 133
 134 //____________________________________________________________
 135 void AliTRDAnalysisTaskTP::UserExec(Option_t *)
 136 {
 137   //
 138   // Exec function
 139   //
 140
 141   //AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
 142   fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
 143   if(!fESD){
 144   //cout << "ERROR: fESDs not available " << endl;
 145   return;
 146   }
 147
 148   fESDfriend = dynamic_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
 149   if (!fESDfriend) {
 150   //cout << "ERROR: fESDfriends not available " << endl;
 151   return;
 152   }
 153
 154   TLinearFitter fitter(2, "pol1");
 155   TLinearFitter fitterz(2, "pol1");
 156
 157   // track cuts
 158   int tpc = 0; // require tpc
 159   int ptu = 0; // require certain pt's (magnetic field and tpc presumably on)
 160
 161   const Float_t kMaxDelta       = 1;
 162   const Float_t kMinNcl         = 60;
 163   const Float_t kMinPtLow       = 0.2;
 164   const Float_t kMinNclLow      = 100;
 165   const Float_t kMinPt0         = 2;
 166   const Float_t kMinPt          = 0;
 167   UInt_t status = AliESDtrack::kTRDrefit;
 168   if (tpc) status |= AliESDtrack::kTPCrefit;
 169
 170   const Float_t kMinRadius2  = 2*2;
 171   const Float_t kMaxRadius2  = 400*400;
 172   const Float_t kDeadSpace   = 4;
 173   const Float_t kTan = TMath::Tan(10*TMath::DegToRad());
 174   Int_t ntracks = fESD->GetNumberOfTracks();
 175   const AliTrackPointArray *array=0;
 176   AliTrackPointArray *tmpArray = 0;
 177   // trdarray contains all trd points in this event, for duplication detection
 178   AliTrackPointArray *trdarray = new AliTrackPointArray(1000);
 179   int ntrdarray = 0;
 180
 181   for (Int_t itrack=0; itrack < ntracks; itrack++) { //track loop
 182
 183     AliESDtrack * track = fESD->GetTrack(itrack);
 184     fNtracks++;
 185     if (!track) continue;
 186
 187     if (track->GetP() < kMinPt) continue;
 188     if (track->GetKinkIndex(0)!=0) continue;
 189     if (tpc) if (track->GetTPCNcls()<kMinNcl) continue;
 190     if (ptu) if (track->GetP() < kMinPtLow) continue;
 191     if (ptu) if (track->GetP() < kMinPt0 && track->GetTPCNcls()<kMinNclLow) continue;
 192
 193     AliESDfriendTrack *friendtrack = fESDfriend->GetTrack(itrack);
 194     if (!friendtrack){
 195     continue;
 196     }
 197
 198     array = friendtrack->GetTrackPointArray();
 199     if (!array) continue;
 200     Int_t npoints = array->GetNPoints();
 201     if (tmpArray) delete tmpArray;
 202     tmpArray = new AliTrackPointArray(npoints);
 203     Int_t current = 0;
 204     int ntpc = 0;
 205     int ntrd = 0;
 206     for (Int_t ipoint=0; ipoint<npoints; ipoint++){
 207
 208       AliTrackPoint p;
 209       array->GetPoint(p, ipoint);
 210
 211       UShort_t volid = array->GetVolumeID()[ipoint];
 212       Int_t iModule;
 213       AliGeomManager::ELayerID layer = AliGeomManager::VolUIDToLayer(volid,iModule);
 214       if ((layer < AliGeomManager::kFirstLayer) || (layer >= AliGeomManager::kLastLayer)) continue;
 215
 216       if ((iModule >= AliGeomManager::LayerSize(layer)) || (iModule < 0)) continue;
 217
 218       Float_t r2 = p.GetX()*p.GetX()+p.GetY()*p.GetY();
 219       if ( r2<kMinRadius2 || r2 > kMaxRadius2 ) continue;
 220
 221
 222
 223       if (layer>=AliGeomManager::kTPC1 && layer<=AliGeomManager::kTPC2){
 224         if (p.GetCov()[0]<0 || p.GetCov()[3]<0 ||  p.GetCov()[5]<0) continue;
 225
 226         AliTrackPoint& plocal = p.MasterToLocal();
 227         Double_t ylocal  = plocal.GetY();
 228         Double_t zlocal  = plocal.GetZ();
 229         Double_t xlocal  = plocal.GetX();
 230         Float_t edgey = TMath::Abs(plocal.GetX()*kTan);
 231         Int_t nclose=0;
 232         fitter.ClearPoints();
 233         fitterz.ClearPoints();
 234         for (Int_t jpoint=ipoint-7; jpoint<=ipoint+7; jpoint++){
 235           if (jpoint<0 || jpoint>=npoints) continue;
 236           if (ipoint==jpoint) continue;
 237           UShort_t volidL = array->GetVolumeID()[jpoint];
 238           if (volidL!=volid) continue;
 239           AliTrackPoint pc;
 240           array->GetPoint(pc, jpoint);
 241           AliTrackPoint &pcl=  pc.MasterToLocal();
 242           Double_t dx = pcl.GetX()-xlocal;
 243           fitter.AddPoint(&dx,pcl.GetY(),1);
 244           fitterz.AddPoint(&dx,pcl.GetZ(),1);
 245           nclose++;
 246         }
 247         if (nclose<6) continue;
 248
 249         fitter.Eval();
 250         fitterz.Eval();
 251         Double_t fity =fitter.GetParameter(0);
 252         Double_t fitz =fitterz.GetParameter(0);
 253         if (TMath::Abs(ylocal-fity)>kMaxDelta) continue;
 254         if (TMath::Abs(zlocal-fitz)>kMaxDelta) continue;
 255         if (TMath::Abs(fity)>edgey-kDeadSpace) continue;
 256         ntpc++;
 257       }
 258
 259       if (layer>=AliGeomManager::kTRD1 && layer<=AliGeomManager::kTRD6){
 260
 261         trdarray->AddPoint(ntrdarray++,&p);
 262         fModpop->Fill(iModule,layer);
 263         ntrd++;
 264       }
 265       tmpArray->AddPoint(current,&p);
 266       current++;
 267     }
 268     if (ntpc < 100) continue;
 269     if (ntrd < 4) continue;
 270     if (fArray) delete fArray;
 271     fArray = new AliTrackPointArray(current);
 272     for (Int_t ipoint=0; ipoint<current; ipoint++){
 273       AliTrackPoint p;
 274       tmpArray->GetPoint(p, ipoint);
 275       fArray->AddPoint(ipoint,&p);
 276     }
 277     fNAcceptedTracks++;
 278     fTree->Fill();
 279   }
 280   delete trdarray;
 281   fNevents++;
 282   PostData(1,fTree);
 283   PostData(2,fArrHists);
 284 }
 285
 286 //____________________________________________________________
 287 void AliTRDAnalysisTaskTP::Terminate(Option_t */*option*/)
 288 {
 289   //
 290   // Terminate
 291   //
 292
 293 }
 294
 295