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 }