TPC/AnalyzeESDtracks.C

   1 //
   2 // Process ESD tracks  -
   3 // Extract TPC tracks  - write them to tree
   4 //
   5 /*
   6   .L AnalyzeESDtracks.C+
   7   .L AliGenInfo.C+
   8   AnalyzeESDtracks(567);   // process tracks
   9   // Tracks are written to the file "TPCtracks.root"
  10   // Now yo can analyze it
  11   TFile fesd("AliESDs.root");
  12   TTree * treeE = (TTree*)fesd.Get("esdTree");
  13   TFile f("TPCtracks.root")
  14   TTree * tree =(TTree*)f.Get("Tracks");
  15   AliComparisonDraw comp;
  16   comp->fTree = tree;
  17
  18 */
  19
  20 //
  21
  22 /*
  23 .L AnalyzeESDtracks.C+
  24 .L AliGenInfo.C+
  25
  26 //AnalyzeESDtracks(567);
  27
  28
  29 TFile fesd("AliESDs.root");
  30 TTree * treeE = (TTree*)fesd.Get("esdTree");
  31 TFile f("TPCtracks.root")
  32 TTree * tree =(TTree*)f.Get("Tracks");
  33 AliComparisonDraw comp;
  34 comp->fTree = tree;
  35
  36 TFile fs("TPCsignal.root");
  37 TTree *treeB =(TTree*)fs.Get("SignalB");
  38 TTree *treeN =(TTree*)fs.Get("SignalN");
  39 TTree *treeS =(TTree*)fs.Get("Signal");
  40 TTree *treef =(TTree*)fs.Get("Fit");
  41
  42
  43 FitSignals(treeB,"Max-Median>100&&RMS06<2.5")
  44 TFile ffit("FitSignal.root");
  45 TTree * treeF = (TTree*)ffit->Get("Fit");
  46
  47
  48 TChain chaincl("TreeR","TreeR")
  49 chaincl.Add("TPC.RecPoints.root/Event0/TreeR")
  50 chaincl.Add("TPC.RecPoints1.root/Event1/TreeR")
  51 chaincl.Add("TPC.RecPoints2.root/Event2/TreeR")
  52 chaincl.Add("TPC.RecPoints3.root/Event3/TreeR")
  53
  54 */
  55
  56
  57
  58 #include "TObject.h"
  59 #include "TFile.h"
  60 #include "TTree.h"
  61 #include "TChain.h"
  62 #include "TBranch.h"
  63 #include "TTreeStream.h"
  64 #include "TEventList.h"
  65 #include "TCut.h"
  66 #include "TFitter.h"
  67 #include  "TMatrixD.h"
  68
  69 #include  "AliLog.h"
  70 #include  "AliMagF.h"
  71
  72 #include  "AliESD.h"
  73 #include  "AliESDfriend.h"
  74 #include  "AliESDtrack.h"
  75 #include  "AliTracker.h"
  76 #include  "AliTPCseed.h"
  77 #include  "AliTPCclusterMI.h"
  78 #include  "TTreeStream.h"
  79 #include  "TF1.h"
  80 #include  "TGraph.h"
  81 #include  "AliSignalProcesor.h"
  82 #include  "TCanvas.h"
  83
  84
  85 void FitSignals(TTree * treeB, TCut cut="Max-Median>150&&RMS06<2&&abs(Median-Mean09)<0.5");
  86
  87
  88 void AnalyzeESDtracks(Int_t run){
  89   //
  90   // output redirect
  91   TTreeSRedirector *  pcstream = new TTreeSRedirector("TPCtracks.root");
  92   TTreeSRedirector &cstream = *pcstream;
  93   //
  94   // dummy magnetic field
  95   AliMagF mag("aaa","aaa",1,1,10);
  96   AliTracker::SetFieldMap(&mag,kTRUE);
  97   TFile f("AliESDs.root");
  98   TTree * tree =(TTree*)f.Get("esdTree");
  99   AliESD * esd =0;
 100   tree->SetBranchAddress("ESD",&esd);
 101   AliESDfriend *evf=0;
 102   tree->AddFriend("esdFriendTree","AliESDfriends.root");
 103   tree->SetBranchAddress("ESDfriend",&evf);
 104
 105   Int_t nevents = tree->GetEntries();
 106   TClonesArray *clusters = new TClonesArray("AliTPCclusterMI",160);
 107   for (Int_t irow=0; irow<160; irow++){
 108     new ((*clusters)[irow])  AliTPCclusterMI;   // iitial dummy clusters
 109   }
 110
 111   for (Int_t ievent=0; ievent<nevents; ievent++){
 112     tree->GetEntry(ievent);
 113     if (!esd) continue;
 114     if (!evf) continue;
 115     esd->SetESDfriend(evf); //Attach the friend to the ESD
 116     for (Int_t itrack =0; itrack<esd->GetNumberOfTracks(); itrack++){
 117       // Int_t itrack = 0;
 118       if (esd->GetTrack(itrack)->GetFriendTrack()==0) continue;
 119       AliESDtrack * etrack = esd->GetTrack(itrack);
 120       AliESDfriendTrack * ftrack = (AliESDfriendTrack *)esd->GetTrack(itrack)->GetFriendTrack();
 121       AliTPCseed * seed =  (AliTPCseed*)(ftrack->GetCalibObject(0));
 122       if (!seed) continue;
 123       for (Int_t irow=0; irow<160; irow++){
 124         if (seed->GetClusterFast(irow)){
 125           AliTPCclusterMI * cl = new ((*clusters)[irow])  AliTPCclusterMI(*(seed->GetClusterFast(irow)));
 126           cl->SetLabel(itrack,0);
 127         }
 128         else{
 129           AliTPCclusterMI * cl = (AliTPCclusterMI*)clusters->At(irow);
 130           cl->SetX(0); cl->SetY(0); cl->SetZ(0); cl->SetQ(0); cl->SetLabel(-1,0);
 131         }
 132       }
 133       Float_t dEdx = seed->GetdEdx();
 134       Float_t dEdxI = seed->CookdEdx(0.05,0.6,0,77);
 135       Float_t dEdxO = seed->CookdEdx(0.05,0.6,78,155);
 136       Int_t ncl = seed->GetNumberOfClusters();
 137       cstream<<"Tracks"<<
 138         "Run="<<run<<
 139         "Ncl="<<ncl<<
 140         "Event="<<ievent<<
 141         "dEdx="<<dEdx<<
 142         "dEdxI="<<dEdxI<<
 143         "dEdxO="<<dEdxO<<
 144         "Track.="<<seed<<
 145         "Etrack.="<<etrack<<
 146         "Cl.="<<clusters<<
 147         "\n";
 148     }
 149   }
 150   delete pcstream;
 151   //
 152   // Fit signal part
 153   //
 154   TFile fs("TPCsignal.root");
 155   TTree *treeB =(TTree*)fs.Get("SignalB");
 156   //  FitSignals(treeB,"Max-Median>150&&RMS06<2.5");
 157 }
 158
 159
 160 void FitSignals(TTree * treeB, TCut cut){
 161   AliSignalProcesor proc;
 162   TF1 * f1 = proc.GetAsymGauss();
 163   TTreeSRedirector cstream("FitSignal.root");
 164   TFile *f = cstream.GetFile();
 165
 166   char lname[100];
 167   sprintf(lname,"Fit%s", cut.GetTitle());
 168   TEventList *list = new TEventList(lname,lname);
 169   sprintf(lname,">>Fit%s", cut.GetTitle());
 170   treeB->Draw(lname,cut);
 171   treeB->SetEventList(list);
 172   for (Int_t ievent=0; ievent<list->GetN(); ievent++){
 173     char ename[100];
 174     sprintf(ename,"Fit%d", ievent);
 175     Double_t nsample = treeB->Draw("Graph.fY-Mean09:Graph.fX","","",1,ievent);
 176     Double_t * signal  = treeB->GetV1();
 177     Double_t * time  = treeB->GetV2();
 178     Double_t maxpos =0;
 179     Double_t max = 0;
 180     for (Int_t ipos = 0; ipos<nsample; ipos++){
 181       if (signal[ipos]>max){
 182         max    = signal[ipos];
 183         maxpos = ipos;
 184       }
 185     }
 186
 187     Int_t first = TMath::Max(maxpos-10,0.);
 188     Int_t last  = TMath::Min(maxpos+60, nsample);
 189     //
 190     f->cd();
 191     TH1F his(ename,ename,last-first,first,last);
 192     for (Int_t ipos=0; ipos<last-first; ipos++){
 193       his.SetBinContent(ipos+1,signal[ipos+first]);
 194     }
 195     treeB->Draw("Sector:Row:Pad","","",1,ievent);
 196     Double_t sector = treeB->GetV1()[0];
 197     Double_t row    = treeB->GetV2()[0];
 198     Double_t pad    = treeB->GetV3()[0];
 199     //    TGraph  graph(last-first,&time[first],&signal[first]);
 200     f1->SetParameters(0.75*max,maxpos,1.1,0.8,0.25,0.2);
 201     //    TH1F * his = (TH1F*)graph.GetHistogram();
 202     his.Fit(f1,"q");
 203     his.Write(ename);
 204     gPad->Clear();
 205     his.Draw();
 206     gPad->Update();
 207     Double_t params[6];
 208     for (Int_t ipar=0; ipar<6; ipar++) params[ipar] = f1->GetParameters()[ipar];
 209     Double_t chi2 = TFitter::GetFitter()->Chisquare(6,params);
 210     TMatrixD cov(6,6);
 211     cov.SetMatrixArray(TFitter::GetFitter()->GetCovarianceMatrix());
 212     //
 213     // tail cancellation
 214     //
 215     Double_t x0[1000];
 216     Double_t x1[1000];
 217     Double_t x2[1000];
 218     for (Int_t ipos=0; ipos<last-first; ipos++){
 219       x0[ipos] = signal[ipos+first];
 220     }
 221     proc.TailCancelationALTRO1(x0,x1,0.85*0.339,0.09,last-first);
 222     proc.TailCancelationALTRO1(x1,x2,0.85,0.789,last-first);
 223     //
 224     sprintf(ename,"Cancel1_%d", ievent);
 225     TH1F his1(ename,ename,last-first,first,last);
 226     for (Int_t ipos=0; ipos<last-first; ipos++){
 227       his1.SetBinContent(ipos+1,x1[ipos]);
 228     }
 229     his1.Write(ename);
 230     sprintf(ename,"Cancel2_%d", ievent);
 231     TH1F his2(ename,ename,last-first,first,last);
 232     for (Int_t ipos=0; ipos<last-first; ipos++){
 233       his2.SetBinContent(ipos+1,x1[ipos]);
 234     }
 235     f1->SetParameters(0.75*max,maxpos,1.1,0.8,0.25,0.2);
 236     his2.Fit(f1,"q");
 237     his2.Write(ename);
 238     Double_t params2[6];
 239     for (Int_t ipar=0; ipar<6; ipar++) params2[ipar] = f1->GetParameters()[ipar];
 240     Double_t chi22 = TFitter::GetFitter()->Chisquare(6,params2);
 241     TMatrixD cov2(6,6);
 242     cov2.SetMatrixArray(TFitter::GetFitter()->GetCovarianceMatrix());
 243
 244     TGraph gr0(last-first, &time[first],x0);
 245     TGraph gr1(last-first, &time[first],x1);
 246     TGraph gr2(last-first, &time[first],x2);
 247     //
 248     cstream<<"Fit"<<
 249       "Sector="<<sector<<
 250       "Row="<<row<<
 251       "Pad="<<pad<<
 252       "First="<<first<<
 253       "Max="<<max<<
 254       "MaxPos="<<maxpos<<
 255       "chi2="<<chi2<<
 256       "chi22="<<chi22<<
 257       "Cov="<<&cov<<
 258       "Cov2="<<&cov2<<
 259       "gr0.="<<&gr0<<
 260       "gr1.="<<&gr1<<
 261       "gr2.="<<&gr2<<
 262       "p0="<<params[0]<<
 263       "p1="<<params[1]<<
 264       "p2="<<params[2]<<
 265       "p3="<<params[3]<<
 266       "p4="<<params[4]<<
 267       "p5="<<params[5]<<
 268       "p02="<<params2[0]<<
 269       "p12="<<params2[1]<<
 270       "p22="<<params2[2]<<
 271       "p32="<<params2[3]<<
 272       "p42="<<params2[4]<<
 273       "p52="<<params2[5]<<
 274       "\n";
 275     //    delete his;
 276   }
 277
 278 }
 279
 280
 281
 282 TChain *MakeChainCL(Int_t first, Int_t last){
 283   TChain *chaincl = new TChain("TreeR","TreeR");
 284   //
 285   char fname[100];
 286   for (Int_t i=first;i<last; i++){
 287     if (i>0) sprintf(fname,"TPC.RecPoints%d.root/Event%d/TreeR",i,i);
 288     if (i==0) sprintf(fname,"TPC.RecPoints.root/Event%d/TreeR",i);
 289     chaincl->Add(fname);
 290   }
 291   return chaincl;
 292 }
 293
 294 TTree* GetTree(Int_t ievent){
 295   char fname[100];
 296   char tname[100];
 297   if (ievent>0) sprintf(fname,"TPC.RecPoints%d.root",ievent);
 298   if (ievent==0) sprintf(fname,"TPC.RecPoints.root");
 299   sprintf(tname,"Event%d/TreeR",ievent);
 300   TFile * f  = new TFile(fname);
 301   TTree * tree = (TTree*)f->Get(tname);
 302   return tree;
 303
 304 }