Compilation on Windows/Cygwin
[u/mrichter/AliRoot.git] / PWG0 / TPC / AliROCESDAnalysisSelector.cxx
CommitLineData
df71af87 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
1d7991a5 18//
19// This class analyses TPC cosmics data from the ESD and the ESDfriend
df71af87 20//
1d7991a5 21// Authors: Jan.Fiete.Grosse-Oetringhaus@cern.ch, Claus.Jorgensen@cern.ch
df71af87 22//
df71af87 23
24#include "AliROCESDAnalysisSelector.h"
25
26#include <AliLog.h>
27#include <AliESD.h>
28#include <AliESDfriend.h>
dc89d87e 29#include <../TPC/AliTPCclusterMI.h>
30#include <../TPC/AliTPCseed.h>
df71af87 31
df71af87 32#include <TFile.h>
9ecad4f3 33#include <TMath.h>
df71af87 34#include <TTree.h>
35#include <TCanvas.h>
899625a7 36#include <TSystem.h>
9ecad4f3 37#include <TObjArray.h>
899625a7 38#include <TTimeStamp.h>
2d9e89d4 39
40#include "TPC/AliTPCClusterHistograms.h"
df71af87 41
899625a7 42extern TSystem* gSystem;
43
df71af87 44ClassImp(AliROCESDAnalysisSelector)
45
46AliROCESDAnalysisSelector::AliROCESDAnalysisSelector() :
47 AliSelector(),
9ecad4f3 48 fESDfriend(0),
52a28f6e 49 fObjectsToSave(0),
50 fMinNumberOfRowsIsTrack(0)
df71af87 51{
52 //
53 // Constructor. Initialization of pointers
54 //
9ecad4f3 55 fMinNumberOfRowsIsTrack = 5;
56
57 fObjectsToSave = new TObjArray();
9cc7192c 58
1d7991a5 59 for (Int_t i=0; i<kTPCHists; i++)
9cc7192c 60 fClusterHistograms[i] = 0;
df71af87 61}
62
63AliROCESDAnalysisSelector::~AliROCESDAnalysisSelector()
64{
65 //
66 // Destructor
67 //
68}
69
70void AliROCESDAnalysisSelector::SlaveBegin(TTree* tree)
71{
72 //
73
74 AliSelector::SlaveBegin(tree);
df71af87 75}
76
77void 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);
1d7991a5 86
87 printf("Init called %p\n", (void*) fESDfriend);
df71af87 88
89 // Set branch address
1d7991a5 90 if (tree)
91 {
92 tree->SetBranchAddress("ESDfriend", &fESDfriend);
93
6c3d6245 94 tree->SetBranchStatus("*", 0);
95 tree->SetBranchStatus("fTracks.*", 1);
96 tree->SetBranchStatus("fTimeStamp", 1);
1d7991a5 97 //tree->SetBranchStatus("fTracks.fCalibContainer", 0);
6c3d6245 98 }
1d7991a5 99
df71af87 100 if (fESDfriend != 0)
101 AliDebug(AliLog::kInfo, "INFO: Found ESDfriend branch in chain.");
102}
103
104Bool_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 }
dc89d87e 128
129 if (fESD->GetNumberOfTracks() != fESDfriend->GetNumberOfTracks())
130 {
52a28f6e 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()));
dc89d87e 132 return kFALSE;
133 }
134
df71af87 135 fESD->SetESDfriend(fESDfriend);
136
899625a7 137 Int_t flag = ProcessEvent(entry, kFALSE);
138 if (flag == 1)
139 ProcessEvent(entry, kTRUE);
9ecad4f3 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;
dc89d87e 146 fESD = 0;
9ecad4f3 147
148 //delete fESDfriend;
149 //fESDfriend = 0;
150
9ecad4f3 151 return kTRUE;
152}
153
899625a7 154Int_t AliROCESDAnalysisSelector::ProcessEvent(Long64_t entry, Bool_t detailedHistogram)
9ecad4f3 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
899625a7 164 // save maximum 50 objects
9ecad4f3 165 if (detailedHistogram)
899625a7 166 if (fObjectsToSave->GetEntries() > 50)
9ecad4f3 167 return 0;
899625a7 168
9ecad4f3 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
df71af87 176 Int_t nTracks = fESD->GetNumberOfTracks();
177
9cc7192c 178 Int_t nSkippedSeeds = 0;
52a28f6e 179 //Int_t nSkippedTracks = 0;
9ecad4f3 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;
9cc7192c 186
df71af87 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 {
9cc7192c 207 AliDebug(AliLog::kDebug, Form("ERROR: Could not retrieve seed of track %d.", t));
208 nSkippedSeeds++;
df71af87 209 continue;
210 }
211
52a28f6e 212 /*if (!AcceptTrack(seed, fMinNumberOfRowsIsTrack))
9ecad4f3 213 {
214 AliDebug(AliLog::kDebug, Form("INFO: Rejected track %d.", t));
215 nSkippedTracks++;
0b3ccaa2 216 continue;
52a28f6e 217 }*/
9ecad4f3 218
df71af87 219 for (Int_t clusterID = 0; clusterID < 160; clusterID++)
220 {
221 AliTPCclusterMI* cluster = seed->GetClusterPointer(clusterID);
9ecad4f3 222 if (!cluster)
df71af87 223 continue;
df71af87 224
9cc7192c 225 Int_t detector = cluster->GetDetector();
226
899625a7 227 if (detector < 0 || detector >= kTPCSectors)
228 {
229 AliDebug(AliLog::kDebug, Form("We found a cluster from invalid sector %d", detector));
230 continue;
9cc7192c 231 }
483b1eb8 232
9ecad4f3 233 if (!detailedHistogram) {
899625a7 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 }
9ecad4f3 254 } // end of if !detailedHistograms
255 else {
899625a7 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);
9ecad4f3 261 }
df71af87 262 }
9ecad4f3 263
264 for (Int_t i=0; i<kTPCHists; i++)
265 if (fClusterHistograms[i])
899625a7 266 fClusterHistograms[i]->FillTrack(seed);
267
df71af87 268 }
9cc7192c 269
9ecad4f3 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) {
899625a7 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;
9ecad4f3 276 }
277 }
278 }
279 else {
280 for (Int_t i=0; i< kTPCSectors; i++) {
281 if (clusterHistograms[i]) {
899625a7 282 fObjectsToSave->Add(clusterHistograms[i]);
9ecad4f3 283 }
284 }
285 }
1d7991a5 286
9ecad4f3 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);
df71af87 291
9ecad4f3 292 return intToReturn;
df71af87 293}
294
9ecad4f3 295Bool_t AliROCESDAnalysisSelector::AcceptTrack(const AliTPCseed* track, Int_t minRowsIncluded) {
0b3ccaa2 296 //
9ecad4f3 297 // check if the track should be accepted.
0b3ccaa2 298 //
9ecad4f3 299 const Int_t kMinClusters = 5;
0b3ccaa2 300 const Float_t kMinRatio = 0.75;
301 const Float_t kMax1pt = 0.5;
302
9ecad4f3 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;
0b3ccaa2 333
9ecad4f3 334 if (nRowsUsedByTracks < minRowsIncluded)
335 return kFALSE;
336
337 // printf(Form(" TRACK: n clusters = %d, n pad rows = %d \n",track->GetNumberOfClusters(), nRowsUsedByTracks));
338
0b3ccaa2 339 if (track->GetNumberOfClusters()<kMinClusters) return kFALSE;
340 Float_t ratio = track->GetNumberOfClusters()/(track->GetNFoundable()+1.);
341 if (ratio<kMinRatio) return kFALSE;
f4df0e80 342 Float_t mpt = track->GetSigned1Pt();
0b3ccaa2 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
df71af87 352void AliROCESDAnalysisSelector::SlaveTerminate()
353{
354 //
355
1d7991a5 356 if (fOutput)
357 {
358 for (Int_t i=0; i<kTPCHists; i++)
359 if (fClusterHistograms[i])
9cc7192c 360 fOutput->Add(fClusterHistograms[i]);
1d7991a5 361 }
df71af87 362}
363
364void AliROCESDAnalysisSelector::Terminate()
365{
1d7991a5 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"));
6c3d6245 386
387 if (comment)
1d7991a5 388 {
899625a7 389 AliDebug(AliLog::kInfo, Form("INFO: Found comment in input list: %s", comment->GetTitle()));
1d7991a5 390 }
391 else
392 return;
6c3d6245 393
394 TFile* file = TFile::Open(Form("rocESD_%s.root",comment->GetTitle()), "RECREATE");
899625a7 395
1d7991a5 396 for (Int_t i=0; i<kTPCHists; i++)
483b1eb8 397 if (fClusterHistograms[i]) {
9cc7192c 398 fClusterHistograms[i]->SaveHistograms();
52a28f6e 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;
483b1eb8 408 }
9ecad4f3 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
df71af87 426 file->Close();
899625a7 427}