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   AnalyzESDtracks(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 //AnalyzESDtracks(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 AnalyzeESDtracks(Int_t run){
  86   //
  87   // output redirect
  88   TTreeSRedirector  cstream("TPCtracks.root");
  89   //
  90   // dummy magnetic field
  91   AliMagF mag("aaa","aaa",1,1,10);
  92   AliTracker::SetFieldMap(&mag,kTRUE);
  93   TFile f("AliESDs.root");
  94   TTree * tree =(TTree*)f.Get("esdTree");
  95   AliESD * esd =0;
  96   tree->SetBranchAddress("ESD",&esd);
  97   AliESDfriend *evf=0;
  98   tree->AddFriend("esdFriendTree","AliESDfriends.root");
  99   tree->SetBranchAddress("ESDfriend",&evf);
 100
 101   Int_t nevents = tree->GetEntries();
 102   TClonesArray *clusters = new TClonesArray("AliTPCclusterMI",160);
 103   for (Int_t irow=0; irow<160; irow++){
 104     new ((*clusters)[irow])  AliTPCclusterMI;   // iitial dummy clusters
 105   }
 106
 107   for (Int_t ievent=0; ievent<nevents; ievent++){
 108     tree->GetEntry(ievent);
 109     if (!esd) continue;
 110     if (!evf) continue;
 111     esd->SetESDfriend(evf); //Attach the friend to the ESD
 112     for (Int_t itrack =0; itrack<esd->GetNumberOfTracks(); itrack++){
 113       // Int_t itrack = 0;
 114       if (esd->GetTrack(itrack)->GetFriendTrack()==0) continue;
 115       AliESDtrack * etrack = esd->GetTrack(itrack);
 116       AliESDfriendTrack * ftrack = (AliESDfriendTrack *)esd->GetTrack(itrack)->GetFriendTrack();
 117       AliTPCseed * seed =  (AliTPCseed*)(ftrack->GetCalibObject(0));
 118       if (!seed) continue;
 119       for (Int_t irow=0; irow<160; irow++){
 120         if (seed->GetClusterFast(irow)){
 121           AliTPCclusterMI * cl = new ((*clusters)[irow])  AliTPCclusterMI(*(seed->GetClusterFast(irow)));
 122           cl->SetLabel(itrack,0);
 123         }
 124         else{
 125           AliTPCclusterMI * cl = (AliTPCclusterMI*)clusters->At(irow);
 126           cl->SetX(0); cl->SetY(0); cl->SetZ(0); cl->SetQ(0); cl->SetLabel(-1,0);
 127         }
 128       }
 129       Float_t dEdx = seed->GetdEdx();
 130       Float_t dEdxI = seed->CookdEdx(0.05,0.6,0,77);
 131       Float_t dEdxO = seed->CookdEdx(0.05,0.6,78,155);
 132       Int_t ncl = seed->GetNumberOfClusters();
 133       cstream<<"Tracks"<<
 134         "Run="<<run<<
 135         "Ncl="<<ncl<<
 136         "Event="<<ievent<<
 137         "dEdx="<<dEdx<<
 138         "dEdxI="<<dEdxI<<
 139         "dEdxO="<<dEdxO<<
 140         "Track.="<<seed<<
 141         "Etrack.="<<etrack<<
 142         "Cl.="<<clusters<<
 143         "\n";
 144     }
 145   }
 146   //
 147   // Fit signal part
 148   //
 149   TFile fs("TPCsignal.root");
 150   TTree *treeB =(TTree*)fs.Get("SignalB");
 151   FitSignals(treeB,"Max-Median>150&&RMS06<2.5");
 152 }
 153
 154
 155 void FitSignals(TTree * treeB, TCut cut="Max-Median>150&&RMS06<2&&abs(Median-Mean09)<0.5"){
 156   AliSignalProcesor proc;
 157   TF1 * f1 = proc.GetAsymGauss();
 158   TTreeSRedirector cstream("FitSignal.root");
 159   TFile *f = cstream.GetFile();
 160
 161   char lname[100];
 162   sprintf(lname,"Fit%s", cut.GetTitle());
 163   TEventList *list = new TEventList(lname,lname);
 164   sprintf(lname,">>Fit%s", cut.GetTitle());
 165   treeB->Draw(lname,cut);
 166   treeB->SetEventList(list);
 167   for (Int_t ievent=0; ievent<list->GetN(); ievent++){
 168     char ename[100];
 169     sprintf(ename,"Fit%d", ievent);
 170     Double_t nsample = treeB->Draw("Graph.fY-Mean09:Graph.fX","","",1,ievent);
 171     Double_t * signal  = treeB->GetV1();
 172     Double_t * time  = treeB->GetV2();
 173     Double_t maxpos =0;
 174     Double_t max = 0;
 175     for (Int_t ipos = 0; ipos<nsample; ipos++){
 176       if (signal[ipos]>max){
 177         max    = signal[ipos];
 178         maxpos = ipos;
 179       }
 180     }
 181
 182     Int_t first = TMath::Max(maxpos-10,0.);
 183     Int_t last  = TMath::Min(maxpos+60, nsample);
 184     //
 185     f->cd();
 186     TH1F his(ename,ename,last-first,first,last);
 187     for (Int_t ipos=0; ipos<last-first; ipos++){
 188       his.SetBinContent(ipos+1,signal[ipos+first]);
 189     }
 190     treeB->Draw("Sector:Row:Pad","","",1,ievent);
 191     Double_t sector = treeB->GetV1()[0];
 192     Double_t row    = treeB->GetV2()[0];
 193     Double_t pad    = treeB->GetV3()[0];
 194     //    TGraph  graph(last-first,&time[first],&signal[first]);
 195     f1->SetParameters(0.75*max,maxpos,1.1,0.8,0.25,0.2);
 196     //    TH1F * his = (TH1F*)graph.GetHistogram();
 197     his.Fit(f1,"q");
 198     his.Write(ename);
 199     gPad->Clear();
 200     his.Draw();
 201     gPad->Update();
 202     Double_t params[6];
 203     for (Int_t ipar=0; ipar<6; ipar++) params[ipar] = f1->GetParameters()[ipar];
 204     Double_t chi2 = TFitter::GetFitter()->Chisquare(6,params);
 205     TMatrixD cov(6,6);
 206     cov.SetMatrixArray(TFitter::GetFitter()->GetCovarianceMatrix());
 207     //
 208     // tail cancellation
 209     //
 210     Double_t x0[1000];
 211     Double_t x1[1000];
 212     Double_t x2[1000];
 213     for (Int_t ipos=0; ipos<last-first; ipos++){
 214       x0[ipos] = signal[ipos+first];
 215     }
 216     proc.TailCancelationALTRO1(x0,x1,0.85*0.339,0.09,last-first);
 217     proc.TailCancelationALTRO1(x1,x2,0.85,0.789,last-first);
 218     //
 219     sprintf(ename,"Cancel1_%d", ievent);
 220     TH1F his1(ename,ename,last-first,first,last);
 221     for (Int_t ipos=0; ipos<last-first; ipos++){
 222       his1.SetBinContent(ipos+1,x1[ipos]);
 223     }
 224     his1.Write(ename);
 225     sprintf(ename,"Cancel2_%d", ievent);
 226     TH1F his2(ename,ename,last-first,first,last);
 227     for (Int_t ipos=0; ipos<last-first; ipos++){
 228       his2.SetBinContent(ipos+1,x1[ipos]);
 229     }
 230     f1->SetParameters(0.75*max,maxpos,1.1,0.8,0.25,0.2);
 231     his2.Fit(f1,"q");
 232     his2.Write(ename);
 233     Double_t params2[6];
 234     for (Int_t ipar=0; ipar<6; ipar++) params2[ipar] = f1->GetParameters()[ipar];
 235     Double_t chi22 = TFitter::GetFitter()->Chisquare(6,params2);
 236     TMatrixD cov2(6,6);
 237     cov2.SetMatrixArray(TFitter::GetFitter()->GetCovarianceMatrix());
 238
 239     TGraph gr0(last-first, &time[first],x0);
 240     TGraph gr1(last-first, &time[first],x1);
 241     TGraph gr2(last-first, &time[first],x2);
 242     //
 243     cstream<<"Fit"<<
 244       "Sector="<<sector<<
 245       "Row="<<row<<
 246       "Pad="<<pad<<
 247       "First="<<first<<
 248       "Max="<<max<<
 249       "MaxPos="<<maxpos<<
 250       "chi2="<<chi2<<
 251       "chi22="<<chi22<<
 252       "Cov="<<&cov<<
 253       "Cov2="<<&cov2<<
 254       "gr0.="<<&gr0<<
 255       "gr1.="<<&gr1<<
 256       "gr2.="<<&gr2<<
 257       "p0="<<params[0]<<
 258       "p1="<<params[1]<<
 259       "p2="<<params[2]<<
 260       "p3="<<params[3]<<
 261       "p4="<<params[4]<<
 262       "p5="<<params[5]<<
 263       "p02="<<params2[0]<<
 264       "p12="<<params2[1]<<
 265       "p22="<<params2[2]<<
 266       "p32="<<params2[3]<<
 267       "p42="<<params2[4]<<
 268       "p52="<<params2[5]<<
 269       "\n";
 270     //    delete his;
 271   }
 272
 273 }
 274
 275
 276
 277 TChain *MakeChainCL(Int_t first, Int_t last){
 278   TChain *chaincl = new TChain("TreeR","TreeR");
 279   //
 280   char fname[100];
 281   for (Int_t i=first;i<last; i++){
 282     if (i>0) sprintf(fname,"TPC.RecPoints%d.root/Event%d/TreeR",i,i);
 283     if (i==0) sprintf(fname,"TPC.RecPoints.root/Event%d/TreeR",i);
 284     chaincl->Add(fname);
 285   }
 286   return chaincl;
 287 }
 288
 289 TTree* GetTree(Int_t ievent){
 290   char fname[100];
 291   char tname[100];
 292   if (ievent>0) sprintf(fname,"TPC.RecPoints%d.root",ievent);
 293   if (ievent==0) sprintf(fname,"TPC.RecPoints.root");
 294   sprintf(tname,"Event%d/TreeR",ievent);
 295   TFile * f  = new TFile(fname);
 296   TTree * tree = (TTree*)f->Get(tname);
 297   return tree;
 298
 299 }