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 #include  "TRobustEstimator.h"
  69 #include  "TTimeStamp.h"
  70
  71 #include  "AliLog.h"
  72 #include  "AliMagF.h"
  73
  74 #include  "AliESD.h"
  75 #include  "AliESDfriend.h"
  76 #include  "AliESDtrack.h"
  77 #include  "AliTracker.h"
  78 #include  "AliTPCseed.h"
  79 #include  "AliTPCclusterMI.h"
  80 #include  "AliTPCParamSR.h"
  81 #include  "AliTPCROC.h"
  82
  83
  84 #include  "TTreeStream.h"
  85 #include  "TF1.h"
  86 #include  "TGraph.h"
  87 #include  "AliSignalProcesor.h"
  88 #include  "TCanvas.h"
  89
  90
  91 void FitSignals(TTree * treeB, TCut cut="Max-Median>150&&RMS06<2&&abs(Median-Mean09)<0.5", Int_t maxS=100);
  92 void LaserCalib(TTreeSRedirector & cstream, TTree * chain, Float_t tmin, Float_t tmax, Float_t fraction=0.7);
  93
  94 void AnalyzeESDtracks(Int_t run){
  95   //
  96   // output redirect
  97   TTreeSRedirector *  pcstream = new TTreeSRedirector("TPCtracks.root");
  98   TTreeSRedirector &cstream = *pcstream;
  99   //
 100   // dummy magnetic field
 101   AliMagF mag("aaa","aaa",1,1,10);
 102   AliTracker::SetFieldMap(&mag,kTRUE);
 103   TFile f("AliESDs.root");
 104   TTree * tree =(TTree*)f.Get("esdTree");
 105   AliESD * esd =0;
 106   tree->SetBranchAddress("ESD",&esd);
 107   AliESDfriend *evf=0;
 108   tree->AddFriend("esdFriendTree","AliESDfriends.root");
 109   tree->SetBranchAddress("ESDfriend",&evf);
 110
 111   Int_t nevents = tree->GetEntries();
 112   TClonesArray *clusters = new TClonesArray("AliTPCclusterMI",160);
 113   for (Int_t irow=0; irow<160; irow++){
 114     new ((*clusters)[irow])  AliTPCclusterMI;   // iitial dummy clusters
 115   }
 116
 117   for (Int_t ievent=0; ievent<nevents; ievent++){
 118     tree->GetEntry(ievent);
 119     if (!esd) continue;
 120     if (!evf) continue;
 121     esd->SetESDfriend(evf); //Attach the friend to the ESD
 122     for (Int_t itrack =0; itrack<esd->GetNumberOfTracks(); itrack++){
 123       // Int_t itrack = 0;
 124       if (esd->GetTrack(itrack)->GetFriendTrack()==0) continue;
 125       AliESDtrack * etrack = esd->GetTrack(itrack);
 126       AliESDfriendTrack * ftrack = (AliESDfriendTrack *)esd->GetTrack(itrack)->GetFriendTrack();
 127       AliTPCseed * seed =  (AliTPCseed*)(ftrack->GetCalibObject(0));
 128       if (!seed) continue;
 129       for (Int_t irow=0; irow<160; irow++){
 130         if (seed->GetClusterFast(irow)){
 131           AliTPCclusterMI * cl = new ((*clusters)[irow])  AliTPCclusterMI(*(seed->GetClusterFast(irow)));
 132           cl->SetLabel(itrack,0);
 133         }
 134         else{
 135           AliTPCclusterMI * cl = (AliTPCclusterMI*)clusters->At(irow);
 136           cl->SetX(0); cl->SetY(0); cl->SetZ(0); cl->SetQ(0); cl->SetLabel(-1,0);
 137         }
 138       }
 139       Float_t dEdx = seed->GetdEdx();
 140       Float_t dEdxI = seed->CookdEdx(0.05,0.6,0,77);
 141       Float_t dEdxO = seed->CookdEdx(0.05,0.6,78,155);
 142       Int_t ncl = seed->GetNumberOfClusters();
 143       cstream<<"Tracks"<<
 144         "Run="<<run<<
 145         "Ncl="<<ncl<<
 146         "Event="<<ievent<<
 147         "dEdx="<<dEdx<<
 148         "dEdxI="<<dEdxI<<
 149         "dEdxO="<<dEdxO<<
 150         "Track.="<<seed<<
 151         "Etrack.="<<etrack<<
 152         "Cl.="<<clusters<<
 153         "\n";
 154     }
 155   }
 156   delete pcstream;
 157   //
 158   // Fit signal part
 159   //
 160   TFile fs("TPCsignal.root");
 161   TTree *treeB =(TTree*)fs.Get("SignalB");
 162   //
 163   // Fit central electrode part
 164   //
 165   TTreeSRedirector * pcestream = new  TTreeSRedirector("TimeRoot.root");
 166   TTree * treece = (TTree*)fs.Get("Signalce");
 167   if (tree) {
 168     LaserCalib(*pcestream, treece, 800,1000, 0.7);
 169     delete pcestream;
 170   }
 171   FitSignals(treeB,"Max-Median>150&&RMS06<1.0&&RMS09<1.5&&abs(Median-Mean09)<0.2&&abs(Mean06-Mean09)<0.2",1000);
 172   //
 173 }
 174
 175
 176 void FitSignals(TTree * treeB, TCut cut, Int_t max){
 177   AliSignalProcesor proc;
 178   TF1 * f1 = proc.GetAsymGauss();
 179   TTreeSRedirector cstream("FitSignal.root");
 180   TFile *f = cstream.GetFile();
 181
 182   char lname[100];
 183   sprintf(lname,"Fit%s", cut.GetTitle());
 184   TEventList *list = new TEventList(lname,lname);
 185   sprintf(lname,">>Fit%s", cut.GetTitle());
 186   treeB->Draw(lname,cut);
 187   treeB->SetEventList(list);
 188   Int_t nFits=0;
 189   for (Int_t ievent=0; ievent<list->GetN(); ievent++){
 190     if (nFits>max) break;
 191     if (nFits%50==0) printf("%d\n",nFits);
 192     char ename[100];
 193     sprintf(ename,"Fit%d", ievent);
 194     Double_t nsample = treeB->Draw("Graph.fY-Mean09:Graph.fX","","",1,ievent);
 195     Double_t * signal  = treeB->GetV1();
 196     Double_t * time  = treeB->GetV2();
 197     Double_t maxpos =0;
 198     Double_t max = 0;
 199     for (Int_t ipos = 0; ipos<nsample; ipos++){
 200       if (signal[ipos]>max){
 201         max    = signal[ipos];
 202         maxpos = ipos;
 203       }
 204     }
 205
 206     Int_t first = TMath::Max(maxpos-10,0.);
 207     Int_t last  = TMath::Min(maxpos+60, nsample);
 208     //
 209     f->cd();
 210     TH1F his(ename,ename,last-first,first,last);
 211     for (Int_t ipos=0; ipos<last-first; ipos++){
 212       his.SetBinContent(ipos+1,signal[ipos+first]);
 213     }
 214     treeB->Draw("Sector:Row:Pad","","",1,ievent);
 215     Double_t sector = treeB->GetV1()[0];
 216     Double_t row    = treeB->GetV2()[0];
 217     Double_t pad    = treeB->GetV3()[0];
 218     //    TGraph  graph(last-first,&time[first],&signal[first]);
 219     f1->SetParameters(0.75*max,maxpos,1.1,0.8,0.25,0.2);
 220     //    TH1F * his = (TH1F*)graph.GetHistogram();
 221     his.Fit(f1,"q");
 222     his.Write(ename);
 223     gPad->Clear();
 224     his.Draw();
 225     gPad->Update();
 226     Double_t params[6];
 227     for (Int_t ipar=0; ipar<6; ipar++) params[ipar] = f1->GetParameters()[ipar];
 228     Double_t chi2 = TFitter::GetFitter()->Chisquare(6,params);
 229     TMatrixD cov(6,6);
 230     cov.SetMatrixArray(TFitter::GetFitter()->GetCovarianceMatrix());
 231     //
 232     // tail cancellation
 233     //
 234     Double_t x0[1000];
 235     Double_t x1[1000];
 236     Double_t x2[1000];
 237     for (Int_t ipos=0; ipos<last-first; ipos++){
 238       x0[ipos] = signal[ipos+first];
 239     }
 240     proc.TailCancelationALTRO1(x0,x1,0.85*0.339,0.09,last-first);
 241     proc.TailCancelationALTRO1(x1,x2,0.85,0.789,last-first);
 242     //
 243     sprintf(ename,"Cancel1_%d", ievent);
 244     TH1F his1(ename,ename,last-first,first,last);
 245     for (Int_t ipos=0; ipos<last-first; ipos++){
 246       his1.SetBinContent(ipos+1,x1[ipos]);
 247     }
 248     his1.Write(ename);
 249     sprintf(ename,"Cancel2_%d", ievent);
 250     TH1F his2(ename,ename,last-first,first,last);
 251     for (Int_t ipos=0; ipos<last-first; ipos++){
 252       his2.SetBinContent(ipos+1,x1[ipos]);
 253     }
 254     f1->SetParameters(0.75*max,maxpos,1.1,0.8,0.25,0.2);
 255     his2.Fit(f1,"q");
 256     his2.Write(ename);
 257     Double_t params2[6];
 258     for (Int_t ipar=0; ipar<6; ipar++) params2[ipar] = f1->GetParameters()[ipar];
 259     Double_t chi22 = TFitter::GetFitter()->Chisquare(6,params2);
 260     TMatrixD cov2(6,6);
 261     cov2.SetMatrixArray(TFitter::GetFitter()->GetCovarianceMatrix());
 262
 263     TGraph gr0(last-first, &time[first],x0);
 264     TGraph gr1(last-first, &time[first],x1);
 265     TGraph gr2(last-first, &time[first],x2);
 266     //
 267     cstream<<"Fit"<<
 268       "Sector="<<sector<<
 269       "Row="<<row<<
 270       "Pad="<<pad<<
 271       "First="<<first<<
 272       "Max="<<max<<
 273       "MaxPos="<<maxpos<<
 274       "chi2="<<chi2<<
 275       "chi22="<<chi22<<
 276       "Cov="<<&cov<<
 277       "Cov2="<<&cov2<<
 278       "gr0.="<<&gr0<<
 279       "gr1.="<<&gr1<<
 280       "gr2.="<<&gr2<<
 281       "p0="<<params[0]<<
 282       "p1="<<params[1]<<
 283       "p2="<<params[2]<<
 284       "p3="<<params[3]<<
 285       "p4="<<params[4]<<
 286       "p5="<<params[5]<<
 287       "p02="<<params2[0]<<
 288       "p12="<<params2[1]<<
 289       "p22="<<params2[2]<<
 290       "p32="<<params2[3]<<
 291       "p42="<<params2[4]<<
 292       "p52="<<params2[5]<<
 293       "\n";
 294     //    delete his;
 295     nFits++;
 296   }
 297
 298 }
 299
 300
 301 void LaserCalib(TTreeSRedirector & cstream, TTree * chain, Float_t tmin, Float_t tmax, Float_t fraction){
 302   //
 303   //
 304   //
 305   const Double_t kMaxDelta=10;
 306   AliTPCParamSR param;
 307   param.Update();
 308   TFile fce("TPCsignal.root");
 309   TTree   * treece =(TTree*)fce.Get("Signalce");
 310   if (chain) treece=chain;
 311   //
 312   TBranch * brsector  = treece->GetBranch("Sector");
 313   TBranch * brpad     = treece->GetBranch("Pad");
 314   TBranch * brrow     = treece->GetBranch("Row");
 315   TBranch * brTimeStamp = treece->GetBranch("TimeStamp");
 316   //
 317   TBranch * brtime    = treece->GetBranch("Time");
 318   TBranch * brrms     = treece->GetBranch("RMS06");
 319   TBranch * brmax     = treece->GetBranch("Max");
 320   TBranch * brsum     = treece->GetBranch("Qsum");
 321
 322   Int_t sector, pad, row=0;
 323   Double_t time=0, rms=0, qMax=0, qSum=0;
 324   UInt_t  timeStamp=0;
 325   brsector->SetAddress(&sector);
 326   brrow->SetAddress(&row);
 327   brpad->SetAddress(&pad);
 328   brTimeStamp->SetAddress(&timeStamp);
 329
 330   brtime->SetAddress(&time);
 331   brrms->SetAddress(&rms);
 332   brmax->SetAddress(&qMax);
 333   brsum->SetAddress(&qSum);
 334
 335
 336   brsector->GetEntry(0);
 337   //
 338   Int_t firstSector  = sector;
 339   Int_t lastSector   = sector;
 340   Int_t fentry = 0;
 341   Int_t lentry = 0;
 342   //
 343   // find time offset for differnt events
 344   //
 345   Int_t count = 0;
 346   Double_t padTimes[500000];
 347   TRobustEstimator restim;
 348   Double_t meanS[72], sigmaS[72];
 349   Int_t   firstS[72], lastS[72];
 350   Double_t   sectorsID[72];
 351   for (Int_t isector=0;isector<72; isector++){
 352     firstS[isector]=-1;
 353     lastS[isector] =-1;
 354   };
 355   TH1F  hisT("hisT","hisT",100,tmin,tmax);
 356   treece->Draw("Time>>hisT","");
 357   Float_t cbin = hisT.GetBinCenter(hisT.GetMaximumBin());
 358
 359   for (Int_t ientry=0; ientry<treece->GetEntriesFast(); ientry++){
 360     treece->GetEvent(ientry);
 361     //
 362     if (sector!=lastSector && sector==firstSector){
 363       //if (sector!=lastSector){
 364       lentry = ientry;
 365       TTimeStamp stamp(timeStamp);
 366       stamp.Print();
 367       printf("\nEvent\t%d\tFirst\t%d\tLast\t%d\t%d\n",count, fentry, lentry, lentry-fentry);
 368       //
 369       //
 370       Int_t ngood=0;
 371       for (Int_t ientry2=fentry; ientry2<lentry; ientry2++){
 372         //      brtime->GetEvent(ientry2);
 373         //  brsector->GetEvent(ientry2);
 374         treece->GetEvent(ientry2);
 375         if (time>tmin&&time<tmax && TMath::Abs(time-cbin)<kMaxDelta){
 376           padTimes[ngood]=time;
 377           ngood++;
 378           if (firstS[sector]<0)  firstS[sector]= ngood;
 379           if (firstS[sector]>=0) lastS[sector] = ngood;
 380         }
 381       }
 382       //
 383       //
 384       Double_t mean,sigma;
 385       restim.EvaluateUni(ngood,padTimes,mean, sigma,int(float(ngood)*fraction));
 386       printf("Event\t%d\t%f\t%f\n",count, mean, sigma);
 387       for (Int_t isector=0; isector<72; isector++){
 388         sectorsID[isector]=sector;
 389         if (firstS[isector]>=0 &&lastS[isector]>=0 && lastS[isector]>firstS[isector] ){
 390           Int_t ngoodS = lastS[isector]-firstS[isector];
 391           restim.EvaluateUni(ngoodS, &padTimes[firstS[isector]],meanS[isector],
 392                              sigmaS[isector],int(float(ngoodS)*fraction));
 393         }
 394       }
 395       TGraph  graphM(72,sectorsID,meanS);
 396       TGraph  graphS(72,sectorsID,sigmaS);
 397       cstream<<"TimeS"<<
 398         "CBin="<<cbin<<
 399         "Event="<<count<<
 400         "GM="<<&graphM<<
 401         "GS="<<&graphS<<
 402         "\n";
 403
 404
 405       for (Int_t ientry2=fentry; ientry2<lentry-1; ientry2++){
 406         treece->GetEvent(ientry2);
 407         Double_t x      = param.GetPadRowRadii(sector,row);
 408         Int_t    maxpad = AliTPCROC::Instance()->GetNPads(sector,row);
 409         Double_t y = (pad - 2.5 - 0.5*maxpad)*param.GetPadPitchWidth(sector);
 410         Double_t alpha = TMath::DegToRad()*(10.+20.*(sector%18));
 411         Double_t gx = x*TMath::Cos(alpha)-y*TMath::Sin(alpha);
 412         Double_t gy = y*TMath::Cos(alpha)+x*TMath::Sin(alpha);
 413
 414         Int_t npadS = lastS[sector]-firstS[sector];
 415         cstream<<"Time"<<
 416           "Event="<<count<<
 417           "TimeStamp="<<timeStamp<<
 418           "CBin="<<cbin<<
 419           "x="<<x<<
 420           "y="<<y<<
 421           "gx="<<gx<<
 422           "gy="<<gy<<
 423           "Sector="<<sector<<
 424           "Row="<<row<<
 425           "Pad="<<pad<<
 426           "Time="<<time<<
 427           "RMS="<<rms<<
 428           "Time0="<<mean<<
 429           "Sigma0="<<sigma<<
 430           "TimeS0="<<meanS[sector]<<
 431           "SigmaS0="<<sigmaS[sector]<<
 432           "npad0="<<ngood<<
 433           "npadS="<<npadS<<
 434           "Max="<<qMax<<
 435           "Sum="<<qSum<<
 436           "\n";
 437       }
 438       treece->GetEvent(ientry);
 439       fentry = ientry;
 440       count++;
 441       for (Int_t isector=0;isector<72; isector++){
 442         firstS[isector]=-1;
 443         lastS[isector] =-1;
 444       }
 445     }
 446     lastSector=sector;
 447   }
 448 }
 449
 450
 451
 452
 453 TChain *MakeChainCL(Int_t first, Int_t last){
 454   TChain *chaincl = new TChain("TreeR","TreeR");
 455   //
 456   char fname[100];
 457   for (Int_t i=first;i<last; i++){
 458     if (i>0) sprintf(fname,"TPC.RecPoints%d.root/Event%d/TreeR",i,i);
 459     if (i==0) sprintf(fname,"TPC.RecPoints.root/Event%d/TreeR",i);
 460     chaincl->Add(fname);
 461   }
 462   return chaincl;
 463 }
 464
 465 TTree* GetTree(Int_t ievent){
 466   char fname[100];
 467   char tname[100];
 468   if (ievent>0) sprintf(fname,"TPC.RecPoints%d.root",ievent);
 469   if (ievent==0) sprintf(fname,"TPC.RecPoints.root");
 470   sprintf(tname,"Event%d/TreeR",ievent);
 471   TFile * f  = new TFile(fname);
 472   TTree * tree = (TTree*)f->Get(tname);
 473   return tree;
 474
 475 }