TRD/Macros/AliTRDpointMaker.C

   1 #define PointMaker_cxx
   2 // The class definition in esdTree.h has been generated automatically
   3 // by the ROOT utility TTree::MakeSelector(). This class is derived
   4 // from the ROOT class TSelector. For more information on the TSelector
   5 // framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.
   6
   7 // The following methods are defined in this file:
   8 //    Begin():        called everytime a loop on the tree starts,
   9 //                    a convenient place to create your histograms.
  10 //    SlaveBegin():   called after Begin(), when on PROOF called only on the
  11 //                    slave servers.
  12 //    Process():      called for each event, in this function you decide what
  13 //                    to read and fill your histograms.
  14 //    SlaveTerminate: called at the end of the loop on the tree, when on PROOF
  15 //                    called only on the slave servers.
  16 //    Terminate():    called at the end of the loop on the tree,
  17 //                    a convenient place to draw/fit your histograms.
  18 //
  19
  20 // To do: rewrite as AliAnalysisTaskSE.(like initiated by Alex at some point)
  21
  22 #include <TROOT.h>
  23 #include <TStyle.h>
  24 #include <TLinearFitter.h>
  25 #include "AliAlignObjParams.h"
  26 #include "AliTrackPointArray.h"
  27 #include "AliLog.h"
  28 #include "PointMaker.h"
  29 //=============================================================================
  30 PointMaker::PointMaker(char *outfil) :
  31   TSelector(),
  32   fChain(0),
  33   fESD(0),
  34   //  fESDfriend(0),
  35   fFile(0),
  36   fTree(0),
  37   fArray(0),
  38   fNevents(0),
  39   fNtracks(0),
  40   fNAcceptedTracks(0),
  41   fOutfil(outfil)
  42 {
  43   // Constructor. Initialization of pointers
  44 }
  45 //=============================================================================
  46 PointMaker::~PointMaker() {
  47   // Remove all pointers
  48
  49   //  delete fESD;
  50
  51   // histograms are in the output list and deleted when the output
  52   // list is deleted by the TSelector dtor
  53   delete fFile;
  54 }
  55 //=============================================================================
  56 void PointMaker::Begin(TTree *)
  57 {
  58   // The Begin() function is called at the start of the query.
  59   // When running with PROOF Begin() is only called on the client.
  60   // The tree argument is deprecated (on PROOF 0 is passed).
  61
  62   TString option = GetOption();
  63
  64   // create output file and tree  - TEMPORARY -until root trees are merged in the memory
  65   fFile = TFile::Open(fOutfil.Data(), "RECREATE",0);
  66   fTree = new TTree("spTree", "Tree with track space point arrays");
  67   fTree->Branch("SP","AliTrackPointArray", &fArray);
  68   printf("Begin called\n");
  69 }
  70 //=============================================================================
  71 void PointMaker::SlaveBegin(TTree * tree)
  72 {
  73   // The SlaveBegin() function is called after the Begin() function.
  74   // When running with PROOF SlaveBegin() is called on each slave server.
  75   // The tree argument is deprecated (on PROOF 0 is passed).
  76
  77   Init(tree);
  78
  79   AliAlignObjParams alobj;  // initialize align obj.
  80   TString option = GetOption();
  81
  82   // histograms to monitor cut efficiency and module population
  83   fCuttra = new TH1D("cuttra","cuttra",20,0.5,20.5);
  84   fCutpoi = new TH1D("cutpoi","cutpoi",20,0.5,20.5);
  85   fModpop = new TH2D("modpop","modpop",90,-0.5,89.5,30,-0.5,29.5);
  86   fModpop->SetXTitle("module nr");
  87   fModpop->SetYTitle("layer nr");
  88   printf("SlaveBegin called\n");
  89 }
  90 //=============================================================================
  91 void PointMaker::Init(TTree *tree)
  92 {
  93   // The Init() function is called when the selector needs to initialize
  94   // a new tree or chain. Typically here the branch addresses of the tree
  95   // will be set. It is normaly not necessary to make changes to the
  96   // generated code, but the routine can be extended by the user if needed.
  97   // Init() will be called many times when running with PROOF.
  98
  99   // Set branch addresses
 100   if (tree == 0) return;
 101   fChain = tree;
 102
 103   /*
 104   fChain->SetBranchAddress("ESD",&fESD);
 105   fChain->SetBranchStatus("ESDfriend*",1);
 106   fChain->SetBranchAddress("ESDfriend.",&fESDfriend);
 107   */
 108
 109   fESD = new AliESDEvent();
 110   fChain->SetBranchStatus("ESDfriend*",1);
 111   fESD->ReadFromTree(fChain);
 112   fESDfriend = (AliESDfriend*)fESD->FindListObject("AliESDfriend");
 113   if(!fESDfriend) fChain->SetBranchAddress("ESDfriend.",&fESDfriend);
 114   printf("Init called\n");
 115 }
 116 //=============================================================================
 117 Bool_t PointMaker::Notify()
 118 {
 119   // The Notify() function is called when a new file is opened. This
 120   // can be either for a new TTree in a TChain or when when a new TTree
 121   // is started when using PROOF. Typically here the branch pointers
 122   // will be retrieved. It is normaly not necessary to make changes
 123   // to the generated code, but the routine can be extended by the
 124   // user if needed.
 125
 126   printf("Notify called\n");
 127   return kTRUE;
 128 }
 129 //=============================================================================
 130 Bool_t PointMaker::IsIdenticalWithOneOf(AliTrackPoint *p, AliTrackPointArray *parray, int nmax) {
 131
 132   // Is the point p identical with one of the points on the list parray?
 133   // This is a fix for aliroot 4-16-Rev-01 (and before) writing some
 134   // spurious uninitialized points.
 135
 136   for (int i=0; i<parray->GetNPoints() && i<nmax; i++) {
 137     AliTrackPoint pa;
 138     parray->GetPoint(pa,i);
 139     //printf("comparing %7.3f with %7.3f\n",p->GetY(),pa.GetY());
 140     if (p->GetResidual(pa,0)<1e-8) return kTRUE;
 141     //printf("different\n");
 142   }
 143   return kFALSE;
 144 }
 145 //=============================================================================
 146 Bool_t PointMaker::Process(Long64_t entry)
 147 {
 148   // The Process() function is called for each entry in the tree (or possibly
 149   // keyed object in the case of PROOF) to be processed. The entry argument
 150   // specifies which entry in the currently loaded tree is to be processed.
 151   // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
 152   // to read either all or the required parts of the data. When processing
 153   // keyed objects with PROOF, the object is already loaded and is available
 154   // via the fObject pointer.
 155   //
 156   // This function should contain the "body" of the analysis. It can contain
 157   // simple or elaborate selection criteria, run algorithms on the data
 158   // of the event and typically fill histograms.
 159
 160   // WARNING when a selector is used with a TChain, you must use
 161   //  the pointer to the current TTree to call GetEntry(entry).
 162   //  The entry is always the local entry number in the current tree.
 163   //  Assuming that fChain is the pointer to the TChain being processed,
 164   //  use fChain->GetTree()->GetEntry(entry).
 165
 166   fChain->GetTree()->GetEntry(entry);
 167   fESD->SetESDfriend(fESDfriend);
 168
 169   TLinearFitter fitter(2, "pol1");
 170   TLinearFitter fitterz(2, "pol1");
 171
 172   // track cuts
 173
 174   int tpc = 0; // require tpc
 175   int ptu = 0; // require certain pt's (magnetic field and tpc presumably on)
 176
 177   const Float_t kMaxDelta       = 1;
 178   const Float_t kMinNcl         = 60;
 179   const Float_t kMinPtLow       = 0.2;
 180   const Float_t kMinNclLow      = 100;
 181   const Float_t kMinPt0         = 2;
 182   UInt_t status = AliESDtrack::kTRDrefit;
 183   if (tpc) status |= AliESDtrack::kTPCrefit;
 184
 185   const Float_t kMinRadius2  = 2*2;
 186   const Float_t kMaxRadius2  = 400*400;
 187   const Float_t kDeadSpace   = 4;
 188   const Float_t kTan = TMath::Tan(10*TMath::DegToRad());
 189   Int_t ntracks = fESD->GetNumberOfTracks();
 190   const AliTrackPointArray *array = 0;
 191   AliTrackPointArray *tmpArray = 0;
 192   // trdarray contains all trd points in this event, for duplication detection
 193   AliTrackPointArray *trdarray = new AliTrackPointArray(1000);
 194   int ntrdarray = 0;
 195
 196   int pr = (!gROOT->IsBatch());
 197   for (Int_t itrack=0; itrack < ntracks; itrack++) {
 198     if (pr) printf("\revent %d   track %d",fNevents,itrack);
 199     fCuttra->Fill(1);
 200     AliESDtrack * track = fESD->GetTrack(itrack);
 201     fNtracks++;
 202     if (!track) continue;
 203     fCuttra->Fill(2);
 204     //if ((track->GetStatus() & status) == 0) continue;
 205     fCuttra->Fill(3);
 206     if (track->GetKinkIndex(0)!=0) continue;
 207     fCuttra->Fill(4);
 208     if (tpc) if (track->GetTPCNcls()<kMinNcl) continue;
 209     fCuttra->Fill(5);
 210     if (ptu) if (track->GetP() < kMinPtLow) continue;
 211     fCuttra->Fill(6);
 212     if (ptu) if (track->GetP() < kMinPt0 && track->GetTPCNcls()<kMinNclLow) continue;
 213     fCuttra->Fill(7);
 214     //
 215     // select points
 216     //
 217     array = track->GetTrackPointArray();
 218     if (!array) continue;
 219     Int_t npoints = array->GetNPoints();
 220     if (tmpArray) delete tmpArray;
 221     tmpArray = new AliTrackPointArray(npoints);
 222     Int_t current = 0;
 223     int ntpc = 0; // number of good TPC points
 224     int ntrd = 0; // number of good TRD points
 225     for (Int_t ipoint=0; ipoint<npoints; ipoint++){
 226       fCutpoi->Fill(1);
 227       AliTrackPoint p;
 228       array->GetPoint(p, ipoint);
 229
 230       // filter buggy points
 231
 232       UShort_t volid = array->GetVolumeID()[ipoint];
 233       Int_t iModule;
 234       AliGeomManager::ELayerID layer = AliGeomManager::VolUIDToLayer(volid,iModule);
 235       if ((layer < AliGeomManager::kFirstLayer) || (layer >= AliGeomManager::kLastLayer)) continue;
 236       fCutpoi->Fill(2);
 237       if ((iModule >= AliGeomManager::LayerSize(layer)) || (iModule < 0)) continue;
 238       fCutpoi->Fill(3);
 239       Float_t r2 = p.GetX()*p.GetX()+p.GetY()*p.GetY();
 240       if ( r2<kMinRadius2 || r2 > kMaxRadius2 ) continue;
 241       fCutpoi->Fill(4);
 242
 243       // TPC clusters handling
 244
 245       if (layer>=AliGeomManager::kTPC1 && layer<=AliGeomManager::kTPC2){
 246         // invalid covariance
 247         if (p.GetCov()[0]<0 || p.GetCov()[3]<0 ||  p.GetCov()[5]<0) continue;
 248         fCutpoi->Fill(7);
 249
 250         // remove edge clusters
 251
 252         AliTrackPoint& plocal = p.MasterToLocal();
 253         Double_t ylocal  = plocal.GetY();
 254         Double_t zlocal  = plocal.GetZ();
 255         Double_t xlocal  = plocal.GetX();
 256         Float_t edgey = TMath::Abs(plocal.GetX()*kTan);
 257         Int_t nclose=0;
 258         fitter.ClearPoints();
 259         fitterz.ClearPoints();
 260         for (Int_t jpoint=ipoint-7; jpoint<=ipoint+7; jpoint++){
 261           if (jpoint<0 || jpoint>=npoints) continue;
 262           if (ipoint==jpoint) continue;
 263           UShort_t volidL = array->GetVolumeID()[jpoint];
 264           if (volidL!=volid) continue;
 265           AliTrackPoint pc;
 266           array->GetPoint(pc, jpoint);
 267           AliTrackPoint &pcl=  pc.MasterToLocal();
 268           Double_t dx = pcl.GetX()-xlocal;
 269           fitter.AddPoint(&dx,pcl.GetY(),1);
 270           fitterz.AddPoint(&dx,pcl.GetZ(),1);
 271           nclose++;
 272         }
 273         if (nclose<6) continue;
 274         fCutpoi->Fill(8);
 275         fitter.Eval();
 276         fitterz.Eval();
 277         Double_t fity =fitter.GetParameter(0);
 278         Double_t fitz =fitterz.GetParameter(0);
 279         if (TMath::Abs(ylocal-fity)>kMaxDelta) continue;  //outlier
 280         fCutpoi->Fill(9);
 281         if (TMath::Abs(zlocal-fitz)>kMaxDelta) continue;  //outlier
 282         fCutpoi->Fill(10);
 283         if (TMath::Abs(fity)>edgey-kDeadSpace) continue;   // remove edge clusters
 284         fCutpoi->Fill(11);
 285         ntpc++;
 286       }
 287       fCutpoi->Fill(12);
 288
 289       // TRD track points handling
 290
 291       if (layer>=AliGeomManager::kTRD1 && layer<=AliGeomManager::kTRD6){
 292         fCutpoi->Fill(14);
 293         if (IsIdenticalWithOneOf(&p,trdarray,ntrdarray)) continue; // bug fix
 294         trdarray->AddPoint(ntrdarray++,&p);
 295         fCutpoi->Fill(15);
 296         fModpop->Fill(iModule,layer);
 297         ntrd++;
 298       }
 299       //printf("event %4d  track %4d  volid %4d  layer %4d  module %4d  %10.3f %10.3f %10.3f\n",
 300       //     (int)entry,itrack,volid,layer,iModule,p.GetX(),p.GetY(),atan2(p.GetY(),p.GetX())*180/3.1416);
 301       tmpArray->AddPoint(current,&p);
 302       current++;
 303       fCutpoi->Fill(16);
 304     }
 305     fCuttra->Fill(8);
 306     if (ntpc < 100) continue;
 307     fCuttra->Fill(9);
 308     if (ntrd < 4) continue;
 309     fCuttra->Fill(10);
 310     if (fArray) delete fArray;
 311     fArray = new AliTrackPointArray(current);
 312     for (Int_t ipoint=0; ipoint<current; ipoint++){
 313       AliTrackPoint p;
 314       tmpArray->GetPoint(p, ipoint);
 315       fArray->AddPoint(ipoint,&p);
 316     }
 317     fNAcceptedTracks++;
 318     fTree->Fill();
 319   }
 320   delete trdarray;
 321   fNevents++;
 322   return kTRUE;
 323 }
 324 //=============================================================================
 325 void PointMaker::SlaveTerminate()
 326 {
 327   // The SlaveTerminate() function is called after all entries or objects
 328   // have been processed. When running with PROOF SlaveTerminate() is called
 329   // on each slave server.
 330
 331   // Add the histograms to the output on each slave server
 332   AliInfo(Form("\nNumber of tracks:\tprocessed\t%d\taccepted\t%d", fNtracks, fNAcceptedTracks));
 333 }
 334 //=============================================================================
 335 void PointMaker::Terminate()
 336 {
 337   // The Terminate() function is the last function to be called during
 338   // a query. It always runs on the client, it can be used to present
 339   // the results graphically or save the results to file.
 340
 341   //  fTree = dynamic_cast<TTree*>(fOutput->FindObject("spTree"));
 342
 343   //  TFile* file = TFile::Open("AliTrackPoints.root", "RECREATE");
 344   fFile->cd();
 345   fTree->Write();
 346   fCuttra->Write();
 347   fCutpoi->Write();
 348   fModpop->Write();
 349   fModpop->Draw();
 350 }
 351 //=============================================================================