TRD/AliTRDCalibTask.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 /////////////////////////////////////////////////////////////////////////////////
  18 //
  19 // AliTRDCalibTask
  20 //
  21 // Offline TRD calibration task
  22 //
  23 // Author:
  24 //   R. Bailhache (R.Bailhache@gsi.de)
  25 //
  26 //////////////////////////////////////////////////////////////////////////////////////
  27
  28
  29
  30 #include "Riostream.h"
  31 #include "TChain.h"
  32 #include "TTree.h"
  33 #include "TFile.h"
  34 #include "TProfile2D.h"
  35 #include "TH2I.h"
  36 #include "TH1F.h"
  37 #include "TH2F.h"
  38 #include "TList.h"
  39 #include "TMath.h"
  40 #include "TObject.h"
  41 #include "TObjArray.h"
  42 #include "TString.h"
  43 #include "TCanvas.h"
  44 #include "TLegend.h"
  45 #include "TStyle.h"
  46 #include "TLine.h"
  47
  48 #include "AliAnalysisTask.h"
  49 #include "AliAnalysisManager.h"
  50
  51 #include "AliESDVertex.h"
  52 #include "AliESDEvent.h"
  53 #include "AliESDfriend.h"
  54 #include "AliESDInputHandler.h"
  55 #include "AliESDtrack.h"
  56 #include "AliESDfriendTrack.h"
  57 #include "AliTRDtrackV1.h"
  58 #include "AliTRDseedV1.h"
  59 #include "AliTRDcluster.h"
  60 #include "AliTRDgeometry.h"
  61 #include "AliESDtrackCuts.h"
  62 #include "AliESDVertex.h"
  63 #include "AliTRDCalDet.h"
  64
  65 #include "AliTRDCalibraFillHisto.h"
  66 #include "AliTRDCalibraVdriftLinearFit.h"
  67
  68 #include "AliTRDcalibDB.h"
  69
  70
  71 #include "AliTRDCalibTask.h"
  72
  73
  74 ClassImp(AliTRDCalibTask)
  75
  76 //________________________________________________________________________
  77   AliTRDCalibTask::AliTRDCalibTask(const char *name)
  78     : AliAnalysisTask(name,""), fESD(0),
  79       fESDfriend(0),
  80       fkEsdTrack(0),
  81       fFriendTrack(0),
  82       fCalibObject(0),
  83       fTrdTrack(0),
  84       fCl(0),
  85       fListHist(0),
  86       fTRDCalibraFillHisto(0),
  87       fNEvents(0),
  88       fNbTRDTrack(0),
  89       fNbTRDTrackOffline(0),
  90       fNbTRDTrackStandalone(0),
  91       fNbTPCTRDtrack(0),
  92       fNbTimeBin(0),
  93       fNbTimeBinOffline(0),
  94       fNbTimeBinStandalone(0),
  95       fNbClusters(0),
  96       fNbClustersOffline(0),
  97       fNbClustersStandalone(0),
  98       fNbTracklets(0),
  99       fNbTrackletsOffline(0),
 100       fNbTrackletsStandalone(0),
 101       fAbsoluteGain(0),
 102       fCH2dSum(0),
 103       fPH2dSum(0),
 104       fCH2dSM(0),
 105       fPH2dSM(0),
 106       fHisto2d(kTRUE),
 107       fVector2d(kFALSE),
 108       fVdriftLinear(kTRUE),
 109       fNbTimeBins(0),
 110       fSelectedTrigger(new TObjArray()),
 111       fRejected(kTRUE),
 112       fEsdTrackCuts(0),
 113       fRequirePrimaryVertex(kFALSE),
 114       fVtxTPC(kFALSE),
 115       fVtxSPD(kFALSE),
 116       fMinNbContributors(0),
 117       fRangePrimaryVertexZ(9999999.0),
 118       fLow(0),
 119       fHigh(30),
 120       fFillZero(kFALSE),
 121       fNormalizeNbOfCluster(kFALSE),
 122       fRelativeScale(0.0),
 123       fMaxCluster(100.0),
 124       fNbMaxCluster(2),
 125       fOfflineTracks(kFALSE),
 126       fStandaloneTracks(kFALSE),
 127       fCalDetGain(0x0),
 128       fMaxEvent(0),
 129       fCounter(0),
 130       fDebug(0)
 131 {
 132   //
 133   // Default constructor
 134   //
 135
 136   fNz[0] = 0;
 137   fNz[1] = 0;
 138   fNz[2] = 0;
 139
 140   fNrphi[0] = 0;
 141   fNrphi[1] = 0;
 142   fNrphi[2] = 0;
 143
 144   // Define input and output slots here
 145   // Input slot #0 works with a TChain
 146   DefineInput(0, TChain::Class());
 147
 148   // Output slot #0 writes into a TList container
 149   DefineOutput(0, TList::Class());
 150
 151
 152 }
 153 //____________________________________________________________________________________
 154 AliTRDCalibTask::~AliTRDCalibTask()
 155 {
 156   //
 157   // AliTRDCalibTask destructor
 158   //
 159
 160   // Pointeur
 161   if(fNEvents) delete fNEvents;
 162   if(fNbTRDTrack) delete fNbTRDTrack;
 163   if(fNbTRDTrackOffline) delete fNbTRDTrackOffline;
 164   if(fNbTRDTrackStandalone) delete fNbTRDTrackStandalone;
 165   if(fNbTPCTRDtrack) delete fNbTPCTRDtrack;
 166   if(fNbTimeBin) delete fNbTimeBin;
 167   if(fNbTimeBinOffline) delete fNbTimeBinOffline;
 168   if(fNbTimeBinStandalone) delete fNbTimeBinStandalone;
 169   if(fNbClusters) delete fNbClusters;
 170   if(fNbClustersOffline) delete fNbClustersOffline;
 171   if(fNbClustersStandalone) delete fNbClustersStandalone;
 172   if(fNbTracklets) delete fNbTracklets;
 173   if(fNbTrackletsOffline) delete fNbTrackletsOffline;
 174   if(fNbTrackletsStandalone) delete fNbTrackletsStandalone;
 175   if(fAbsoluteGain) delete fAbsoluteGain;
 176   if(fCH2dSum) delete fCH2dSum;
 177   if(fPH2dSum) delete fPH2dSum;
 178   if(fCH2dSM) delete fCH2dSM;
 179   if(fPH2dSM) delete fPH2dSM;
 180   if(fCalDetGain) delete fCalDetGain;
 181
 182   if(fSelectedTrigger) {
 183     fSelectedTrigger->Delete();
 184     delete fSelectedTrigger;
 185   }
 186   if(fEsdTrackCuts) {
 187     delete fEsdTrackCuts;
 188   }
 189
 190 }
 191 //________________________________________________________________________
 192 void AliTRDCalibTask::ConnectInputData(Option_t *)
 193 {
 194   // Connect ESD or AOD here
 195   // Called once
 196
 197   TTree* tree = dynamic_cast<TTree*> (GetInputData(0)); //pointer wird "umgecastet" auf anderen Variablentyp
 198   if (!tree) {
 199     //Printf("ERROR: Could not read chain from input slot 0");
 200   } else {
 201
 202     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 203
 204     if (!esdH) {
 205       //Printf("ERROR: Could not get ESDInputHandler");
 206     } else {
 207       fESD = esdH->GetEvent();
 208       esdH->SetActiveBranches("ESDfriend*");
 209       //Printf("*** CONNECTED NEW EVENT ****");
 210     }
 211
 212   }
 213 }
 214 //________________________________________________________________________
 215 void AliTRDCalibTask::CreateOutputObjects()
 216 {
 217   //
 218   // Create the histos
 219   //
 220
 221   // Number of time bins
 222   if(fNbTimeBins==0) {
 223     AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
 224     fNbTimeBins = cal->GetNumberOfTimeBinsDCS();
 225     if(fNbTimeBins <= 0){
 226       AliWarning(Form("No of TimeBins from DB [%d] use default [30]", fNbTimeBins));
 227       fNbTimeBins = 30;
 228     }
 229   }
 230
 231   // instance calibration
 232   fTRDCalibraFillHisto = AliTRDCalibraFillHisto::Instance();
 233   fTRDCalibraFillHisto->SetHisto2d(fHisto2d); // choose to use histograms
 234   fTRDCalibraFillHisto->SetVector2d(fVector2d); // choose to use vectors
 235   fTRDCalibraFillHisto->SetCH2dOn();  // choose to calibrate the gain
 236   fTRDCalibraFillHisto->SetPH2dOn();  // choose to calibrate the drift velocity
 237   fTRDCalibraFillHisto->SetPRF2dOn(); // choose to look at the PRF
 238   fTRDCalibraFillHisto->SetLinearFitterOn(fVdriftLinear); // Other possibility vdrift VDRIFT
 239   fTRDCalibraFillHisto->SetLinearFitterDebugOn(fVdriftLinear); // Other possibility vdrift
 240   for(Int_t k = 0; k < 3; k++){
 241     if(((fNz[k] != 10) && (fNrphi[k] != 10)) && ((fNz[k] != 100) && (fNrphi[k] != 100))) {
 242       fTRDCalibraFillHisto->SetNz(k,fNz[k]);                                    // Mode calibration
 243       fTRDCalibraFillHisto->SetNrphi(k,fNrphi[k]);                             // Mode calibration
 244     }
 245     else {
 246       if((fNz[k] == 100) && (fNrphi[k] == 100))  {
 247         if(fVector2d) AliInfo("The mode all together is not supported by the vector method");
 248         fTRDCalibraFillHisto->SetAllTogether(k);
 249       }
 250       if((fNz[k] == 10) && (fNrphi[k] == 10))  {
 251         if(fVector2d) AliInfo("The mode per supermodule is not supported by the vector method");
 252         fTRDCalibraFillHisto->SetPerSuperModule(k);
 253       }
 254     }
 255   }
 256   // Variables for how to fill
 257   fTRDCalibraFillHisto->SetFillWithZero(fFillZero);
 258   fTRDCalibraFillHisto->SetNormalizeNbOfCluster(fNormalizeNbOfCluster);
 259   fTRDCalibraFillHisto->SetMaxCluster(fMaxCluster);
 260   fTRDCalibraFillHisto->SetNbMaxCluster(fNbMaxCluster);
 261
 262   // Init with 30 timebins
 263   fTRDCalibraFillHisto->Init2Dhistos(fNbTimeBins); // initialise the histos
 264   fTRDCalibraFillHisto->SetNumberClusters(fLow); // At least 11 clusters
 265   fTRDCalibraFillHisto->SetNumberClustersf(fHigh); // At least 11 clusters
 266   fRelativeScale = fTRDCalibraFillHisto->GetRelativeScale(); // Get the relative scale for the gain
 267
 268   // For testing only
 269   if(fDebug > 2) fTRDCalibraFillHisto->SetDebugLevel(1); //debug stuff
 270
 271   // output list
 272   fListHist = new TList();
 273   if(fHisto2d) {
 274     fListHist->Add(fTRDCalibraFillHisto->GetCH2d());
 275     fListHist->Add(fTRDCalibraFillHisto->GetPH2d());
 276     fListHist->Add(fTRDCalibraFillHisto->GetPRF2d());
 277   }
 278   if(fVdriftLinear) fListHist->Add((TObject *) fTRDCalibraFillHisto->GetVdriftLinearFit());
 279   if(fVector2d) fListHist->Add((TObject *) fTRDCalibraFillHisto->GetCalibraVector()); //calibra vector
 280   fNEvents = new TH1I("NEvents","NEvents", 2, 0, 2);
 281   fListHist->Add(fNEvents);
 282
 283   // absolute gain calibration even without AliESDfriend
 284   Int_t nBinsPt = 25;
 285   Double_t minPt = 0.001;
 286   Double_t maxPt = 10.0;
 287
 288   Double_t *binLimLogPt = new Double_t[nBinsPt+1];
 289   Double_t *binLimPt    = new Double_t[nBinsPt+1];
 290   for(Int_t i=0; i<=nBinsPt; i++) binLimLogPt[i]=(Double_t)TMath::Log10(minPt) + (TMath::Log10(maxPt)-TMath::Log10(minPt))/nBinsPt*(Double_t)i ;
 291   for(Int_t i=0; i<=nBinsPt; i++) binLimPt[i]=(Double_t)TMath::Power(10,binLimLogPt[i]);
 292
 293   fAbsoluteGain = new TH2F("AbsoluteGain","AbsoluteGain", 200, 0.0, 700.0, nBinsPt, binLimPt);
 294   fAbsoluteGain->SetYTitle("Momentum at TRD");
 295   fAbsoluteGain->SetXTitle("charge deposit [a.u]");
 296   fAbsoluteGain->SetZTitle("counts");
 297   fAbsoluteGain->SetStats(0);
 298   fAbsoluteGain->Sumw2();
 299   fListHist->Add(fAbsoluteGain);
 300
 301
 302
 303   /////////////////////////////////////////
 304   // First debug level
 305   ///////////////////////////////////////
 306   if(fDebug > 0) {
 307
 308     // Standart with AliESDfriend
 309     fPH2dSM = new TProfile2D("PH2dSM","Nz10Nrphi10"
 310                             ,fNbTimeBins,-0.05,(Double_t)((fNbTimeBins-0.5)/10.0)
 311                            ,18,0,18);
 312     fPH2dSM->SetYTitle("Det/pad groups");
 313     fPH2dSM->SetXTitle("time [#mus]");
 314     fPH2dSM->SetZTitle("<PH> [a.u.]");
 315     fPH2dSM->SetStats(0);
 316     //
 317     fCH2dSM = new TH2I("CH2dSM","Nz10Nrphi10",50,0,300,18,0,18);
 318     fCH2dSM->SetYTitle("Det/pad groups");
 319     fCH2dSM->SetXTitle("charge deposit [a.u]");
 320     fCH2dSM->SetZTitle("counts");
 321     fCH2dSM->SetStats(0);
 322     fCH2dSM->Sumw2();
 323     //
 324     fPH2dSum = new TProfile2D("PH2dSum","Nz100Nrphi100"
 325                             ,fNbTimeBins,-0.05,(Double_t)((fNbTimeBins-0.5)/10.0)
 326                             ,1,0,1);
 327     fPH2dSum->SetYTitle("Det/pad groups");
 328     fPH2dSum->SetXTitle("time [#mus]");
 329     fPH2dSum->SetZTitle("<PH> [a.u.]");
 330     fPH2dSum->SetStats(0);
 331     //
 332     fCH2dSum = new TH2I("CH2dSum","Nz100Nrphi100",50,0,300,1,0,1);
 333     fCH2dSum->SetYTitle("Det/pad groups");
 334     fCH2dSum->SetXTitle("charge deposit [a.u]");
 335     fCH2dSum->SetZTitle("counts");
 336     fCH2dSum->SetStats(0);
 337     fCH2dSum->Sumw2();
 338
 339     // Add them
 340     fListHist->Add(fPH2dSM);
 341     fListHist->Add(fCH2dSM);
 342     fListHist->Add(fPH2dSum);
 343     fListHist->Add(fCH2dSum);
 344   }
 345
 346   /////////////////////////////////////////
 347   // Second debug level
 348   ///////////////////////////////////////
 349   if(fDebug > 1) {
 350
 351     fNbTRDTrack = new TH1F("TRDTrack","TRDTrack",50,0,50);
 352     fNbTRDTrack->Sumw2();
 353     fNbTRDTrackOffline = new TH1F("TRDTrackOffline","TRDTrackOffline",50,0,50);
 354     fNbTRDTrackOffline->Sumw2();
 355     fNbTRDTrackStandalone = new TH1F("TRDTrackStandalone","TRDTrackStandalone",50,0,50);
 356     fNbTRDTrackStandalone->Sumw2();
 357     fNbTPCTRDtrack = new TH2F("NbTPCTRDtrack","NbTPCTRDtrack",100,0,100,100,0,100);
 358     fNbTPCTRDtrack->Sumw2();
 359     //
 360     fNbTimeBin = new TH1F("NbTimeBin","NbTimeBin",35,0,35);
 361     fNbTimeBin->Sumw2();
 362     fNbTimeBinOffline = new TH1F("NbTimeBinOffline","NbTimeBinOffline",35,0,35);
 363     fNbTimeBinOffline->Sumw2();
 364     fNbTimeBinStandalone = new TH1F("NbTimeBinStandalone","NbTimeBinStandalone",35,0,35);
 365     fNbTimeBinStandalone->Sumw2();
 366     //
 367     fNbClusters = new TH1F("NbClusters","",35,0,35);
 368     fNbClusters->Sumw2();
 369     fNbClustersOffline = new TH1F("NbClustersOffline","",35,0,35);
 370     fNbClustersOffline->Sumw2();
 371     fNbClustersStandalone = new TH1F("NbClustersStandalone","",35,0,35);
 372     fNbClustersStandalone->Sumw2();
 373     //
 374     fNbTracklets = new TH1F("NbTracklets","NbTracklets",540,0.,540.);
 375     fNbTracklets->Sumw2();
 376     fNbTrackletsOffline = new TH1F("NbTrackletsOffline","NbTrackletsOffline",540,0.,540.);
 377     fNbTrackletsOffline->Sumw2();
 378     fNbTrackletsStandalone = new TH1F("NbTrackletsStandalone","NbTrackletsStandalone",540,0.,540.);
 379     fNbTrackletsStandalone->Sumw2();
 380
 381     fListHist->Add(fNbTRDTrack);
 382     fListHist->Add(fNbTRDTrackOffline);
 383     fListHist->Add(fNbTRDTrackStandalone);
 384     fListHist->Add(fNbTPCTRDtrack);
 385
 386     fListHist->Add(fNbTimeBin);
 387     fListHist->Add(fNbTimeBinOffline);
 388     fListHist->Add(fNbTimeBinStandalone);
 389     fListHist->Add(fNbClusters);
 390     fListHist->Add(fNbClustersOffline);
 391     fListHist->Add(fNbClustersStandalone);
 392     fListHist->Add(fNbTracklets);
 393     fListHist->Add(fNbTrackletsOffline);
 394     fListHist->Add(fNbTrackletsStandalone);
 395
 396   }
 397
 398  }
 399  //________________________________________________________________________
 400  void AliTRDCalibTask::Exec(Option_t *)
 401  {
 402    //
 403    // Filling of the histos
 404    //
 405
 406    AliLog::SetGlobalLogLevel(AliLog::kError);
 407
 408    if (!fESD) {
 409      //Printf("ERROR: fESD not available");
 410      return;
 411    }
 412
 413    //printf("Counter %d\n",fCounter);
 414
 415    fCounter++;
 416    if((fMaxEvent != 0) && (fMaxEvent < fCounter)) return;
 417
 418    //printf("Counter %d\n",fCounter);
 419
 420    ///////////////////////////////
 421    // Require a primary vertex
 422    //////////////////////////////
 423    if(fRequirePrimaryVertex) {
 424      const AliESDVertex* vtxESD = 0x0;
 425      if      (fVtxTPC) vtxESD = fESD->GetPrimaryVertexTPC() ;
 426      else if (fVtxSPD) vtxESD = fESD->GetPrimaryVertexSPD() ;
 427      else              vtxESD = fESD->GetPrimaryVertexTracks() ;
 428      if(!vtxESD){
 429        return;
 430      }
 431      Int_t nCtrb = vtxESD->GetNContributors();
 432      if(nCtrb < fMinNbContributors) return;
 433      Double_t zPosition = vtxESD->GetZ();
 434      if(TMath::Abs(zPosition) > fRangePrimaryVertexZ) return;
 435
 436    }
 437
 438    //printf("Primary vertex passed\n");
 439
 440    ///////////////////
 441    // Check trigger
 442    ///////////////////
 443    //printf("Class Fired %s\n",(const char*)fESD->GetFiredTriggerClasses());
 444    Bool_t pass = kTRUE;
 445    Int_t numberOfTriggerSelected = fSelectedTrigger->GetEntriesFast();
 446    //printf("numberofTriggerSelected %d\n",numberOfTriggerSelected);
 447    if(fRejected) {
 448      pass = kTRUE;
 449      for(Int_t k = 0; k < numberOfTriggerSelected; k++){
 450        const TObjString *const obString=(TObjString*)fSelectedTrigger->At(k);
 451        const TString tString=obString->GetString();
 452        if(fESD->IsTriggerClassFired((const char*)tString)) {
 453          pass = kFALSE;
 454        }
 455      }
 456    }
 457    else {
 458      pass = kFALSE;
 459      for(Int_t k = 0; k < numberOfTriggerSelected; k++){
 460        const TObjString *const obString=(TObjString*)fSelectedTrigger->At(k);
 461        const TString tString=obString->GetString();
 462        if(fESD->IsTriggerClassFired((const char*)tString)) {
 463          pass = kTRUE;
 464        }
 465      }
 466    }
 467    if(!pass) return;
 468
 469    //printf("Trigger passed\n");
 470
 471    fNEvents->Fill(1);
 472
 473    // In total
 474    Int_t nbTrdTracks = 0;
 475    // standalone
 476    Int_t nbTrdTracksStandalone = 0;
 477    // offline
 478    Int_t nbTrdTracksOffline = 0;
 479    // TPC
 480    Int_t nbtrackTPC = 0;
 481
 482    Double_t nbTracks = fESD->GetNumberOfTracks();
 483    //printf("Number of tracks %f\n",nbTracks);
 484
 485    if (nbTracks <= 0.0) {
 486
 487      if(fDebug > 1) {
 488        fNbTRDTrack->Fill(nbTrdTracks);
 489        fNbTRDTrackStandalone->Fill(nbTrdTracksStandalone);
 490        fNbTRDTrackOffline->Fill(nbTrdTracksOffline);
 491      }
 492      PostData(0, fListHist);
 493      return;
 494    }
 495
 496    fESDfriend = (AliESDfriend *)fESD->FindListObject("AliESDfriend");
 497    fESD->SetESDfriend(fESDfriend);
 498    if(!fESDfriend) return;
 499
 500    //printf("has friends\n");
 501
 502    ////////////////////////////////////
 503    // Check the number of TPC tracks
 504    ///////////////////////////////////
 505    //printf("Nb of tracks %f\n",nbTracks);
 506    for(Int_t itrk = 0; itrk < nbTracks; itrk++){
 507      // Get ESD track
 508      fkEsdTrack = fESD->GetTrack(itrk);
 509      ULong_t status = fkEsdTrack->GetStatus();
 510      if(status&(AliESDtrack::kTPCout)) nbtrackTPC++;
 511      // Check that the calibration object is here
 512      if(fESDfriend && (fESDfriend->GetTrack(itrk))) {
 513        fFriendTrack = new AliESDfriendTrack(*(fESDfriend->GetTrack(itrk)));
 514        Int_t counteer = 0;
 515        Int_t icalib=0;
 516        while((fCalibObject = (TObject *)(fFriendTrack->GetCalibObject(icalib++)))){
 517          //printf("Name %s\n",fCalibObject->IsA()->GetName());
 518          if(strcmp(fCalibObject->IsA()->GetName(), "AliTRDtrackV1") != 0) continue;
 519          counteer++;
 520        }
 521        //printf("TRDntracklets %d, TRDntrackletsPID %d\n",fkEsdTrack->GetTRDntracklets(),fkEsdTrack->GetTRDntrackletsPID());
 522        if(counteer > 0) {
 523          nbTrdTracks++;
 524          if((status&(AliESDtrack::kTRDout)) && (!(status&(AliESDtrack::kTRDin)))) {
 525            nbTrdTracksStandalone++;
 526          }
 527          if((status&(AliESDtrack::kTRDin))) {
 528            nbTrdTracksOffline++;
 529          }
 530        }
 531        if(fFriendTrack) delete fFriendTrack;
 532      }
 533    }
 534    //printf("Number of TPC tracks %d, TRD %d\n",nbtrackTPC,nbTrdTracks);
 535    if((nbtrackTPC>0) && (nbTrdTracks > (3.0*nbtrackTPC))) pass = kFALSE;
 536
 537    if(fDebug > 1) {
 538
 539      fNbTRDTrack->Fill(nbTrdTracks);
 540      fNbTRDTrackStandalone->Fill(nbTrdTracksStandalone);
 541      fNbTRDTrackOffline->Fill(nbTrdTracksOffline);
 542      fNbTPCTRDtrack->Fill(nbTrdTracks,nbtrackTPC);
 543
 544    }
 545
 546    if(!pass) {
 547      PostData(0, fListHist);
 548      return;
 549    }
 550    //printf("Pass\n");
 551
 552    /////////////////////////////////////
 553    // Loop on AliESDtrack
 554    ////////////////////////////////////
 555    //printf("Nb of tracks %f\n",nbTracks);
 556    for(int itrk=0; itrk < nbTracks; itrk++){
 557
 558      // Get ESD track
 559      fkEsdTrack = fESD->GetTrack(itrk);
 560      if(!fkEsdTrack) continue;
 561
 562      // Quality cuts on the AliESDtrack
 563      if((fEsdTrackCuts) && (!fEsdTrackCuts->IsSelected((AliVParticle *)fkEsdTrack))) {
 564        //printf("Not a good track\n");
 565        continue;
 566      }
 567
 568      // First Absolute gain calibration
 569      Int_t trdNTracklets = (Int_t) fkEsdTrack->GetTRDntracklets();
 570      Int_t trdNTrackletsPID = (Int_t) fkEsdTrack->GetTRDntrackletsPID();
 571      if((trdNTracklets > 0) && (trdNTrackletsPID > 0)) {
 572        for(Int_t iPlane = 0; iPlane < 6; iPlane++){
 573          Double_t slide = fkEsdTrack->GetTRDslice(iPlane);
 574          //printf("Number of slide %d\n",fkEsdTrack->GetNumberOfTRDslices());
 575          Double_t momentum = fkEsdTrack->GetTRDmomentum(iPlane);
 576          //printf("momentum %f, slide %f\n",momentum,slide);
 577          if(slide > 0.0) fAbsoluteGain->Fill(slide*8.0/100.0,momentum);
 578        }
 579      }
 580
 581      // Other cuts
 582      Bool_t good = kTRUE;
 583      Bool_t standalonetrack = kFALSE;
 584      Bool_t offlinetrack = kFALSE;
 585      ULong_t status = fkEsdTrack->GetStatus();
 586
 587      if(!(fESDfriend->GetTrack(itrk)))  continue;
 588
 589      fFriendTrack = new AliESDfriendTrack(*(fESDfriend->GetTrack(itrk)));
 590
 591      //////////////////////////////////////
 592      // Loop on calibration objects
 593      //////////////////////////////////////
 594      Int_t icalib=0;
 595      while((fCalibObject = (TObject *)(fFriendTrack->GetCalibObject(icalib++)))){
 596        if(strcmp(fCalibObject->IsA()->GetName(), "AliTRDtrackV1") != 0) continue;
 597        //printf("Find the calibration object\n");
 598
 599        if((status&(AliESDtrack::kTRDout)) && (!(status&(AliESDtrack::kTRDin)))) {
 600          standalonetrack = kTRUE;
 601        }
 602        if((status&(AliESDtrack::kTRDin))) {
 603          offlinetrack = kTRUE;
 604        }
 605        if(fOfflineTracks){
 606          if(!offlinetrack){
 607            good = kFALSE;
 608          }
 609        }
 610        else if(fStandaloneTracks){
 611          if(!standalonetrack){
 612            good = kFALSE;
 613          }
 614        }
 615
 616        fTrdTrack = (AliTRDtrackV1 *)fCalibObject;
 617        if(good) {
 618          fTRDCalibraFillHisto->UpdateHistogramsV1(fTrdTrack);
 619          //printf("Fill fTRDCalibraFillHisto\n");
 620        }
 621
 622        //////////////////////////////////
 623        // Debug
 624        ////////////////////////////////
 625
 626        if(fDebug > 0) {
 627
 628          //printf("Enter debug\n");
 629
 630          Int_t nbtracklets = 0;
 631
 632          // Check some stuff
 633          Bool_t standalonetracklet = kFALSE;
 634          const AliTRDseedV1 *tracklet = 0x0;
 635          //////////////////////////////////////
 636          // Loop tracklets
 637          /////////////////////////////////////
 638          for(Int_t itr = 0; itr < 6; itr++){
 639
 640            if(!(tracklet = fTrdTrack->GetTracklet(itr))) continue;
 641            if(!tracklet->IsOK()) continue;
 642            nbtracklets++;
 643            standalonetracklet = kFALSE;
 644            if(tracklet->IsStandAlone()) standalonetracklet = kTRUE;
 645
 646            Int_t nbclusters = 0;
 647            Double_t phtb[AliTRDseedV1::kNtb];
 648            memset(phtb, 0, AliTRDseedV1::kNtb*sizeof(Double_t));
 649            Double_t sum = 0.0;
 650            Float_t normalisation = 6.67;
 651            Int_t detector = 0;
 652            Int_t sector = 0;
 653            //Int_t crossrow = 0;
 654
 655            // Check no shared clusters
 656            //for(int icc=AliTRDseedV1::kNtb; icc<AliTRDseedV1::kNclusters; icc++){
 657            //  if((fcl = tracklet->GetClusters(icc)))  crossrow = 1;
 658            // }
 659
 660            // Loop on clusters
 661            for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
 662
 663              if(!(fCl = tracklet->GetClusters(ic))) continue;
 664              nbclusters++;
 665              Int_t time = fCl->GetPadTime();
 666              Float_t ch =  tracklet->GetdQdl(ic);
 667              Float_t qcl = TMath::Abs(fCl->GetQ());
 668              detector = fCl->GetDetector();
 669              // Add the charge if shared cluster
 670              if((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) {
 671                if((fCl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb))) {
 672                  qcl += TMath::Abs(fCl->GetQ());
 673                  //printf("Add the cluster charge\n");
 674                }
 675              }
 676              if((time>-1) && (time<fNbTimeBins)) phtb[time]=qcl;
 677              if((fCalDetGain) && (fCalDetGain->GetValue(detector) > 0.0)) sum += ch*fCalDetGain->GetValue(detector)/normalisation;
 678              else sum += ch/normalisation;
 679
 680              if(fDebug > 1) {
 681                fNbTimeBin->Fill(time);
 682                if(tracklet->IsStandAlone()) fNbTimeBinStandalone->Fill(time);
 683                else fNbTimeBinOffline->Fill(time);
 684              }
 685            }
 686            sector = AliTRDgeometry::GetSector(detector);
 687
 688            if(fDebug > 1) {
 689              fNbTracklets->Fill(detector);
 690              if(tracklet->IsStandAlone()) fNbTrackletsStandalone->Fill(detector);
 691              else fNbTrackletsOffline->Fill(detector);
 692
 693              fNbClusters->Fill(nbclusters);
 694              if(tracklet->IsStandAlone())  fNbClustersStandalone->Fill(nbclusters);
 695              else  fNbClustersOffline->Fill(nbclusters);
 696            }
 697
 698            if(fDebug > 0) {
 699              if((nbclusters > fLow) && (nbclusters < fHigh)){
 700                if(fRelativeScale > 0.0) sum = sum/fRelativeScale;
 701                fCH2dSM->Fill(sum,sector+0.5);
 702                fCH2dSum->Fill(sum,0.5);
 703                Bool_t checknoise = kTRUE;
 704                if(fMaxCluster > 0) {
 705                  if(phtb[0] > fMaxCluster) checknoise = kFALSE;
 706                  if(fNbTimeBins > fNbMaxCluster) {
 707                    for(Int_t k = (fNbTimeBins-fNbMaxCluster); k < fNbTimeBins; k++){
 708                      if(phtb[k] > fMaxCluster) checknoise = kFALSE;
 709                    }
 710                  }
 711                }
 712                if(checknoise) {
 713                  for(int ic=0; ic<fNbTimeBins; ic++){
 714                    if(fFillZero) {
 715                      fPH2dSum->Fill((Double_t)(ic/10.0),0.5,(Double_t)phtb[ic]);
 716                      fPH2dSM->Fill((Double_t)(ic/10.0),sector+0.5,(Double_t)phtb[ic]);
 717                    }
 718                    else {
 719                      if(phtb[ic] > 0.0) {
 720                        fPH2dSum->Fill((Double_t)(ic/10.0),0.0,(Double_t)phtb[ic]);
 721                        fPH2dSM->Fill((Double_t)(ic/10.0),sector+0.5,(Double_t)phtb[ic]);
 722                      }
 723                    }
 724                  }
 725                }
 726              }
 727            }
 728          } // loop on tracklets
 729
 730        } // debug
 731
 732      }// while calibration objects
 733
 734      delete fFriendTrack;
 735
 736    } // loop ESD track
 737
 738    // Post output data
 739    PostData(0, fListHist);
 740  }
 741 //________________________________________________________________________
 742 void AliTRDCalibTask::Terminate(Option_t *)
 743 {
 744   //
 745   // Terminate
 746   //
 747
 748   if(fTRDCalibraFillHisto) fTRDCalibraFillHisto->DestroyDebugStreamer();
 749
 750
 751 }
 752 //_______________________________________________________
 753 Bool_t AliTRDCalibTask::Load(const Char_t *filename)
 754 {
 755   //
 756   // Generic container loader
 757   //
 758
 759   if(!TFile::Open(filename)){
 760     //AliWarning(Form("Couldn't open file %s.", filename));
 761     return kFALSE;
 762   }
 763   TList *o = 0x0;
 764   if(!(o = (TList*)gFile->Get(GetName()))){
 765     //AliWarning("Missing histogram container.");
 766     return kFALSE;
 767   }
 768   fListHist = (TList*)o->Clone(GetName());
 769   gFile->Close();
 770   return kTRUE;
 771 }
 772 //________________________________________________________________________
 773 void AliTRDCalibTask::Plot()
 774 {
 775   //
 776   // Plot the histos stored in the TList
 777   //
 778
 779   if(!fListHist) return;
 780
 781   /////////////////////////////////////
 782   // Take the debug stuff
 783   /////////////////////////////////////
 784
 785   TH1I *nEvents  = (TH1I *) fListHist->FindObject("NEvents");
 786
 787   TH2F *absoluteGain  = (TH2F *) fListHist->FindObject("AbsoluteGain");
 788
 789   TH1F *trdTrack = (TH1F *) fListHist->FindObject("TRDTrack");
 790   TH1F *trdTrackOffline = (TH1F *) fListHist->FindObject("TRDTrackOffline");
 791   TH1F *trdTrackStandalone = (TH1F *) fListHist->FindObject("TRDTrackStandalone");
 792
 793   TH2F *tpctrdTrack = (TH2F *) fListHist->FindObject("NbTPCTRDtrack");
 794
 795   TH1F *nbTimeBin = (TH1F *) fListHist->FindObject("NbTimeBin");
 796   TH1F *nbTimeBinOffline = (TH1F *) fListHist->FindObject("NbTimeBinOffline");
 797   TH1F *nbTimeBinStandalone = (TH1F *) fListHist->FindObject("NbTimeBinStandalone");
 798
 799   TH1F *nbClusters = (TH1F *) fListHist->FindObject("NbClusters");
 800   TH1F *nbClustersOffline = (TH1F *) fListHist->FindObject("NbClustersOffline");
 801   TH1F *nbClustersStandalone = (TH1F *) fListHist->FindObject("NbClustersStandalone");
 802
 803   TH1F *nbTracklets = (TH1F *) fListHist->FindObject("NbTracklets");
 804   TH1F *nbTrackletsOffline = (TH1F *) fListHist->FindObject("NbTrackletsOffline");
 805   TH1F *nbTrackletsStandalone = (TH1F *) fListHist->FindObject("NbTrackletsStandalone");
 806
 807   /////////////////////////////////////
 808   // Take the calibration objects
 809   /////////////////////////////////////
 810
 811   TH2I *ch2d = (TH2I *) fListHist->FindObject("CH2d");
 812   TProfile2D *ph2d = (TProfile2D *) fListHist->FindObject("PH2d");
 813
 814   TH2I *ch2dSum = (TH2I *) fListHist->FindObject("CH2dSum");
 815   TProfile2D *ph2dSum = (TProfile2D *) fListHist->FindObject("PH2dSum");
 816
 817   TH2I *ch2dSM = (TH2I *) fListHist->FindObject("CH2dSM");
 818   TProfile2D *ph2dSM = (TProfile2D *) fListHist->FindObject("PH2dSM");
 819
 820   AliTRDCalibraVdriftLinearFit *linearfit = (AliTRDCalibraVdriftLinearFit *) fListHist->FindObject("AliTRDCalibraVdriftLinearFit");
 821
 822   ////////////////////////////////////////////////
 823   // Add the AliTRDCalibraVdriftLinearFit
 824   ///////////////////////////////////////////////
 825
 826   Int_t first = 0;
 827   TH2S *histolinearfitsum = 0x0;
 828
 829   if(linearfit) {
 830     for(Int_t det = 0; det < 540; det++) {
 831       if(linearfit->GetLinearFitterHisto(det)){
 832         if(TMath::Abs(first)<0.0001){
 833           histolinearfitsum = linearfit->GetLinearFitterHisto(det);
 834           first += 1;
 835         }
 836         else {
 837           histolinearfitsum ->Add(linearfit->GetLinearFitterHisto(det));
 838         }
 839       }
 840     }
 841   }
 842
 843   ///////////////////////////
 844   // Style
 845   //////////////////////////
 846
 847   gStyle->SetPalette(1);
 848   gStyle->SetOptStat(1111);
 849   gStyle->SetOptFit(1111);
 850   gStyle->SetPadBorderMode(0);
 851   gStyle->SetCanvasColor(10);
 852   gStyle->SetPadLeftMargin(0.13);
 853   gStyle->SetPadRightMargin(0.13);
 854
 855   /////////////////////////
 856   // Plot
 857   ////////////////////////
 858
 859  if(nEvents) {
 860
 861     TCanvas *debugEvents = new TCanvas("cNEvents","cNEvents",10,10,510,510);
 862     debugEvents->cd(1);
 863     if(nEvents) nEvents->Draw();
 864
 865   }
 866
 867  if(absoluteGain) {
 868
 869     TCanvas *debugAbsoluteGain = new TCanvas("cAbsoluteGain","cAbsoluteGain",10,10,510,510);
 870     debugAbsoluteGain->cd(1);
 871     if(absoluteGain) absoluteGain->Draw();
 872
 873   }
 874
 875   if(trdTrack || tpctrdTrack) {
 876
 877     TCanvas *debugtrdtpcTrack = new TCanvas("TRDtracktpctrdtrack","TRDtracktpctrdtrack",10,10,510,510);
 878     debugtrdtpcTrack->Divide(2,1);
 879     debugtrdtpcTrack->cd(1);
 880     if(trdTrack) trdTrack->Draw();
 881     if(trdTrackOffline) trdTrackOffline->Draw("same");
 882     if(trdTrackStandalone) trdTrackStandalone->Draw("same");
 883     TLegend *leg = new TLegend(0.4,0.6,0.89,0.89);
 884     if(trdTrack) leg->AddEntry(trdTrack,"All","p");
 885     if(trdTrackOffline) leg->AddEntry(trdTrackOffline,"Offline","p");
 886     if(trdTrackStandalone) leg->AddEntry(trdTrackStandalone,"Standalone","p");
 887     leg->Draw("same");
 888     debugtrdtpcTrack->cd(2);
 889     if(tpctrdTrack) tpctrdTrack->Draw();
 890     TLine *line = new TLine(0.0,0.0,100.0,100.0);
 891     line->Draw("same");
 892
 893   }
 894
 895   if(nbTimeBin || nbTracklets || nbClusters) {
 896
 897     TCanvas *debugTracklets = new TCanvas("TRDtimebintrackletcluster","TRDtimebintrackletcluster",10,10,510,510);
 898     debugTracklets->Divide(3,1);
 899     debugTracklets->cd(1);
 900     if(nbTimeBin) nbTimeBin->Draw();
 901     if(nbTimeBinOffline) nbTimeBinOffline->Draw("same");
 902     if(nbTimeBinStandalone) nbTimeBinStandalone->Draw("same");
 903     TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
 904     if(nbTimeBin) lega->AddEntry(nbTimeBin,"All","p");
 905     if(nbTimeBinOffline) lega->AddEntry(nbTimeBinOffline,"Offline","p");
 906     if(nbTimeBinStandalone) lega->AddEntry(nbTimeBinStandalone,"Standalone","p");
 907     lega->Draw("same");
 908     debugTracklets->cd(2);
 909     if(nbTracklets) nbTracklets->Draw();
 910     if(nbTrackletsOffline) nbTrackletsOffline->Draw("same");
 911     if(nbTrackletsStandalone) nbTrackletsStandalone->Draw("same");
 912     TLegend *legb = new TLegend(0.4,0.6,0.89,0.89);
 913     if(nbTracklets) legb->AddEntry(nbTracklets,"All","p");
 914     if(nbTrackletsOffline) legb->AddEntry(nbTrackletsOffline,"Offline","p");
 915     if(nbTrackletsStandalone) legb->AddEntry(nbTrackletsStandalone,"Standalone","p");
 916     legb->Draw("same");
 917     debugTracklets->cd(3);
 918     if(nbClusters) nbClusters->Draw();
 919     if(nbClustersOffline) nbClustersOffline->Draw("same");
 920     if(nbClustersStandalone) nbClustersStandalone->Draw("same");
 921     TLegend *legc = new TLegend(0.4,0.6,0.89,0.89);
 922     if(nbClusters) legc->AddEntry(nbClusters,"All","p");
 923     if(nbClustersOffline) legc->AddEntry(nbClustersOffline,"Offline","p");
 924     if(nbClustersStandalone) legc->AddEntry(nbClustersStandalone,"Standalone","p");
 925     legc->Draw("same");
 926
 927   }
 928
 929   if(ch2dSum || ph2dSum || histolinearfitsum) {
 930
 931     TCanvas *debugSum = new TCanvas("SumCalibrationObjects","SumCalibrationObjects",10,10,510,510);
 932     debugSum->Divide(3,1);
 933     debugSum->cd(1);
 934     if(ch2dSum) ch2dSum->Draw("lego");
 935     debugSum->cd(2);
 936     if(ph2dSum) ph2dSum->Draw("lego");
 937     debugSum->cd(3);
 938     if(histolinearfitsum) histolinearfitsum->Draw();
 939
 940   }
 941
 942   if(ch2dSM || ph2dSM) {
 943
 944     TCanvas *debugSM = new TCanvas("SMCalibrationObjects","SMCalibrationObjects",10,10,510,510);
 945     debugSM->Divide(2,1);
 946     debugSM->cd(1);
 947     if(ch2dSM) ch2dSM->Draw("lego");
 948     debugSM->cd(2);
 949     if(ph2dSM) ph2dSM->Draw("lego");
 950
 951   }
 952
 953   if(ch2d || ph2d) {
 954
 955     TCanvas *debug = new TCanvas("CalibrationObjects","CalibrationObjects",10,10,510,510);
 956     debug->Divide(2,1);
 957     debug->cd(1);
 958     if(ch2d) ch2d->Draw("lego");
 959     debug->cd(2);
 960     if(ph2d) ph2d->Draw("lego");
 961
 962   }
 963
 964 }
 965
 966