next50 trigger mask in AliHLTGlobalEsdConverterComponent
[u/mrichter/AliRoot.git] / 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 //=============================================================================