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