redoing of TPC vertex when flag is set
[u/mrichter/AliRoot.git] / PWG0 / TPC / AliROCESDAnalysisSelector.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$ */
17
18 // 
19 //  This class analyses TPC cosmics data from the ESD and the ESDfriend
20 //
21 //  Authors: Jan.Fiete.Grosse-Oetringhaus@cern.ch, Claus.Jorgensen@cern.ch
22 //
23
24 #include "AliROCESDAnalysisSelector.h"
25
26 #include <AliLog.h>
27 #include <AliESD.h>
28 #include <AliESDfriend.h>
29 #include <../TPC/AliTPCclusterMI.h>
30 #include <../TPC/AliTPCseed.h>
31
32 #include <TFile.h>
33 #include <TMath.h>
34 #include <TTree.h>
35 #include <TCanvas.h>
36 #include <TSystem.h>
37 #include <TObjArray.h>
38 #include <TTimeStamp.h>
39
40 #include "TPC/AliTPCClusterHistograms.h"
41
42 extern TSystem* gSystem;
43
44 ClassImp(AliROCESDAnalysisSelector)
45
46 AliROCESDAnalysisSelector::AliROCESDAnalysisSelector() :
47   AliSelector(),
48   fESDfriend(0),
49   fObjectsToSave(0),
50   fMinNumberOfRowsIsTrack(0)
51 {
52   //
53   // Constructor. Initialization of pointers
54   //
55   fMinNumberOfRowsIsTrack = 5;
56
57   fObjectsToSave = new TObjArray();
58   
59   for (Int_t i=0; i<kTPCHists; i++)
60     fClusterHistograms[i] = 0;
61 }
62
63 AliROCESDAnalysisSelector::~AliROCESDAnalysisSelector()
64 {
65   //
66   // Destructor
67   //
68 }
69
70 void AliROCESDAnalysisSelector::SlaveBegin(TTree* tree)
71 {
72   //
73   
74   AliSelector::SlaveBegin(tree);
75
76
77 void AliROCESDAnalysisSelector::Init(TTree *tree)
78 {
79   // The Init() function is called when the selector needs to initialize
80   // a new tree or chain. Typically here the branch addresses of the tree
81   // will be set. It is normaly not necessary to make changes to the
82   // generated code, but the routine can be extended by the user if needed.
83   // Init() will be called many times when running with PROOF.
84
85   AliSelector::Init(tree);
86   
87   printf("Init called %p\n", (void*) fESDfriend);
88
89   // Set branch address
90   if (tree) 
91   {
92     tree->SetBranchAddress("ESDfriend", &fESDfriend);
93   
94     tree->SetBranchStatus("*", 0);
95     tree->SetBranchStatus("fTracks.*", 1);
96     tree->SetBranchStatus("fTimeStamp", 1);
97     //tree->SetBranchStatus("fTracks.fCalibContainer", 0);
98   }
99
100   if (fESDfriend != 0)
101     AliDebug(AliLog::kInfo, "INFO: Found ESDfriend branch in chain.");
102 }
103
104 Bool_t AliROCESDAnalysisSelector::Process(Long64_t entry)
105 {
106   //
107   // Implement your analysis here. Do not forget to call the parent class Process by
108   // if (AliSelector::Process(entry) == kFALSE)
109   //   return kFALSE;
110   //
111
112   if (AliSelector::Process(entry) == kFALSE)
113     return kFALSE;
114
115   // Check prerequisites
116   if (!fESD)
117   {
118     AliDebug(AliLog::kError, "ESD branch not available");
119     return kFALSE;
120   }
121
122   // Check prerequisites
123   if (!fESDfriend)
124   {
125     AliDebug(AliLog::kError, "ESDfriend branch not available");
126     return kFALSE;
127   }
128
129   if (fESD->GetNumberOfTracks() != fESDfriend->GetNumberOfTracks())
130   {
131     AliDebug(AliLog::kError, Form("Event %lld: Number of tracks differ between ESD (%d) and ESDfriend (%d)! Skipping event!\n", entry, fESD->GetNumberOfTracks(), fESDfriend->GetNumberOfTracks()));
132     return kFALSE;
133   }  
134
135   fESD->SetESDfriend(fESDfriend);
136
137   Int_t flag = ProcessEvent(entry, kFALSE);
138   if (flag == 1)
139     ProcessEvent(entry, kTRUE);
140
141   // TODO This should not be needed, the TTree::GetEntry() should take care of this, maybe because it has a reference member, to be analyzed
142   // if the ESDfriend is not deleted we get a major memory leak
143   // here the esdfriend seems to be also deleted, very weird behaviour....
144
145   delete fESD;
146   fESD = 0;
147   
148   //delete fESDfriend;
149   //fESDfriend = 0;
150
151   return kTRUE;
152 }
153
154 Int_t AliROCESDAnalysisSelector::ProcessEvent(Long64_t entry, Bool_t detailedHistogram)
155 {
156   //
157   // Looping over tracks and clusters in event and filling histograms 
158   //
159   // - if detailedHistogram = kTRUE special histograms are saved (in fObjectsToSave)
160   // - the method returns 
161   //   1 : if a "flash" is detected something special in this event
162   //   
163
164   // save maximum 50 objects
165   if (detailedHistogram) 
166     if (fObjectsToSave->GetEntries() > 50) 
167       return 0;
168       
169   // for saving single events
170   AliTPCClusterHistograms* clusterHistograms[kTPCSectors];
171   for (Int_t i=0; i<kTPCSectors; i++) 
172     clusterHistograms[i] = 0;  
173
174   Bool_t intToReturn = 0;
175
176   Int_t nTracks = fESD->GetNumberOfTracks();
177   
178   Int_t nSkippedSeeds = 0;
179   //Int_t nSkippedTracks = 0;
180
181   // for "flash" detection
182   Int_t nClusters = 0;
183   Float_t clusterQtotSumVsTime[250];  
184   for (Int_t z=0; z<250; z++)
185     clusterQtotSumVsTime[z] = 0;
186   
187   // loop over esd tracks
188   for (Int_t t=0; t<nTracks; t++)
189   {
190     AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (fESD->GetTrack(t));
191     if (!esdTrack)
192     {
193       AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", t));
194       continue;
195     }
196     
197     AliESDfriendTrack* friendtrack = const_cast<AliESDfriendTrack*> (dynamic_cast<const AliESDfriendTrack*> (esdTrack->GetFriendTrack()));
198     if (!friendtrack)
199     {
200       AliDebug(AliLog::kError, Form("ERROR: Could not retrieve friend of track %d.", t));
201       continue;
202     }
203     
204     const AliTPCseed* seed = dynamic_cast<const AliTPCseed*> (friendtrack->GetCalibObject(0));
205     if (!seed)
206     {
207       AliDebug(AliLog::kDebug, Form("ERROR: Could not retrieve seed of track %d.", t));
208       nSkippedSeeds++;
209       continue;
210     }
211     
212     /*if (!AcceptTrack(seed, fMinNumberOfRowsIsTrack))
213     {
214       AliDebug(AliLog::kDebug, Form("INFO: Rejected track %d.", t));
215       nSkippedTracks++;
216       continue;
217     }*/
218     
219     for (Int_t clusterID = 0; clusterID < 160; clusterID++)
220     {
221       AliTPCclusterMI* cluster = seed->GetClusterPointer(clusterID);
222       if (!cluster) 
223         continue;
224       
225       Int_t detector = cluster->GetDetector();
226       
227       if (detector < 0 || detector >= kTPCSectors) 
228       {
229         AliDebug(AliLog::kDebug, Form("We found a cluster from invalid sector %d", detector));
230         continue;
231       }
232
233       if (!detailedHistogram) {
234         // TODO: find a clever way to handle the time      
235         Int_t time = 0;
236         
237         if (fESD->GetTimeStamp()>1160000000)
238           time = fESD->GetTimeStamp();      
239         
240         if (!fClusterHistograms[detector])
241           fClusterHistograms[detector] = new AliTPCClusterHistograms(detector,"",time,time+5*60*60);
242         
243         if (!fClusterHistograms[detector+kTPCSectors])
244           fClusterHistograms[detector+kTPCSectors] = new AliTPCClusterHistograms(detector,"",time,time+5*60*60, kTRUE);
245         
246         fClusterHistograms[detector]->FillCluster(cluster, time);
247         fClusterHistograms[detector+kTPCSectors]->FillCluster(cluster, time);
248         
249         Int_t z = Int_t(cluster->GetZ()); 
250         if (z>=0 && z<250) {
251           nClusters++;
252           clusterQtotSumVsTime[z] += cluster->GetQ();
253         }
254       } // end of if !detailedHistograms
255       else {
256         // if we need the detailed histograms for this event
257         if (!clusterHistograms[detector])
258           clusterHistograms[detector] = new AliTPCClusterHistograms(detector, Form("flash_entry%d", entry));
259         
260         clusterHistograms[detector]->FillCluster(cluster);
261       }
262     }
263     
264     for (Int_t i=0; i<kTPCHists; i++) 
265       if (fClusterHistograms[i]) 
266         fClusterHistograms[i]->FillTrack(seed);
267
268   }
269   
270   // check if there's a very large q deposit ("flash")
271   if (!detailedHistogram) {
272     for (Int_t z=0; z<250; z++) {
273       if (clusterQtotSumVsTime[z] > 150000) {
274         printf(Form("  \n   -> Entry %lld sum of clusters at time %d is %f, ESD timestamp: %s (%d) \n \n", entry, z, clusterQtotSumVsTime[z], TTimeStamp(fESD->GetTimeStamp()).AsString(), fESD->GetTimeStamp()));
275         intToReturn = 1;
276       }
277     }
278   }
279   else {
280     for (Int_t i=0; i< kTPCSectors; i++) {
281       if (clusterHistograms[i]) {
282         fObjectsToSave->Add(clusterHistograms[i]);
283       }
284     }    
285   }
286
287 //   if (nSkippedSeeds > 0)
288 //     printf("WARNING: The seed was not found for %d out of %d tracks.\n", nSkippedSeeds, nTracks);
289 //   if (nSkippedTracks > 0)
290 //     printf("INFO: Rejected %d out of %d tracks.\n", nSkippedTracks, nTracks);
291    
292   return intToReturn;
293 }
294
295 Bool_t AliROCESDAnalysisSelector::AcceptTrack(const AliTPCseed* track, Int_t minRowsIncluded) {
296   //
297   // check if the track should be accepted.
298   //
299   const Int_t   kMinClusters = 5;
300   const Float_t kMinRatio    = 0.75;
301   const Float_t kMax1pt      = 0.5;
302
303   Int_t  nRowsUsedByTracks = 0;
304   Bool_t rowIncluded[96];
305   
306   Float_t totalQtot = 0;
307   Int_t   nClusters = 0;
308
309   for(Int_t r=0; r<96; r++) 
310     rowIncluded[r] = kFALSE;
311   
312   for (Int_t clusterID = 0; clusterID < 160; clusterID++) {
313     AliTPCclusterMI* cluster = track->GetClusterPointer(clusterID);
314     
315     if (!cluster) 
316       continue;
317     
318     Float_t qTot =   cluster->GetQ();    
319     
320     nClusters++;
321     totalQtot += qTot;
322
323     if (!rowIncluded[cluster->GetRow()]) {
324       nRowsUsedByTracks++;
325       rowIncluded[cluster->GetRow()] = kTRUE;
326     }
327   }
328
329   Float_t meanQtot = totalQtot/nClusters;
330
331   if (meanQtot<70)
332     return kFALSE;
333
334   if (nRowsUsedByTracks < minRowsIncluded)
335     return kFALSE;
336   
337   //  printf(Form("    TRACK: n clusters = %d,  n pad rows = %d \n",track->GetNumberOfClusters(), nRowsUsedByTracks));
338   
339   if (track->GetNumberOfClusters()<kMinClusters) return kFALSE;
340   Float_t ratio = track->GetNumberOfClusters()/(track->GetNFoundable()+1.);
341   if (ratio<kMinRatio) return kFALSE;
342   Float_t mpt = track->GetSigned1Pt();
343   if (TMath::Abs(mpt)>kMax1pt) return kFALSE;
344
345   //if (TMath::Abs(track->GetZ())>240.) return kFALSE;
346   //if (TMath::Abs(track->GetZ())<10.) return kFALSE;
347   //if (TMath::Abs(track->GetTgl())>0.03) return kFALSE;
348   
349   return kTRUE;
350 }
351
352 void AliROCESDAnalysisSelector::SlaveTerminate()
353 {
354   //
355   
356   if (fOutput)
357   {
358     for (Int_t i=0; i<kTPCHists; i++)
359       if (fClusterHistograms[i])
360         fOutput->Add(fClusterHistograms[i]);
361   }
362
363
364 void AliROCESDAnalysisSelector::Terminate()
365 {
366   // 
367   // read the objects from the output list and write them to a file
368   // the filename is modified by the object comment passed in the tree info or input list
369   //
370
371   if (fOutput)
372   {  
373     fOutput->Print();
374         
375     for (Int_t i=0; i<kTPCSectors; i++)
376       fClusterHistograms[i] = dynamic_cast<AliTPCClusterHistograms*> (fOutput->FindObject(AliTPCClusterHistograms::FormDetectorName(i, kFALSE)));
377     for (Int_t i=0; i<kTPCSectors; i++)
378       fClusterHistograms[kTPCSectors+i] = dynamic_cast<AliTPCClusterHistograms*> (fOutput->FindObject(AliTPCClusterHistograms::FormDetectorName(i, kTRUE)));
379   }
380   
381   TNamed* comment = 0;
382   if (fTree && fTree->GetUserInfo())
383     comment = dynamic_cast<TNamed*>(fTree->GetUserInfo()->FindObject("comment"));
384   if (!comment && fInput)
385     comment = dynamic_cast<TNamed*>(fInput->FindObject("comment"));
386
387   if (comment)
388   {
389     AliDebug(AliLog::kInfo, Form("INFO: Found comment in input list: %s", comment->GetTitle()));
390   }
391   else
392     return;
393
394   TFile* file = TFile::Open(Form("rocESD_%s.root",comment->GetTitle()), "RECREATE");
395
396   for (Int_t i=0; i<kTPCHists; i++)
397     if (fClusterHistograms[i]) {
398       fClusterHistograms[i]->SaveHistograms();
399      // TCanvas* c = fClusterHistograms[i]->DrawHistograms(comment->GetTitle());
400                         //TString dir;
401                         //dir.Form("WWW/%s", comment->GetTitle(), c->GetName());
402                         //gSystem->mkdir(dir, kTRUE);
403       //c->SaveAs(Form("%s/plots_%s_%s.eps",dir.Data(),comment->GetTitle(),c->GetName()));
404       //c->SaveAs(Form("%s/plots_%s_%s.gif",dir.Data(),comment->GetTitle(),c->GetName()));
405
406       //c->Close();
407       //delete c;
408     }
409
410   gDirectory->mkdir("saved_objects");
411   gDirectory->cd("saved_objects");
412
413   for (Int_t i=0; i<fObjectsToSave->GetSize(); i++) {
414     if (fObjectsToSave->At(i)) {
415       AliTPCClusterHistograms* clusterHistograms = dynamic_cast<AliTPCClusterHistograms*> (fObjectsToSave->At(i));
416       if (clusterHistograms)
417         clusterHistograms->SaveHistograms();
418       else
419         fObjectsToSave->At(i)->Write();
420     }
421   }
422
423   gDirectory->cd("../");
424
425
426   file->Close();
427 }