2 /**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 /////////////////////////////////////////////////////////////////////////////////
21 // Task to run the calibration offline.
23 // R. Bailhache (rbailhache@ikf.uni-frankfurt.de, R.Bailhache@gsi.de)
25 //////////////////////////////////////////////////////////////////////////////////
28 #include "Riostream.h"
31 #include "TProfile2D.h"
39 #include "TObjArray.h"
43 #include "TGraphErrors.h"
45 #include "AliTRDrecoTask.h"
46 #include "AliAnalysisManager.h"
48 #include "AliESDInputHandler.h"
49 #include "AliTRDtrackV1.h"
50 #include "AliTRDseedV1.h"
51 #include "AliTRDcluster.h"
52 #include "info/AliTRDtrackInfo.h"
53 #include "AliTRDcalibDB.h"
55 #include "AliTRDCalibraFillHisto.h"
56 #include "AliTRDCalibraFit.h"
57 #include "AliTRDCalibraVdriftLinearFit.h"
58 #include "AliTRDCalibraMode.h"
59 #include "AliTRDCalibraVector.h"
60 #include "./Cal/AliTRDCalPad.h"
61 #include "./Cal/AliTRDCalDet.h"
65 #include "AliTRDcalibration.h"
68 ClassImp(AliTRDcalibration)
70 //________________________________________________________________________
71 AliTRDcalibration::AliTRDcalibration()
76 ,fTRDCalibraFillHisto(0)
78 ,fNbTRDTrackOffline(0)
79 ,fNbTRDTrackStandalone(0)
81 ,fNbTRDTrackletOffline(0)
82 ,fNbTRDTrackletStandalone(0)
84 ,fNbTimeBinOffline(0x0)
85 ,fNbTimeBinStandalone(0x0)
87 ,fNbClustersOffline(0)
88 ,fNbClustersStandalone(0)
102 ,fnormalizeNbOfCluster(kFALSE)
104 ,fOfflineTracks(kFALSE)
105 ,fStandaloneTracks(kFALSE)
106 ,fCompressPerDetector(kFALSE)
109 ,fPostProcess(kFALSE)
115 for(Int_t k = 0; k < 3; k++)
123 AliTRDcalibration::AliTRDcalibration(char* name)
124 :AliTRDrecoTask(name, "Calibration on tracks")
128 ,fTRDCalibraFillHisto(0)
130 ,fNbTRDTrackOffline(0)
131 ,fNbTRDTrackStandalone(0)
133 ,fNbTRDTrackletOffline(0)
134 ,fNbTRDTrackletStandalone(0)
136 ,fNbTimeBinOffline(0x0)
137 ,fNbTimeBinStandalone(0x0)
139 ,fNbClustersOffline(0)
140 ,fNbClustersStandalone(0)
149 ,fVdriftLinear(kTRUE)
154 ,fnormalizeNbOfCluster(kFALSE)
156 ,fOfflineTracks(kFALSE)
157 ,fStandaloneTracks(kFALSE)
158 ,fCompressPerDetector(kFALSE)
161 ,fPostProcess(kFALSE)
167 for(Int_t k = 0; k < 3; k++)
175 //________________________________________________________________________
176 AliTRDcalibration::~AliTRDcalibration()
178 // Default destructor
180 if(fNbTRDTrack) delete fNbTRDTrack;
181 if(fNbTRDTrackOffline) delete fNbTRDTrackOffline;
182 if(fNbTRDTrackStandalone) delete fNbTRDTrackStandalone;
183 if(fNbTRDTracklet) delete fNbTRDTracklet;
184 if(fNbTRDTrackletOffline) delete fNbTRDTrackletOffline;
185 if(fNbTRDTrackletStandalone) delete fNbTRDTrackletStandalone;
186 if(fNbTimeBin) delete fNbTimeBin;
187 if(fNbTimeBinOffline) delete fNbTimeBinOffline;
188 if(fNbTimeBinStandalone) delete fNbTimeBinStandalone;
189 if(fNbClusters) delete fNbClusters;
190 if(fNbClustersOffline) delete fNbClustersOffline;
191 if(fNbClustersStandalone) delete fNbClustersStandalone;
192 if(fPHSM) delete fPHSM;
193 if(fCHSM) delete fCHSM;
194 if(fPHSum) delete fPHSum;
195 if(fCHSum) delete fCHSum;
196 if(fDetSum) delete fDetSum;
197 if(fDetSumVector) delete fDetSumVector;
198 if(fGraph){fGraph->Delete(); delete fGraph;}
199 if(fArrayCalib){fArrayCalib->Delete(); delete fArrayCalib;}
202 //________________________________________________________________________
203 void AliTRDcalibration::UserCreateOutputObjects()
205 // Create output objects
207 // Number of time bins
209 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
210 fNbTimeBins = cal->GetNumberOfTimeBinsDCS();
211 if(fNbTimeBins <= 0){
212 AliWarning(Form("No of TimeBins from DB [%d] use default [30]", fNbTimeBins));
217 // instance calibration: what to calibrate
218 fTRDCalibraFillHisto = AliTRDCalibraFillHisto::Instance();
219 fTRDCalibraFillHisto->SetHisto2d(fHisto2d); // choose to use histograms
220 fTRDCalibraFillHisto->SetVector2d(fVector2d); // choose to use vectors
221 fTRDCalibraFillHisto->SetCH2dOn(); // choose to calibrate the gain
222 fTRDCalibraFillHisto->SetPH2dOn(); // choose to calibrate the drift velocity
223 fTRDCalibraFillHisto->SetPRF2dOn(); // choose to look at the PRF
224 fTRDCalibraFillHisto->SetLinearFitterOn(fVdriftLinear); // Other possibility vdrift VDRIFT
225 fTRDCalibraFillHisto->SetLinearFitterDebugOn(fVdriftLinear); // Other possibility vdrift
227 // segmentation (should be per default the max and add at the end)
228 for(Int_t k = 0; k < 3; k++){
229 if(((fNz[k] != 10) && (fNrphi[k] != 10)) && ((fNz[k] != 100) && (fNrphi[k] != 100))) {
230 fTRDCalibraFillHisto->SetNz(k,fNz[k]); // Mode calibration
231 fTRDCalibraFillHisto->SetNrphi(k,fNrphi[k]); // Mode calibration
234 if((fNz[k] == 100) && (fNrphi[k] == 100)) {
235 if(fVector2d) AliInfo("The mode all together is not supported by the vector method");
236 fTRDCalibraFillHisto->SetAllTogether(k);
238 if((fNz[k] == 10) && (fNrphi[k] == 10)) {
239 if(fVector2d) AliInfo("The mode per supermodule is not supported by the vector method");
240 fTRDCalibraFillHisto->SetPerSuperModule(k);
246 fTRDCalibraFillHisto->SetDebugLevel(DebugLevel()); //debug stuff
249 fTRDCalibraFillHisto->Init2Dhistos(fNbTimeBins); // initialise the histos
252 fTRDCalibraFillHisto->SetNumberClusters(flow); // At least flow clusters
253 fTRDCalibraFillHisto->SetNumberClustersf(fhigh); // The more fhigh clusters
254 fTRDCalibraFillHisto->SetFillWithZero(ffillZero); // Fill zeros
255 fTRDCalibraFillHisto->SetNormalizeNbOfCluster(fnormalizeNbOfCluster); // For iterations
257 // Add them to the container
258 fContainer = new TObjArray();
260 fContainer->Add(fTRDCalibraFillHisto->GetCH2d()); //TH2I
261 fContainer->Add(fTRDCalibraFillHisto->GetPH2d()); //TProfile2D
262 fContainer->Add(fTRDCalibraFillHisto->GetPRF2d()); //TProfile2D
264 if(fVdriftLinear) fContainer->Add(fTRDCalibraFillHisto->GetVdriftLinearFit()); // Other drift velocity
265 if(fVector2d) fContainer->Add(fTRDCalibraFillHisto->GetCalibraVector()); //calibra vector
269 // Init the debug histos
270 fNbTRDTrack = new TH1F("TRDTrack","TRDTrack",500,0,500);
271 fNbTRDTrack->Sumw2();
272 fNbTRDTrackOffline = new TH1F("TRDTrackOffline","TRDTrackOffline",500,0,500);
273 fNbTRDTrackOffline->Sumw2();
274 fNbTRDTrackStandalone = new TH1F("TRDTrackStandalone","TRDTrackStandalone",500,0,500);
275 fNbTRDTrackStandalone->Sumw2();
277 fNbTRDTracklet = new TH1F("TRDTracklet","TRDTracklet",540,0.,540.);
278 fNbTRDTracklet->Sumw2();
279 fNbTRDTrackletOffline = new TH1F("TRDTrackletOffline","TRDTrackletOffline",540,0.,540.);
280 fNbTRDTrackletOffline->Sumw2();
281 fNbTRDTrackletStandalone = new TH1F("TRDTrackletStandalone","TRDTrackletStandalone",540,0.,540.);
282 fNbTRDTrackletStandalone->Sumw2();
284 fNbTimeBin = new TH1F("TimeBin","TimeBin",35,0,35);
286 fNbTimeBinOffline = new TH1F("TimeBinOffline","TimeBinOffline",35,0,35);
287 fNbTimeBinOffline->Sumw2();
288 fNbTimeBinStandalone = new TH1F("TimeBinStandalone","TimeBinStandalone",35,0,35);
289 fNbTimeBinStandalone->Sumw2();
291 fNbClusters = new TH1F("NbClusters","",35,0,35);
292 fNbClusters->Sumw2();
293 fNbClustersOffline = new TH1F("NbClustersOffline","",35,0,35);
294 fNbClustersOffline->Sumw2();
295 fNbClustersStandalone = new TH1F("NbClustersStandalone","",35,0,35);
296 fNbClustersStandalone->Sumw2();
298 fPHSM = new TProfile2D("PH2dSM","Nz10Nrphi10"
299 ,fNbTimeBins,-0.05,(Double_t)((fNbTimeBins-0.5)/10.0)
301 fPHSM->SetYTitle("Det/pad groups");
302 fPHSM->SetXTitle("time [#mus]");
303 fPHSM->SetZTitle("<PH> [a.u.]");
306 fCHSM = new TH2I("CH2dSM","Nz10Nrphi10",50,0,300,18,0,18);
307 fCHSM->SetYTitle("Det/pad groups");
308 fCHSM->SetXTitle("charge deposit [a.u]");
309 fCHSM->SetZTitle("counts");
313 fPHSum = new TProfile2D("PH2dSum","Nz100Nrphi100"
314 ,fNbTimeBins,-0.05,(Double_t)((fNbTimeBins-0.5)/10.0)
316 fPHSum->SetYTitle("Det/pad groups");
317 fPHSum->SetXTitle("time [#mus]");
318 fPHSum->SetZTitle("<PH> [a.u.]");
321 fCHSum = new TH2I("CH2dSum","Nz100Nrphi100",50,0,300,1,0,1);
322 fCHSum->SetYTitle("Det/pad groups");
323 fCHSum->SetXTitle("charge deposit [a.u]");
324 fCHSum->SetZTitle("counts");
329 fContainer->Add(fNbTRDTrack);
330 fContainer->Add(fNbTRDTrackOffline);
331 fContainer->Add(fNbTRDTrackStandalone);
332 fContainer->Add(fNbTRDTracklet);
333 fContainer->Add(fNbTRDTrackletOffline);
334 fContainer->Add(fNbTRDTrackletStandalone);
335 fContainer->Add(fNbTimeBin);
336 fContainer->Add(fNbTimeBinOffline);
337 fContainer->Add(fNbTimeBinStandalone);
338 fContainer->Add(fNbClusters);
339 fContainer->Add(fNbClustersOffline);
340 fContainer->Add(fNbClustersStandalone);
341 fContainer->Add(fPHSM);
342 fContainer->Add(fCHSM);
343 fContainer->Add(fPHSum);
344 fContainer->Add(fCHSum);
348 PostData(1, fContainer);
351 //________________________________________________________________________
352 void AliTRDcalibration::UserExec(Option_t *)
355 // Execute function where the reference data are filled
361 Int_t nbTrdTracks = 0;
363 Int_t nbTrdTracksStandalone = 0;
365 Int_t nbTrdTracksOffline = 0;
369 // Loop on track in the event
371 //printf("Total of %d\n",fTracks->GetEntriesFast());
372 for(Int_t itrk=0; itrk < fTracks->GetEntriesFast(); itrk++){
374 //printf("itrk %d\n",itrk);
376 fTrackInfo = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
377 ftrdTrack = fTrackInfo->GetTrack();
378 if(!ftrdTrack) continue;
382 fTRDCalibraFillHisto->UpdateHistogramsV1(ftrdTrack);
386 Bool_t standalonetracklet = kFALSE;
387 const AliTRDseedV1 *tracklet = 0x0;
389 // Loop on tracklet in the event
391 for(Int_t itr = 0; itr < 6; itr++){
392 //printf("itr %d\n",itr);
393 if(!(tracklet = ftrdTrack->GetTracklet(itr))) continue;
394 if(!tracklet->IsOK()) continue;
396 if(tracklet->IsStandAlone()) standalonetracklet = kTRUE;
397 Int_t nbclusters = 0;
399 Double_t phtb[AliTRDseedV1::kNtb];
400 memset(phtb, 0, AliTRDseedV1::kNtb*sizeof(Double_t));
404 Float_t normalisation = 6.67;
408 // Check no shared clusters
409 for(int icc=AliTRDseedV1::kNtb; icc<AliTRDseedV1::kNclusters; icc++){
410 if((fcl = tracklet->GetClusters(icc))) crossrow = 1;
413 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
414 if(!(fcl = tracklet->GetClusters(ic))) continue;
416 Int_t time = fcl->GetPadTime();
417 Float_t ch = tracklet->GetdQdl(ic);
418 Float_t qcl = TMath::Abs(fcl->GetQ());
419 detector = fcl->GetDetector();
420 if((time>-1) && (time<fNbTimeBins)) phtb[time]=qcl;
421 sum += ch/normalisation;
422 fNbTimeBin->Fill(time);
423 if(tracklet->IsStandAlone()) fNbTimeBinStandalone->Fill(time);
424 else fNbTimeBinOffline->Fill(time);
426 sector = AliTRDgeometry::GetSector(detector);
428 fNbTRDTracklet->Fill(detector);
429 if(tracklet->IsStandAlone()) fNbTRDTrackletStandalone->Fill(detector);
430 else fNbTRDTrackletOffline->Fill(detector);
432 fNbClusters->Fill(nbclusters);
433 if(tracklet->IsStandAlone()) fNbClustersStandalone->Fill(nbclusters);
434 else fNbClustersOffline->Fill(nbclusters);
436 if((nbclusters > flow) && (nbclusters < fhigh)){
437 fCHSM->Fill(sum/20.0,sector+0.5);
438 fCHSum->Fill(sum/20.0,0.5);
439 for(int ic=0; ic<fNbTimeBins; ic++){
441 fPHSum->Fill((Double_t)(ic/10.0),0.5,(Double_t)phtb[ic]);
442 fPHSM->Fill((Double_t)(ic/10.0),sector+0.5,(Double_t)phtb[ic]);
446 fPHSum->Fill((Double_t)(ic/10.0),0.0,(Double_t)phtb[ic]);
447 fPHSM->Fill((Double_t)(ic/10.0),sector+0.5,(Double_t)phtb[ic]);
454 if(standalonetracklet) nbTrdTracksStandalone++;
455 else nbTrdTracksOffline++;
464 fNbTRDTrack->Fill(nbTrdTracks);
465 fNbTRDTrackStandalone->Fill(nbTrdTracksStandalone);
466 fNbTRDTrackOffline->Fill(nbTrdTracksOffline);
471 //________________________________________________________________________
472 void AliTRDcalibration::Terminate(Option_t *)
474 // Draw result to the screen
475 // Called once at the end of the query
477 //printf("terminate\n");
479 if(fTRDCalibraFillHisto) fTRDCalibraFillHisto->DestroyDebugStreamer();
482 //________________________________________________________
483 Bool_t AliTRDcalibration::GetRefFigure(Int_t ifig)
486 // Draw filled histos
489 gStyle->SetPalette(1);
490 gStyle->SetOptStat(1111);
491 gStyle->SetPadBorderMode(0);
492 gStyle->SetCanvasColor(10);
493 gStyle->SetPadLeftMargin(0.13);
494 gStyle->SetPadRightMargin(0.13);
496 if(!fContainer) return kFALSE;
500 TCanvas *c0 = new TCanvas("c0","c0",10,10,510,510);
501 TLegend *legNbTrack = new TLegend(.7, .7, .98, .98);
502 legNbTrack->SetBorderSize(1);
506 if(!(h = (TH1F *)fContainer->FindObject("TRDTrack"))) break;
507 if(!(ho = (TH1F *)fContainer->FindObject("TRDTrackOffline"))) break;
508 if(!(hs = (TH1F *)fContainer->FindObject("TRDTrackStandalone"))) break;
516 legNbTrack->AddEntry(h, "all", "p");
517 legNbTrack->AddEntry(ho, "offline", "p");
518 legNbTrack->AddEntry(hs, "standalone", "p");
519 legNbTrack->Draw("same");
523 TLegend *legNbTracklet = new TLegend(.7, .7, .98, .98);
524 legNbTracklet->SetBorderSize(1);
528 if(!(h = (TH1F *)fContainer->FindObject("TRDTracklet"))) break;
529 if(!(ho = (TH1F *)fContainer->FindObject("TRDTrackletOffline"))) break;
530 if(!(hs = (TH1F *)fContainer->FindObject("TRDTrackletStandalone"))) break;
534 legNbTracklet->AddEntry(h, "all", "p");
535 legNbTracklet->AddEntry(ho, "offline", "p");
536 legNbTracklet->AddEntry(hs, "standalone", "p");
537 legNbTracklet->Draw("same");
544 TLegend *legNbTimeBin = new TLegend(.7, .7, .98, .98);
545 legNbTimeBin->SetBorderSize(1);
549 if(!(h = (TH1F *)fContainer->FindObject("TimeBin"))) break;
550 if(!(ho = (TH1F *)fContainer->FindObject("TimeBinOffline"))) break;
551 if(!(hs = (TH1F *)fContainer->FindObject("TimeBinStandalone"))) break;
555 legNbTimeBin->AddEntry(h, "all", "p");
556 legNbTimeBin->AddEntry(ho, "offline", "p");
557 legNbTimeBin->AddEntry(hs, "standalone", "p");
558 legNbTimeBin->Draw("same");
565 TLegend *legNbClusters = new TLegend(.7, .7, .98, .98);
566 legNbClusters->SetBorderSize(1);
570 if(!(h = (TH1F *)fContainer->FindObject("NbClusters"))) break;
571 if(!(ho = (TH1F *)fContainer->FindObject("NbClustersOffline"))) break;
572 if(!(hs = (TH1F *)fContainer->FindObject("NbClustersStandalone"))) break;
576 legNbClusters->AddEntry(h, "all", "p");
577 legNbClusters->AddEntry(ho, "offline", "p");
578 legNbClusters->AddEntry(hs, "standalone", "p");
579 legNbClusters->Draw("same");
587 if(!(h = (TProfile2D *)fContainer->FindObject("PH2dSum"))) break;
588 TH1D *projh = h->ProjectionX("projh",1,1,"e");
596 if(!(h = (TH2I *)fContainer->FindObject("CH2dSum"))) break;
597 TH1D *projh = h->ProjectionX("projhh",1,1,"e");
605 AliInfo("Histo was not filled!");
609 if(!(h = (TProfile2D *)fContainer->FindObject("PH2d"))) break;
615 AliInfo("Histo was not filled!");
619 if(!(h = (TH2I *)fContainer->FindObject("CH2d"))) break;
625 AliInfo("Histo was not filled!");
629 if(!(h = (TProfile2D *)fContainer->FindObject("PRF2d"))) break;
635 AliInfo("vector was not filled!");
638 AliTRDCalibraVector *v = 0x0;
639 TGraphErrors *vdet = 0x0;
640 if(!(v = (AliTRDCalibraVector *)fContainer->FindObject("AliTRDCalibraVector"))) break;
641 Int_t detectormax = -1;
643 if(!v->FindTheMaxEntries(1,detectormax,groupmax)) break;
644 if(!(vdet = v->ConvertVectorPHTGraphErrors((Int_t)detectormax,groupmax,"plotPH2dVector"))) break;
645 Int_t nbeentries = 0;
646 TH1F *ko = v->CorrectTheError(vdet,nbeentries);
648 AliInfo(Form("There are %d entries in the detector %d and group %d",nbeentries,detectormax,groupmax));
653 AliInfo("vector was not filled!");
656 AliTRDCalibraVector *v = 0x0;
658 if(!(v = (AliTRDCalibraVector *)fContainer->FindObject("AliTRDCalibraVector"))) break;
659 Int_t detectormax = -1;
661 if(!v->FindTheMaxEntries(0,detectormax,groupmax)) break;
662 if(!(vdet = v->ConvertVectorCHHisto((Int_t)detectormax,groupmax,"plotCH2dVector"))) break;
664 AliInfo(Form("The detectormax and groupmax are %d and %d",detectormax,groupmax));
669 AliInfo("vector was not filled!");
672 AliTRDCalibraVector *v = 0x0;
673 TGraphErrors *vdet = 0x0;
674 if(!(v = (AliTRDCalibraVector *)fContainer->FindObject("AliTRDCalibraVector"))) break;
675 Int_t detectormax = -1;
677 Int_t nbeentries = 0;
678 if(!v->FindTheMaxEntries(2,detectormax,groupmax)) break;
679 if(!(vdet = v->ConvertVectorPRFTGraphErrors((Int_t)detectormax,groupmax,"plotPRF2dVector"))) break;
680 TH1F *ko = v->CorrectTheError(vdet,nbeentries);
682 AliInfo(Form("The detectormax and groupmax are %d and %d",detectormax,groupmax));
687 AliInfo("vdrift linear was not filled!");
690 AliTRDCalibraVdriftLinearFit *h = 0x0;
691 TH2S *hdetector = 0x0;
692 if(!(h = (AliTRDCalibraVdriftLinearFit *)fContainer->FindObject("AliTRDCalibraVdriftLinearFit"))) break;
693 Double_t entries[540];
694 for(Int_t k = 0; k < 540; k++){
697 if(!(hdetector = (TH2S *)h->GetLinearFitterHisto(k,kFALSE))) continue;
698 entries[k] = hdetector->GetEntries();
700 Double_t max = -10.0;
701 Double_t detectormax = -1;
702 for(Int_t k = 0; k < 540; k++){
703 if(entries[k] > max) {
709 if((TMath::Abs(max) <= 0.001) || (detectormax <0.0) || (detectormax >=540.0)) break;
710 if(!(hdetector = (TH2S *)h->GetLinearFitterHisto((Int_t)detectormax,kFALSE))) break;
711 AliInfo(Form("The detector with the maximum of entries is %f",detectormax));
717 if(!PostProcess()) break;
719 TGraph *fgain = (TGraph *) fGraph->At(0);
724 case kVdriftT0Factor:{
726 if(!PostProcess()) break;
728 TCanvas *c = new TCanvas("c","c",10,10,510,510);
730 TGraph *fvd = (TGraph *) fGraph->At(1);
735 TGraph *ft0 = (TGraph *) fGraph->At(2);
742 case kVdriftLorentzAngleFactor:{
744 AliInfo("vdrift linear was not filled!");
748 if(!PostProcess()) break;
750 TCanvas *c = new TCanvas("c","c",10,10,510,510);
752 TGraph *fvdl = (TGraph *) fGraph->At(3);
757 TGraph *flal = (TGraph *) fGraph->At(4);
766 if(!PostProcess()) break;
768 TGraph *fprf = (TGraph *) fGraph->At(5);
778 //________________________________________________________________________
779 Bool_t AliTRDcalibration::PostProcess()
782 // Fit the filled histos
783 // Put the calibration object in fArrayCalib
784 // 0 and 1 AliTRDCalDet and AliTRDCalPad gain
785 // 2 and 3 AliTRDCalDet and AliTRDCalPad vdrift PH
786 // 4 and 5 AliTRDCalDet and AliTRDCalPad timeoffset PH
787 // 6 AliTRDCalPad PRF
788 // 7 and 8 AliTRDCalDet and AliTRDCalPad vdrift second
789 // 9 and 10 AliTRDCalDet and AliTRDCalPad lorentz angle second
793 fArrayCalib = new TObjArray(11);
794 fArrayCalib->SetOwner();
802 fGraph = new TObjArray(6);
810 Bool_t storage[3] = {kFALSE,kFALSE,kFALSE};
812 // Objects for fitting
813 AliTRDCalibraFit *calibra = AliTRDCalibraFit::Instance();
814 calibra->SetDebugLevel(2); // 0 rien, 1 fitvoir, 2 debug files, 3 one detector
818 Printf("ERROR: list not available");
822 if(fHisto2d && fVector2d) AliInfo("We will only look at histos. Set fHisto2d off if you don't want");
823 AliTRDCalibraVector *calibraVector = 0x0;
824 if(fVector2d) calibraVector = (AliTRDCalibraVector *) fContainer->FindObject("CalibraVector");
828 Bool_t pass = kFALSE;
829 AliTRDCalibraVector *vvect = 0x0;
831 TH2I *histogain = (TH2I *) fContainer->FindObject("CH2d");
833 histogain->SetDirectory(0);
834 calibra->SetMinEntries(20);
835 if(fCompressPerDetector){
836 if(AddStatsPerDetector(histogain)) pass = calibra->AnalyseCH(fCHSum);
838 else pass = calibra->AnalyseCH(histogain);
842 if(fVector2d && calibraVector) {
843 calibra->SetMinEntries(20);
844 if(fCompressPerDetector){
845 if(!(vvect = calibraVector->AddStatsPerDetectorCH())) return kFALSE;
846 pass = calibra->AnalyseCH(vvect);
848 else pass = calibra->AnalyseCH(calibraVector);
854 Int_t nbtg = 6*4*18*((Int_t) ((AliTRDCalibraMode *)calibra->GetCalibraMode())->GetDetChamb0(0))
855 + 6* 18*((Int_t) ((AliTRDCalibraMode *)calibra->GetCalibraMode())->GetDetChamb2(0));
856 Int_t nbfit = calibra->GetNumberFit(); //Number of fits
857 Int_t nbE = calibra->GetNumberEnt(); //Number of detector mit entries
860 (nbfit >= 0.001*nbE))
862 // create the cal objects
863 calibra->PutMeanValueOtherVectorFit(1,kTRUE);
864 TObjArray object = calibra->GetVectorFit();
865 AliTRDCalDet *objgaindet = calibra->CreateDetObjectGain(&object);
866 TObject *objgainpad = calibra->CreatePadObjectGain();
868 fArrayCalib->AddAt(objgaindet,0);
869 fArrayCalib->AddAt(objgainpad,1);
873 if(FillGraphIndex(&object,graph)){
874 fGraph->AddAt(graph,0);
876 }//if(enough statistics?)
877 calibra->ResetVectorFit();
882 // VDRIFT average pulse height
886 TProfile2D *histodriftvelocity = (TProfile2D *) fContainer->FindObject("PH2d");
887 if(histodriftvelocity) {
888 histodriftvelocity->SetDirectory(0);
889 calibra->SetMinEntries(20*20);
890 if(fCompressPerDetector){
891 if(AddStatsPerDetector(histodriftvelocity,1)) {
892 pass = calibra->AnalysePH(fDetSumVector);
895 else pass = calibra->AnalysePH(histodriftvelocity);
899 if(fVector2d && calibraVector) {
900 calibra->SetMinEntries(20*20);
901 if(fCompressPerDetector){
902 if(!(vvect = calibraVector->AddStatsPerDetectorPH())) return kFALSE;
903 pass = calibra->AnalysePH(vvect);
905 else pass = calibra->AnalysePH(calibraVector);
910 Int_t nbtg = 6*4*18*((Int_t) ((AliTRDCalibraMode *)calibra->GetCalibraMode())->GetDetChamb0(1))
911 + 6* 18*((Int_t) ((AliTRDCalibraMode *)calibra->GetCalibraMode())->GetDetChamb2(1));
912 Int_t nbfit = calibra->GetNumberFit();
913 Int_t nbE = calibra->GetNumberEnt();
916 (nbfit >= 0.001*nbE))
918 // create the cal objects
919 calibra->PutMeanValueOtherVectorFit(1,kTRUE);
920 calibra->PutMeanValueOtherVectorFit2(1,kTRUE);
921 TObjArray object = calibra->GetVectorFit();
922 AliTRDCalDet *objdriftvelocitydet = calibra->CreateDetObjectVdrift(&object,kTRUE);
923 TObject *objdriftvelocitypad = calibra->CreatePadObjectVdrift();
924 TObjArray objectt = calibra->GetVectorFit2();
925 AliTRDCalDet *objtime0det = calibra->CreateDetObjectT0(&object,kTRUE);
926 TObject *objtime0pad = calibra->CreatePadObjectT0();
928 fArrayCalib->AddAt(objdriftvelocitydet,2);
929 fArrayCalib->AddAt(objdriftvelocitypad,3);
931 fArrayCalib->AddAt(objtime0det,4);
932 fArrayCalib->AddAt(objtime0pad,5);
935 if(FillGraphIndex(&object,graph)){
936 fGraph->AddAt(graph,1);
938 TGraph *graphh = 0x0;
939 if(FillGraphIndex(&objectt,graphh)){
940 fGraph->AddAt(graphh,2);
942 }//if(enough statistics)
943 calibra->ResetVectorFit();
952 TProfile2D *histoprf = (TProfile2D *) fContainer->FindObject("PRF2d");
954 histoprf->SetDirectory(0);
955 calibra->SetMinEntries(600);
956 if(fCompressPerDetector){
957 if(AddStatsPerDetector(histoprf,2)) pass = calibra->AnalysePRFMarianFit(fDetSumVector);
959 else pass = calibra->AnalysePRFMarianFit(histoprf);
963 if(fVector2d && calibraVector) {
964 calibra->SetMinEntries(600);
965 if(fCompressPerDetector){
966 if(!(vvect =calibraVector->AddStatsPerDetectorPRF())) return kFALSE;
967 pass = calibra->AnalysePRFMarianFit(vvect);
969 else pass = calibra->AnalysePRFMarianFit(calibraVector);
974 Int_t nbtg = 6*4*18*((Int_t) ((AliTRDCalibraMode *)calibra->GetCalibraMode())->GetDetChamb0(2))
975 + 6* 18*((Int_t) ((AliTRDCalibraMode *)calibra->GetCalibraMode())->GetDetChamb2(2));
976 Int_t nbfit = calibra->GetNumberFit();
977 Int_t nbE = calibra->GetNumberEnt();
980 (nbfit >= 0.001*nbE)) {
981 TObjArray object = calibra->GetVectorFit();
982 TObject *objPRFpad = calibra->CreatePadObjectPRF(&object);
984 fArrayCalib->AddAt(objPRFpad,6);
987 if(FillGraphIndex(&object,graph)){
988 fGraph->AddAt(graph,5);
991 calibra->ResetVectorFit();
998 AliTRDCalibraVdriftLinearFit *vlinearfit = (AliTRDCalibraVdriftLinearFit *) fContainer->FindObject("LinearVdriftFit");
1000 calibra->SetMinEntries(20*20);
1001 if(calibra->AnalyseLinearFitters(vlinearfit)) {
1003 Int_t nbtg = 6*4*18*((Int_t) ((AliTRDCalibraMode *)calibra->GetCalibraMode())->GetDetChamb0(2))
1004 + 6* 18*((Int_t) ((AliTRDCalibraMode *)calibra->GetCalibraMode())->GetDetChamb2(2));
1005 Int_t nbfit = calibra->GetNumberFit();
1006 Int_t nbE = calibra->GetNumberEnt();
1007 // enough statistics
1009 (nbfit >= 0.001*nbE))
1011 // create the cal objects
1012 calibra->PutMeanValueOtherVectorFit(1,kTRUE);
1013 calibra->PutMeanValueOtherVectorFit2(1,kTRUE);
1014 TObjArray object = calibra->GetVectorFit();
1015 AliTRDCalDet *objdriftvelocitydet = calibra->CreateDetObjectVdrift(&object,kTRUE);
1016 TObject *objdriftvelocitypad = calibra->CreatePadObjectVdrift();
1017 TObjArray objectt = calibra->GetVectorFit2();
1018 AliTRDCalDet *objtime0det = calibra->CreateDetObjectT0(&object,kTRUE);
1019 TObject *objtime0pad = calibra->CreatePadObjectT0();
1021 fArrayCalib->AddAt(objdriftvelocitydet,7);
1022 fArrayCalib->AddAt(objdriftvelocitypad,8);
1024 fArrayCalib->AddAt(objtime0det,9);
1025 fArrayCalib->AddAt(objtime0pad,10);
1027 TGraph *graph = 0x0;
1028 if(FillGraphIndex(&object,graph)){
1029 fGraph->AddAt(graph,3);
1031 TGraph *graphh = 0x0;
1032 if(FillGraphIndex(&objectt,graphh)){
1033 fGraph->AddAt(graphh,4);
1035 }//if(enough statistics)
1037 calibra->ResetVectorFit();
1041 fPostProcess = kTRUE;
1047 //________________________________________________________________________
1048 Bool_t AliTRDcalibration::FillGraphIndex(const TObjArray *vectora,TGraph *graph) const
1051 // Fill one value (the first one) per detector
1054 Int_t loop = (Int_t) vectora->GetEntriesFast();
1056 AliInfo("The Vector Fit is not complete!");
1062 for (Int_t k = 0; k < loop; k++) {
1063 if(!vectora->At(k)){
1068 x[k] = ((AliTRDCalibraFit::AliTRDFitInfo *) vectora->At(k))->GetDetector();
1069 y[k] = ((Float_t *)((AliTRDCalibraFit::AliTRDFitInfo *) vectora->At(k))->GetCoef())[0];
1072 if(!graph) graph = new TGraph(540,&x[0],&y[0]);
1075 new(graph) TGraph(540,&x[0],&y[0]);
1081 //________________________________________________________________________
1082 Bool_t AliTRDcalibration::AddStatsPerDetector(const TH2I *ch)
1085 // Add statistic per detector
1088 AliTRDCalibraMode calibMode = AliTRDCalibraMode();
1089 const char *name = ch->GetTitle();
1090 if(!SetNzFromTObject(name,0,&calibMode)) return 0x0;
1091 if(!SetNrphiFromTObject(name,0,&calibMode)) return 0x0;
1092 if(((calibMode.GetNz(0) == 100) && (calibMode.GetNrphi(0) == 100)) || ((calibMode.GetNz(0) == 10) && (calibMode.GetNrphi(0) == 10))) return kFALSE;
1094 Int_t nybins = ch->GetNbinsY();// groups number
1095 Int_t nxbins = ch->GetNbinsX();// number of bins X
1096 TAxis *xaxis = ch->GetXaxis();
1097 Double_t lowedge = xaxis->GetBinLowEdge(1);
1098 Double_t upedge = xaxis->GetBinUpEdge(nxbins);
1100 // number per chamber 2
1101 calibMode.ModePadCalibration(2,0);
1102 calibMode.ModePadFragmentation(0,2,0,0);
1103 calibMode.SetDetChamb2(0);
1104 Int_t perChamber2 = (Int_t) calibMode.GetDetChamb2(0);
1106 // number per other chamber
1107 calibMode.ModePadCalibration(0,0);
1108 calibMode.ModePadFragmentation(0,0,0,0);
1109 calibMode.SetDetChamb0(0);
1110 Int_t perChamber0 = (Int_t) calibMode.GetDetChamb0(0);
1112 if(nybins != (6*18*perChamber2+6*4*18*perChamber0)) return kFALSE;
1115 TString nname((const char *)ch->GetName());
1116 nname += "PerDetector";
1117 TString title("Nz");
1121 if(!fCHSum) fCHSum = new TH2I((const char *)nname,(const char *)title
1122 ,nxbins,lowedge,upedge,540,0,540);
1125 new(fCHSum) TH2I((const Char_t *) nname,(const char *)title
1126 ,nxbins,lowedge,upedge,540,0,540);
1128 fCHSum->SetYTitle("Detector number");
1129 fCHSum->SetXTitle("charge deposit [a.u]");
1130 fCHSum->SetZTitle("counts");
1131 fCHSum->SetStats(0);
1136 for(Int_t det = 0; det < 540; det++){
1138 Int_t numberofgroup = 0;
1139 if(AliTRDgeometry::GetStack(det) == 2) numberofgroup = perChamber2;
1140 else numberofgroup = perChamber0;
1141 TH1I *projch = (TH1I *) ch->ProjectionX("projch",counter+1,counter+numberofgroup,(Option_t *)"e");
1142 projch->SetDirectory(0);
1144 for(Int_t nx = 0; nx <= nxbins; nx++) {
1145 fCHSum->SetBinContent(nx,det+1,projch->GetBinContent(nx));
1146 fCHSum->SetBinError(nx,det+1,projch->GetBinError(nx));
1149 counter += numberofgroup;
1158 //_____________________________________________________________________________________________________________________
1159 Bool_t AliTRDcalibration::AddStatsPerDetector(const TProfile2D *ph,Int_t i)
1162 // Add statistic per detector
1165 AliTRDCalibraMode calibMode = AliTRDCalibraMode();
1166 const char *name = ph->GetTitle();
1167 //printf("name %s\n",name);
1168 if(!SetNzFromTObject(name,0,&calibMode)) return kFALSE;
1169 if(!SetNrphiFromTObject(name,0,&calibMode)) return kFALSE;
1170 if(((calibMode.GetNz(0) == 100) && (calibMode.GetNrphi(0) == 100)) || ((calibMode.GetNz(0) == 10) && (calibMode.GetNrphi(0) == 10))) return kFALSE;
1171 //printf("Found mode Mz %d, Nrphi %d\n",calibMode.GetNz(0),calibMode.GetNrphi(0));
1174 Int_t nybins = ph->GetNbinsY();// groups number
1175 Int_t nxbins = ph->GetNbinsX();// number of bins X
1176 TAxis *xaxis = ph->GetXaxis();
1177 Double_t lowedge = xaxis->GetBinLowEdge(1);
1178 Double_t upedge = xaxis->GetBinUpEdge(nxbins);
1180 // number per chamber 2
1181 calibMode.ModePadCalibration(2,0);
1182 calibMode.ModePadFragmentation(0,2,0,0);
1183 calibMode.SetDetChamb2(0);
1184 Int_t perChamber2 = (Int_t) calibMode.GetDetChamb2(0);
1186 // number per other chamber
1187 calibMode.ModePadCalibration(0,0);
1188 calibMode.ModePadFragmentation(0,0,0,0);
1189 calibMode.SetDetChamb0(0);
1190 Int_t perChamber0 = (Int_t) calibMode.GetDetChamb0(0);
1192 if(nybins != (6*18*perChamber2+6*4*18*perChamber0)) return kFALSE;
1195 TString nbname((const char *)ph->GetName());
1196 nbname += "PerDetectorVector";
1197 if(!fDetSumVector) fDetSumVector = new AliTRDCalibraVector();
1199 fDetSumVector->~AliTRDCalibraVector();
1200 new(fDetSumVector) AliTRDCalibraVector();
1203 fDetSumVector->SetTimeMax(nxbins);
1206 fDetSumVector->SetNumberBinPRF(nxbins);
1207 fDetSumVector->SetPRFRange(TMath::Abs(lowedge));
1209 fDetSumVector->SetDetCha0(i,1);
1210 fDetSumVector->SetDetCha2(i,1);
1211 fDetSumVector->SetNzNrphi(i,0,0);
1213 Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
1214 fDetSumVector->SetNbGroupPRF(nbg);
1218 TString nname((const char *)ph->GetName());
1219 nname += "PerDetector";
1220 TString title("Nz");
1224 if(!fDetSum) fDetSum = new TH2D((const char *)nname,(const Char_t *) title
1225 ,nxbins,lowedge,upedge,540,0,540);
1228 new(fDetSum) TH2D((const Char_t *) nname,(const Char_t *) title
1229 ,nxbins,lowedge,upedge,540,0,540);
1231 fDetSum->SetYTitle("Detector number");
1232 fDetSum->SetXTitle(xaxis->GetTitle());
1233 fDetSum->SetStats(0);
1237 for(Int_t det = 0; det < 540; det++){
1239 Int_t numberofgroup = 0;
1240 if(AliTRDgeometry::GetStack(det) == 2) numberofgroup = perChamber2;
1241 else numberofgroup = perChamber0;
1243 for(Int_t nx = 1; nx <= nxbins; nx++) {
1245 Double_t entries = 0.0;
1246 Double_t sumw2 = 0.0;
1247 Double_t sumw = 0.0;
1249 for(Int_t k = counter+1; k <= (counter+numberofgroup); k++){
1250 Int_t binnumber = ph->GetBin(nx,k);
1251 entries += ph->GetBinEntries(binnumber);
1252 sumw2 += (ph->GetBinError(binnumber)*ph->GetBinError(binnumber)+ph->GetBinContent(binnumber)*ph->GetBinContent(binnumber))*ph->GetBinEntries(binnumber);
1253 sumw += ph->GetBinContent(binnumber)*ph->GetBinEntries(binnumber);
1256 Double_t mean = 0.0;
1257 if(entries > 0.0) mean = sumw/entries;
1258 Double_t squaremean = 0.0;
1259 if(entries > 0.0) squaremean = sumw2/entries;
1260 Double_t errorf = squaremean - mean*mean;
1261 Double_t error = 0.0;
1262 if(entries > 0.0) error = TMath::Sqrt(TMath::Abs(errorf)/entries);
1264 fDetSum->SetBinContent(nx,det+1,mean);
1265 fDetSum->SetBinError(nx,det+1,error);
1267 if(i==1) fDetSumVector->FillVectorPH(det,0,nx-1,(Int_t)entries,(Float_t)mean,(Float_t)squaremean);
1268 if(i==2) fDetSumVector->FillVectorPRF(det,0,nx-1,(Int_t)entries,(Float_t)mean,(Float_t)squaremean);
1272 counter += numberofgroup;
1280 //_____________________________________________________________________________
1281 Bool_t AliTRDcalibration::SetNrphiFromTObject(const char *name, Int_t i, AliTRDCalibraMode *calibMode) const
1284 // Set the granularity from object
1287 const Char_t *patternrphi0 = "Nrphi0";
1288 const Char_t *patternrphi1 = "Nrphi1";
1289 const Char_t *patternrphi2 = "Nrphi2";
1290 const Char_t *patternrphi3 = "Nrphi3";
1291 const Char_t *patternrphi4 = "Nrphi4";
1292 const Char_t *patternrphi5 = "Nrphi5";
1293 const Char_t *patternrphi6 = "Nrphi6";
1296 const Char_t *patternrphi10 = "Nrphi10";
1297 const Char_t *patternrphi100 = "Nrphi100";
1298 const Char_t *patternz10 = "Nz10";
1299 const Char_t *patternz100 = "Nz100";
1302 if ((strstr(name,patternrphi100)) && (strstr(name,patternz100))) {
1303 calibMode->SetAllTogether(i);
1306 if ((strstr(name,patternrphi10)) && (strstr(name,patternz10))) {
1307 calibMode->SetPerSuperModule(i);
1311 if (strstr(name,patternrphi0)) {
1312 calibMode->SetNrphi(i ,0);
1315 if (strstr(name,patternrphi1)) {
1316 calibMode->SetNrphi(i, 1);
1319 if (strstr(name,patternrphi2)) {
1320 calibMode->SetNrphi(i, 2);
1323 if (strstr(name,patternrphi3)) {
1324 calibMode->SetNrphi(i, 3);
1327 if (strstr(name,patternrphi4)) {
1328 calibMode->SetNrphi(i, 4);
1331 if (strstr(name,patternrphi5)) {
1332 calibMode->SetNrphi(i, 5);
1335 if (strstr(name,patternrphi6)) {
1336 calibMode->SetNrphi(i, 6);
1340 calibMode->SetNrphi(i ,0);
1344 //_____________________________________________________________________________
1345 Bool_t AliTRDcalibration::SetNzFromTObject(const char *name, Int_t i, AliTRDCalibraMode *calibMode) const
1348 // Set fNz[i] of the AliTRDCalibraFit::Instance()
1349 // corresponding to the given TObject
1353 const Char_t *patternz0 = "Nz0";
1354 const Char_t *patternz1 = "Nz1";
1355 const Char_t *patternz2 = "Nz2";
1356 const Char_t *patternz3 = "Nz3";
1357 const Char_t *patternz4 = "Nz4";
1359 const Char_t *patternrphi10 = "Nrphi10";
1360 const Char_t *patternrphi100 = "Nrphi100";
1361 const Char_t *patternz10 = "Nz10";
1362 const Char_t *patternz100 = "Nz100";
1364 if ((strstr(name,patternrphi100)) && (strstr(name,patternz100))) {
1365 calibMode->SetAllTogether(i);
1368 if ((strstr(name,patternrphi10)) && (strstr(name,patternz10))) {
1369 calibMode->SetPerSuperModule(i);
1372 if (strstr(name,patternz0)) {
1373 calibMode->SetNz(i, 0);
1376 if (strstr(name,patternz1)) {
1377 calibMode->SetNz(i ,1);
1380 if (strstr(name,patternz2)) {
1381 calibMode->SetNz(i ,2);
1384 if (strstr(name,patternz3)) {
1385 calibMode->SetNz(i ,3);
1388 if (strstr(name,patternz4)) {
1389 calibMode->SetNz(i ,4);
1393 calibMode->SetNz(i ,0);
1396 //____________________________________________________________________________________________________
1397 Int_t AliTRDcalibration::GetNumberOfGroupsPRF(const char* nametitle) const
1400 // Get numberofgroupsprf
1404 const Char_t *pattern0 = "Ngp0";
1405 const Char_t *pattern1 = "Ngp1";
1406 const Char_t *pattern2 = "Ngp2";
1407 const Char_t *pattern3 = "Ngp3";
1408 const Char_t *pattern4 = "Ngp4";
1409 const Char_t *pattern5 = "Ngp5";
1410 const Char_t *pattern6 = "Ngp6";
1413 if (strstr(nametitle,pattern0)) {
1416 if (strstr(nametitle,pattern1)) {
1419 if (strstr(nametitle,pattern2)) {
1422 if (strstr(nametitle,pattern3)) {
1425 if (strstr(nametitle,pattern4)) {
1428 if (strstr(nametitle,pattern5)) {
1431 if (strstr(nametitle,pattern6)){