Simple example macro to make some analyzis for the ESD tracks - see comments (Not...

diff --git a/TPC/AnalyzeESDtracks.C b/TPC/AnalyzeESDtracks.C
new file mode 100644 (file)
index 0000000..8cbf2fa
--- /dev/null
+++ b/TPC/AnalyzeESDtracks.C
@@ -0,0 +1,265 @@
+//
+// Process ESD tracks  - 
+// Extract TPC tracks  - write them to tree
+//
+/*
+  .L AnalyzeESDtracks.C+
+  .L AliGenInfo.C+
+  AnalyzESDtracks(567);   // process tracks
+  // Tracks are written to the file "TPCtracks.root"
+  // Now yo can analyze it
+  TFile fesd("AliESDs.root");
+  TTree * treeE = (TTree*)fesd.Get("esdTree");
+  TFile f("TPCtracks.root")
+  TTree * tree =(TTree*)f.Get("Tracks");
+  AliComparisonDraw comp;
+  comp->fTree = tree;
+
+*/
+
+// 
+
+/*
+.L AnalyzeESDtracks.C+
+.L AliGenInfo.C+
+
+//AnalyzESDtracks(567);
+
+
+TFile fesd("AliESDs.root");
+TTree * treeE = (TTree*)fesd.Get("esdTree");
+TFile f("TPCtracks.root")
+TTree * tree =(TTree*)f.Get("Tracks");
+AliComparisonDraw comp;
+comp->fTree = tree;
+
+TFile fs("TPCsignal.root");
+TTree *treeB =(TTree*)fs.Get("SignalB");
+TTree *treeN =(TTree*)fs.Get("SignalN");
+TTree *treeS =(TTree*)fs.Get("Signal");
+TTree *treef =(TTree*)fs.Get("Fit");
+
+
+FitSignals(treeB,"Max-Median>100&&RMS06<2.5")
+TFile ffit("FitSignal.root");
+TTree * treeF = (TTree*)ffit->Get("Fit");
+
+
+TChain chaincl("TreeR","TreeR")
+chaincl.Add("TPC.RecPoints.root/Event0/TreeR")
+chaincl.Add("TPC.RecPoints1.root/Event1/TreeR")
+chaincl.Add("TPC.RecPoints2.root/Event2/TreeR")
+chaincl.Add("TPC.RecPoints3.root/Event3/TreeR")
+
+*/
+
+
+
+#include "TObject.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TBranch.h"
+#include "TTreeStream.h"
+#include "TEventList.h"
+#include "TCut.h"
+#include "TFitter.h"
+#include  "TMatrixD.h"
+
+#include  "AliLog.h"
+#include  "AliMagF.h"
+
+#include  "AliESD.h"
+#include  "AliESDfriend.h"
+#include  "AliESDtrack.h"
+#include  "AliTracker.h"
+#include  "AliTPCseed.h"
+#include  "AliTPCclusterMI.h"
+#include  "TTreeStream.h"
+#include  "TF1.h"
+#include  "TGraph.h"
+#include  "AliSignalProcesor.h"
+#include  "TCanvas.h"
+
+
+void AnalyzESDtracks(Int_t run){
+  //
+  // output redirect 
+  TTreeSRedirector  cstream("TPCtracks.root");
+  //
+  // dummy magnetic field
+  AliMagF mag("aaa","aaa",1,1,10);
+  AliTracker::SetFieldMap(&mag,kTRUE);
+  TFile f("AliESDs.root");
+  TTree * tree =(TTree*)f.Get("esdTree"); 
+  AliESD * esd =0;
+  tree->SetBranchAddress("ESD",&esd);
+  AliESDfriend *evf=0;
+  tree->AddFriend("esdFriendTree","AliESDfriends.root");
+  tree->SetBranchAddress("ESDfriend",&evf);
+
+  Int_t nevents = tree->GetEntries();
+  TClonesArray *clusters = new TClonesArray("AliTPCclusterMI",160);
+  for (Int_t irow=0; irow<160; irow++){
+    new ((*clusters)[irow])  AliTPCclusterMI;   // iitial dummy clusters
+  }
+  
+  for (Int_t ievent=0; ievent<nevents; ievent++){
+    tree->GetEntry(ievent);
+    if (!esd) continue;
+    if (!evf) continue;
+    esd->SetESDfriend(evf); //Attach the friend to the ESD      
+    for (Int_t itrack =0; itrack<esd->GetNumberOfTracks(); itrack++){
+      // Int_t itrack = 0;      
+      if (esd->GetTrack(itrack)->GetFriendTrack()==0) continue;
+      AliESDtrack * etrack = esd->GetTrack(itrack);
+      AliESDfriendTrack * ftrack = (AliESDfriendTrack *)esd->GetTrack(itrack)->GetFriendTrack();
+      AliTPCseed * seed =  (AliTPCseed*)(ftrack->GetCalibObject(0));
+      if (!seed) continue;  
+      for (Int_t irow=0; irow<160; irow++){
+       if (seed->GetClusterFast(irow)){
+         AliTPCclusterMI * cl = new ((*clusters)[irow])  AliTPCclusterMI(*(seed->GetClusterFast(irow)));
+         cl->SetLabel(itrack,0);
+       }
+       else{
+         AliTPCclusterMI * cl = (AliTPCclusterMI*)clusters->At(irow);
+         cl->SetX(0); cl->SetY(0); cl->SetZ(0); cl->SetQ(0); cl->SetLabel(-1,0);
+       }
+      }
+      Float_t dEdx = seed->GetdEdx();
+      Float_t dEdxI = seed->CookdEdx(0.05,0.6,0,77);
+      Float_t dEdxO = seed->CookdEdx(0.05,0.6,78,155);
+      Int_t ncl = seed->GetNumberOfClusters();
+      cstream<<"Tracks"<<
+       "Run="<<run<<
+       "Event="<<ievent<<
+       "dEdx="<<dEdx<<
+       "dEdxI="<<dEdxI<<
+       "dEdxO="<<dEdxO<<
+       "Track.="<<seed<<
+       "Etrack.="<<etrack<<
+       "Cl.="<<clusters<<
+       "\n";
+    }  
+  }
+}
+
+
+void FitSignals(TTree * treeB, TCut cut="Max-Median>150&&RMS06<2&&abs(Median-Mean09)<0.5"){
+  AliSignalProcesor proc;
+  TF1 * f1 = proc.GetAsymGauss();
+  TTreeSRedirector cstream("FitSignal.root");
+  TFile *f = cstream.GetFile();
+
+  char lname[100];
+  sprintf(lname,"Fit%s", cut.GetTitle());
+  TEventList *list = new TEventList(lname,lname);
+  sprintf(lname,">>Fit%s", cut.GetTitle());
+  treeB->Draw(lname,cut);
+  treeB->SetEventList(list);
+  for (Int_t ievent=0; ievent<list->GetN(); ievent++){
+    char ename[100];
+    sprintf(ename,"Fit%d", ievent);
+    Double_t nsample = treeB->Draw("Graph.fY-Mean09:Graph.fX","","",1,ievent);
+    Double_t * signal  = treeB->GetV1();
+    Double_t * time  = treeB->GetV2();
+    Double_t maxpos =0;
+    Double_t max = 0;
+    for (Int_t ipos = 0; ipos<nsample; ipos++){
+      if (signal[ipos]>max){
+       max    = signal[ipos];
+       maxpos = ipos;
+      }
+    }
+
+    Int_t first = TMath::Max(maxpos-10,0.);
+    Int_t last  = TMath::Min(maxpos+60, nsample);
+    //
+    f->cd();
+    TH1F his(ename,ename,last-first,first,last);
+    for (Int_t ipos=0; ipos<last-first; ipos++){
+      his.SetBinContent(ipos+1,signal[ipos+first]);
+    }
+    treeB->Draw("Sector:Row:Pad","","",1,ievent);
+    Double_t sector = treeB->GetV1()[0];
+    Double_t row    = treeB->GetV2()[0];
+    Double_t pad    = treeB->GetV3()[0];
+    //    TGraph  graph(last-first,&time[first],&signal[first]);
+    f1->SetParameters(0.75*max,maxpos,1.1,0.8,0.25,0.2);
+    //    TH1F * his = (TH1F*)graph.GetHistogram();
+    his.Fit(f1,"q");
+    his.Write(ename);
+    gPad->Clear();
+    his.Draw();
+    gPad->Update();
+    Double_t params[6];
+    for (Int_t ipar=0; ipar<6; ipar++) params[ipar] = f1->GetParameters()[ipar];
+    Double_t chi2 = TFitter::GetFitter()->Chisquare(6,params);
+    TMatrixD cov(6,6);
+    cov.SetMatrixArray(TFitter::GetFitter()->GetCovarianceMatrix());
+    //
+    // tail cancellation
+    //
+    Double_t x0[1000];
+    Double_t x1[1000];
+    Double_t x2[1000];
+    for (Int_t ipos=0; ipos<last-first; ipos++){
+      x0[ipos] = signal[ipos+first];
+    }
+    proc.TailCancelationALTRO1(x0,x1,0.85*0.339,0.09,last-first);
+    proc.TailCancelationALTRO1(x1,x2,0.85,0.789,last-first);
+    //
+    sprintf(ename,"Cancel1_%d", ievent);
+    TH1F his1(ename,ename,last-first,first,last);
+    for (Int_t ipos=0; ipos<last-first; ipos++){
+      his1.SetBinContent(ipos+1,x1[ipos]);
+    }
+    his1.Write(ename);
+    sprintf(ename,"Cancel2_%d", ievent);
+    TH1F his2(ename,ename,last-first,first,last);
+    for (Int_t ipos=0; ipos<last-first; ipos++){
+      his2.SetBinContent(ipos+1,x1[ipos]);
+    }
+    f1->SetParameters(0.75*max,maxpos,1.1,0.8,0.25,0.2);
+    his2.Fit(f1,"q");
+    his2.Write(ename);
+    Double_t params2[6];
+    for (Int_t ipar=0; ipar<6; ipar++) params2[ipar] = f1->GetParameters()[ipar];
+    Double_t chi22 = TFitter::GetFitter()->Chisquare(6,params2);    
+    TMatrixD cov2(6,6);
+    cov2.SetMatrixArray(TFitter::GetFitter()->GetCovarianceMatrix());
+
+    TGraph gr0(last-first, &time[first],x0);
+    TGraph gr1(last-first, &time[first],x1);
+    TGraph gr2(last-first, &time[first],x2);
+    //
+    cstream<<"Fit"<<
+      "Sector="<<sector<<
+      "Row="<<row<<
+      "Pad="<<pad<<
+      "First="<<first<<
+      "Max="<<max<<
+      "MaxPos="<<maxpos<<
+      "chi2="<<chi2<<
+      "chi22="<<chi22<<
+      "Cov="<<&cov<<
+      "Cov2="<<&cov2<<
+      "gr0.="<<&gr0<<
+      "gr1.="<<&gr1<<
+      "gr2.="<<&gr2<<
+      "p0="<<params[0]<<
+      "p1="<<params[1]<<
+      "p2="<<params[2]<<
+      "p3="<<params[3]<<
+      "p4="<<params[4]<<
+      "p5="<<params[5]<<
+      "p02="<<params2[0]<<
+      "p12="<<params2[1]<<
+      "p22="<<params2[2]<<
+      "p32="<<params2[3]<<
+      "p42="<<params2[4]<<
+      "p52="<<params2[5]<<
+      "\n";
+    //    delete his;
+  }
+
+}